Anthropogenic landscape change involves habitat loss and shifts in configuration, often increasing patch number and isolation. Ecological specialization is hypothesized to mediate species vulnerability, with specialists sensitive to habitat amount and generalists potentially utilizing structural heterogeneity. We tested this scale-dependent differential response across seven Pampean coastal dune sites in Argentina, focusing on two sympatric lizards: the endangered habitat specialist Liolaemus multimaculatus and the habitat generalist Liolaemus wiegmannii. We assessed relative lizard abundance using fixed-route visual encounter transects and modeled species abundance association with landscape metrics using Generalized Linear Mixed Models across three nested spatial scales (100-m, 300-m, and 500-m). Results revealed a clear dichotomy driven by specialization. The specialist L. multimaculatus's abundance was explained exclusively by habitat composition metrics, showing a strong positive association with the total area of active dune habitat (500-m scale) and a consistent negative relationship with the total area of low-quality habitats (exotic forestation, semi-fixed dune, and beach) across all three scales. This indicates its vulnerability stems primarily from habitat loss. In contrast, the generalist L. wiegmannii's abundance was explained exclusively by configuration metrics at the smallest scale (100-m), showing a positive association with the number and complexity of semi-fixed dune patches. This suggests the generalist benefits from fragmentation-induced heterogeneity. Our findings confirm that ecological specialization critically dictates species vulnerability and scale-dependent responses to landscape change. Conservation efforts for the vulnerable L. multimaculatus must prioritize maintaining the amount of core active dune habitat.
Anthropized landscapes are complex mosaics whose composition and configuration are widely recognized as critical factors influencing biodiversity (Duflot et al., 2017; Martin et al., 2019). While landscape composition is determined by the amount of different land cover (Li and Reynolds, 1995; Bennett et al., 2006), configuration refers to the spatial arrangement, number, and isolation of the patches (Lindenmayer and Hobbs, 2007). Land-use change, a primary driver of global biodiversity loss, alters both the landscape’s composition and configuration at the landscape level (Bogaert et al., 2011; Ramirez-Arce et al., 2022). This transformation often simultaneously reduces available habitat area, increases patch isolation, and modifies functional connectivity, making it difficult to assess which landscape attribute drives a species' response (Mimet et al., 2016). Crucially, habitat fragmentation per se refers to changes in spatial configuration (e.g., patch number or isolation) that occur independent of any change in total amount of different land cover (Fahrig, 2017). Although the detrimental effects of habitat loss are broadly accepted, the effect of fragmentation per se remains debated due to differences in conceptual frameworks and analytical models (Arroyo-Rodríguez et al., 2016; Tan et al., 2023; dos Anjos et al., 2025). Therefore, effective ecological conservation requires a multi-scale understanding of how animal populations respond to the independent effects of habitat loss (composition metrics) and fragmentation (configuration metrics) induced by land-use change (McGarigal, 1995; Mimet et al., 2016; Galán-Acedo et al., 2024).
Two ecological conceptual frameworks posit different mechanisms driving species dynamics in these landscapes: the Habitat Amount Hypothesis contends that species abundance and richness decline proportionally with local habitat loss (Fahrig, 2013); conversely, the effect of fragment size and isolation may determine dynamics through demographic effects that cause a disproportionate reduction (Westaway et al., 2024). Small vertebrates exhibit differential vulnerability to landscape change, a phenomenon modulated by life history traits and ecological specializations (Cruz-Salazar et al., 2016; Amburgey et al., 2021). Species with broad niches (generalists) demonstrate less spatial constraint than specialists, which are highly sensitive to local habitat features and can perceive the surrounding matrix as a dispersal barrier (Swihart et al., 2003). Consequently, generalists can often utilize non-matrix habitat or even view the matrix as habitat, suggesting that matrix quality can mediate species responses (Devictor et al., 2008; Fahrig, 2017). In contrast, specialist species exhibit strong responses to habitat amount and patch features, and greater isolation in fragmented landscapes (D’Eon et al., 2002). This variability necessitates multi-scale analyses, as species responses differ across local, patch-level, and broader landscape configurations (da Silva Carneiro et al., 2022).
For reptiles, species abundance is often more strongly influenced by landscape composition than configuration (Jellinek et al., 2004; McAlpine et al., 2015; Ryberg and Fitzgerald, 2016). While disaggregation of suitable patches and patch shape irregularity may have negative effects on small-bodied, low-vagility reptiles, the extent of suitable patch cover is generally the strongest driver of their distribution (Ghosh and Basu, 2020; Prieto-Ramírez, 2023; Díaz et al., 2024). However, habitat specialization exacerbates vulnerability; specialists experience disproportionately adverse effects compared to generalists (Keinath et al., 2017), as generalists capitalize on edge habitats and matrix resources, dominating fragmented landscapes, while specialists are confined to high-quality patches (Russildi et al., 2016). Some reptiles also respond to configuration metrics like patch shape complexity (e.g., perimeter-to-area ratio), which may advantage generalists requiring resources in multiple land covers (Nimmo et al., 2019; Mulhall et al., 2022).
In the Pampas coastal region of Argentina, anthropogenic activities have significantly altered the composition and configuration of the Eastern Dune Barrier landscape, with human impact increasing dramatically over the last three decades (Austrich et al., 2021). This habitat modification threatens two sympatric, sand-dwelling lizard species with contrasting life history strategies: Liolaemus multimaculatus, an endangered habitat specialist restricted to open sandy grasslands (Kacoliris et al., 2016; Vega et al., 2019), and the habitat generalist Liolaemus wiegmannii, which inhabits from xerophilic forests to dune grasslands with low-to-moderate anthropogenic impact (Ramirez Pinilla, 1991; Martori et al., 1998; Stellatelli et al., 2015a, b). Previous studies have shown that while the presence of L. multimaculatus is positively associated with the size of active dune patches, and negatively associated with fragmentation (Block et al., 2023), the presence of L. wiegmannii is negatively associated with the aggregation and complexity of semi-fixed dune patches and with the coverage of exotic forestation (Block et al., 2016).
Despite this body of work, the influence of multi-scale landscape fragmentation on the relative abundance of Liolaemus lizards remains a critical knowledge gap. Specifically, how individual species' responses vary across configuration and compositional attributes depending on the spatial scale warrants investigation. The primary objective of this study was to assess the effects of habitat composition and configuration on the abundance patterns of L. multimaculatus and L. wiegmannii across the Southern Pampean Coast and to identify the key landscape metrics driving these spatial dynamics at small (100-m), medium (300-m), and large (500-m) spatial scales. We examined the associations between lizard abundance and metrics of habitat composition (amount of specific habitat types) and configuration (mean patch size, number of patches, edge longitude, and patch shape complexity). Based on the principles of the species-specific response framework, we formulated two primary, non-mutually exclusive hypotheses. Hypothesis 1 (Habitat Amount): We hypothesized that ecological specialization dictates a differential response, with the specialist L. multimaculatus being primarily affected by the total available area of its core habitat, while the generalist L. wiegmannii tolerates variability in habitat amount. Prediction 1: The abundance of L. multimaculatus will show a strong, positive association with the amount (coverage) of suitable active dune habitat across multiple spatial scales, and changes in habitat amount will explain a greater proportion of its abundance variation than configuration metrics, providing support for the Habitat Amount Hypothesis. Conversely, L. wiegmannii abundance is predicted to be less sensitive to the total coverage of its preferred semi-fixed dune habitat. Hypothesis 2 (Fragmentation per se): Landscape configuration effects, or fragmentation per se, will be contingent upon the species' specialization, with the narrow-niched specialist being negatively impacted by habitat pattern due to increased barrier effects, and the generalist benefiting from the structural heterogeneity characteristic of fragmentation per se. Prediction 2: L. multimaculatus abundance will be negatively associated with configuration metrics indicative of fragmentation (e.g., high patch density, high edge longitude, small patches), reflecting a strong barrier effect from the matrix and reduced functional connectivity. In contrast, L. wiegmannii abundance will show a null or positive association with metrics of patch complexity and number (fragmentation), reflecting its inherent flexibility and ability to capitalize on landscape heterogeneity and potential for resource complementation.
Materials and methodsStudy areaThe study was conducted within the temperate coastal Pampa dune landscape of Buenos Aires province, Argentina, focusing specifically on the Eastern Dune Barrier (Fig. 1), which extends for approximately 192 linear km (Bértola et al., 2021). Sampling sites were established across a longitudinal extent of approximately 130 km, encompassing representative sectors of this dune system, including, from north to south: 1-Punta Médanos (36°53.00'30.96"S, 56°40.00'50.90"W); 2-Mar de las Pampas (37°19.00'21.01"S, 57°2.00' 44.58"W); 3-Mar Azul (37°20.00'37.31"S, 57°3.00'37.84"W); 4-Faro Querandí Nature Reserve (37°27.00'11.38"S, 57°6.00'43.22"W); 5-Parque Atlántico Mar Chiquito UNESCO–MAB Biosphere Reserve (37°44.00'16.12"S, 57°24.00'53.16"W); 6-Balneario Parque Mar Chiquito (37°45.00'36.20"S, 57°26.00'4.26"W); 7-La Caleta (37°47.00'26.42"S, 57°27.00'40.90"W; Fig. 1). Liolaemus multimaculatus was sampled at sites 1, 4, 5, 6, and 7; while L. wiegmannii was sampled at sites 1, 2, 3, and 4. The dune system exhibits a typical geomorphological zoning pattern from sea to upland, encompassing upper beach, active fore-dunes, active inland dunes, interdune depressions, and semi-fixed dunes (Wiedemann and Pickart, 2004). This region experiences a temperate oceanic climate, characterized by humid weather, an average maximum summer temperature of 33 °C, an average winter temperature of 9 °C, and mean annual precipitation ranging from 800 to 1000 mm (Aliaga et al., 2017; Burgos and Vidal, 1995; NASA, 2000–2016). It is noteworthy that these coastal dunes represent one of the last remnants of native Pampean grasslands, though significant structural alteration, including the replacement of native plants with exotic trees associated with human settlements, has impacted the ecosystem (Alberio and Comparatore, 2014; Yezzi et al., 2021).
Location of sampled sites of Liolaemus multimaculatus and Liolaemus wiegmannii in Pampas coastal dunes of Argentina. The grey area in the left map demarcates the Pampas Region in South America. The detailed map at the right shows the surveyed sites in the Eastern Dune Barrier. Numbers 1–7 indicate sampling sites. References: 1—Punta Médanos; 2—Mar de las Pampas; 3—Mar Azul; 4— Faro Querandí Nature Reserve; 5— Parque Atlántico Mar Chiquito UNESCO–MAB Biosphere Reserve; 6—Balneario Parque Mar Chiquito; 7—La Caleta.
Field surveys for the two focal species, L. multimaculatus and L. wiegmannii, were conducted over three consecutive activity seasons: January-March 2021, October 2021-March 2022, and October 2022-March 2023, covering the peak of the lizards' annual activity in the temperate Pampas (Vega et al., 2019). We employed the line transect methodology (100 m long, 6 m wide), a standardized technique recommended for comparing abundances across broad geographical scales (Smolensky and Fitzgerald, 2011). A total of 119 fixed-route transects were established across seven study sites, 60 dedicated to L. multimaculatus and 59 to L. wiegmannii. Transects were distributed randomly, with the number of transects per site proportional to its area (see Supplementary Table S1) and were separated by a minimum distance of 100 m to ensure spatial independence, considering the mean home ranges of L. multimaculatus (mean = 219 m2) and L. wiegmannii (mean = 39 m2) (Stellatelli et al., 2016a). Lizard counts were performed monthly on each transect during the hours of maximal daily activity (09:00 to 18:00 h) (Block et al., 2016), resulting in 12–15 surveys per transect and a cumulative effort of 1,608 temporal transect-replicates (807 for L. multimaculatus and 801 for L. wiegmannii) over the three-year period. All surveys were conducted by the same researcher (J.E. Dajil) and restricted to periods with adequate weather conditions for lizard activity, characterized by sunny days, low winds, and suitable environmental air temperatures (around 22 °C; Block et al., 2013). To minimize detection bias, the search time was adjusted proportionally to the vegetation cover, and only lizards observed within 3 m of the observer were recorded, consistent with previous work demonstrating uniform detection rates within this radius across coastal dune vegetation (Stellatelli et al., 2015b). To formally address the inherent challenges of imperfect detection—where failure to detect a species does not imply its absence (Weir et al., 2005)—we utilized site occupancy models to simultaneously estimate the probability of detection and true occupancy (MacKenzie et al., 2002; Mazerolle et al., 2007).
Satellite image processing and landscape variablesA Geographic Information System (GIS) was developed using the open-source software QGIS (version 3.14.0-Pi) to characterize the coastal dune landscape. Habitat classification was performed via supervised classification on a Landsat 8 Operational Land Imager satellite image (30 m resolution) acquired in February 2022 (Path: 223; Row: 086; UTM/WGS-84; available from USGS Earth Explorer). The classification utilized the Semi-Automatic Classification Plugin (SCP), an open-source QGIS tool, employing the Spectral Angle Mapper algorithm (Congedo, 2021). Urban areas, roads, and the beach were pre-digitized as vector shapefiles, rasterized, and subsequently integrated into the final classification. The resulting thematic map (Fig. 1) delineated eight land cover categories: (1) Active Dunes (sparsely vegetated upper beach, active foredunes, and active inland dunes); (2) Semi-Fixed Dunes (continuous natural grassland, clump herbs, and shrubs); (3) Interdune Depressions (flood-associated vegetation); (4) Forested Dunes (exotic trees); (5) Urban Areas; (6) Beaches; (7) Water Bodies; and (8) Roads. Classification accuracy was assessed using an error matrix analysis and the Kappa index (Campbell, 2002), yielding an overall accuracy of 91.71% and a Kappa statistic of 0.88. The decision to use a single, cloud-free image from February 2022 for all landscape analyses (spanning 2021–2023 surveys) was a methodological choice made to prioritize internal consistency and minimize classification error. A preliminary visual assessment of the study area (2021–2023) using Google Earth Engine suggested that major landscape alterations due to new urbanization or large-scale afforestation were minimal and likely imperceptible during the sampling period. By focusing classification efforts on the single, best-quality image from the central year of the temporally narrow sampling period, we ensured maximal spectral consistency across all survey transects, thereby prioritizing the comparability of resulting landscape metrics over the capture of minor, short-term temporal dynamics. This approach provided the most robust and consistent landscape dataset for characterizing the habitat attributes relevant to this study.
Landscape metrics were obtained from the thematic map of the Eastern Dune Barrier using FRAGSTATS software (v4.0; McGarigal et al., 2012). We employed a multi-scale approach to ensure that predictor variables were measured at scales commensurate with the organism's known ecological processes, as recommended by Elith and Leathwick (2009). We measured the landscape metrics within circular windows centered on each transect, creating a nested series of three spatial scales with radii of 100 m, 300 m, and 500 m. This selection of scales was explicitly tied to the known spatial ecology of the target species. The smallest buffer, 100 m, was substantially larger than the reported maximum linear movement distance for both species (L. multimaculatus mean linear distance: 3–28 m; L. wiegmannii mean linear distance: 0.7–10 m) and fully encompasses their reported home range variation (L. multimaculatus: 13 to 570 m2; L. wiegmannii: 6 to 79 m2; Stellatelli et al., 2016a, b). Furthermore, spatial resolution of the satellite image (30 × 30 m) included the mean size of the home range estimated for both species (Block et al., 2016, 2023). The two larger scales, 300 m and 500 m, were chosen to capture the potential influence of landscape structure on broader spatial extensions, representing ecologically relevant scales beyond immediate daily activities that might affect colonization, gene flow, and long-term population persistence. To reflect species-specific microhabitat preferences, composition metrics for L. multimaculatus were calculated focusing on active dunes, while for L. wiegmannii, semi-fixed dunes were considered (Vega et al., 2019). From each landscape unit (100, 300, and 500 m radius), composition metrics included AFD, total area of forested dunes; AAD, total area of active dunes; ASD, total area of semi-fixed dunes; AID, total area of interdune depressions; ABE, total area of the beach. Configuration metrics included TAE, total length of the active dune edge; TSE, total length of the semi-fixed dune edge; SAD, geometric complexity of active dune patches; SSD; geometric complexity of semi-fixed dune patches; NPA, number of patches of active dune; NPS, number of patches of semi-fixed dunes; XAD, mean patch size of active dunes; XSD, mean patch size of semi-fixed dunes; VAD, coefficient of variation of the area of active dune patches; VSD, coefficient of variation of the area of semi-fixed dune patches; dUR, distance to urbanization; dRO, distance to road (see Supplementary Table S2).
Statistical analysisTo account for imperfect detection and provide a measure of detection probability, we implemented occupancy models (MacKenzie et al., 2002) using the “occu” function in the R package “unmarked” (Fiske and Chandler, 2011). Abundance data were converted to binary detection histories (1 = one or more individuals detected; 0 = no individuals detected) for each of the 12–15 repeated visits per transect. These models simultaneously yielded estimates for two fundamental parameters: (1) the true site occupancy probability (ψ), defined as the proportion of transects occupied by the species, and (2) the detection probability (P), defined as the probability of detecting the species during a single visit conditional on presence. The estimated single-visit detection probability was subsequently used to calculate the cumulative detection probability across the entire sampling period (Kellner et al., 2023).
We assessed the associations between lizard relative abundance and landscape metrics using Generalized Linear Mixed Models (GLMMs) with a Poisson error structure and a log-link function (Crawley, 2012). Preliminary analysis, comparing zero-inflated Poisson models (using “glmmTMB”) with standard Poisson models, indicated that the more parsimonious non-zero-inflated models provided a better fit to our data and were thus retained. All models incorporated random factors to account for the hierarchical structure of the sampling design: transect identity nested within sites (Bolker et al., 2009). To address the temporal autocorrelation arising from repeated measures, all models included an autoregressive structure of order 1 [ar1(0 + as.factor(time) | ID)], where “ID” identifies the individual transect and “time” represents the temporal order of the observation. This structure explicitly models the temporal dependence, assuming that the correlation between observations decreases exponentially with time, which is a suitable assumption given our monthly sampling intervals. The term "0 + as.factor(time)" allows for the estimation of a single autocorrelation parameter (rho) across all transects, while the nesting via ID ensures the temporal structure is applied to each transect (Bolker et al., 2009).The landscape metrics, obtained at radii of 100, 300, and 500 m for each transect, were included as fixed effects (Supplementary Table S4). Before model construction, multicollinearity among predictors was initially explored using the Spearman rank correlation test (Zar, 1999). We constructed global models combining non-correlated landscape metrics (Spearman's r < 0.5 and/or p ≥ 0.05; Supplementary Tables S5a-f) and ensuring a Variance Inflation Factor (VIF) of less than 5 (Supplementary Table S4B). This procedure yielded a total of 110 candidate global models across all species and scales: L. multimaculatus: 22 models at 500 m, 18 at 300 m, and 6 at 100 m; L. wiegmannii: 19 models at 500 m, 25 at 300 m, and 20 at 100 m. The model with the lowest Akaike Information Criterion corrected for small sample sizes (AICc) value was selected for each species-scale combination. This resulted in a total of six retained global models (one per species and spatial scale) for subsequent analysis (Supplementary Table S4). Selection of all sub-models (combinations of explanatory variables within the selected global models) was based on the AICc, as it is designed to optimize predictive capacity, which aligns with our focus on ecological prediction and spatial pattern identification (Vrieze, 2012; Aho et al., 2014). We reported model uncertainty using ΔAICc (where a value <2 indicates similar performance) and Akaike weights (the probability that a given model is the best approximating model) using the “MuMIn” package (Barton, 2022). Multi-model inference was performed via AICc weight-based model averaging to calculate final parameter estimates and their 95% Confidence Intervals (CIs). Overlap of the CI with zero was used to assess the strength of support for each predictor (Burnham and Anderson, 2002). All analyses were conducted in R (R Core Team, 2023), with models run using the “glmmTMB” package (Brooks et al., 2017) and model assumptions rigorously checked using the “DHARMa” package (Hartig, 2016; Supplementary Fig. S1).
ResultsThe mean occupancy probability (ψ) for L. multimaculatus was 0.84, with a cumulative detection probability of 0.98 and a single-visit detection probability (P) of 0.25 (95% confidence interval [CI]: 0.22–0.29). Spatial variation in the relative abundance of L. multimaculatus was evident, with mean transect abundance ranging from 0.08 (95% CI: 0.05–0.13) individuals per transect at Punta Médanos to the highest values of 0.28 (95% CI: 0.17–0.45) individuals per transect at Parque Atlántico Mar Chiquito UNESCO–MAB Biosphere Reserve and 0.28 (95% CI: 0.21–0.38) individuals per transect at Faro Querandí Nature Reserve (Fig. 2). Abundance variation was exclusively associated with habitat composition metrics across all three analyzed spatial scales. At the smallest landscape scale (100 m radius), two models explained L. multimaculatus abundance (Supplementary Table S3), both showing a significant negative association with the total area of semi-fixed dunes (ASD) and beach (ABE) (Table 1; Fig. 3). At the intermediate scale (300 m radius), six models were significant (Supplementary Table S3), with total area of forested dunes (AFD) showing a significant negative association with abundance (Table 1; Fig. 3). At the largest scale (500 m radius), two models were significant (Supplementary Table S3), and abundance was positively associated with the total area of active dunes (AAD) and negatively associated with ASD, AFD, and ABE (Table 1; Fig. 3).
Estimated relative abundance per transect for the specialist Liolaemus multimaculatus and the generalist Liolaemus wiegmannii across seven sampling sites in the Pampas coastal dunes of Argentina. Circles represent the mean estimated abundance, and vertical bars indicate the 95% confidence intervals. The numbers above the bars indicate the total number of transect-replicates. Numbers 1–7 indicate sampling sites. References: 1—Punta Médanos; 2—Mar de las Pampas; 3—Mar Azul; 4— Faro Querandí Nature Reserve; 5— Parque Atlántico Mar Chiquito UNESCO–MAB Biosphere Reserve; 6—Balneario Parque Mar Chiquito; 7—La Caleta.
Coefficient estimates (±SE) from generalized linear mixed models (GLMM) describing the landscape metrics factors affecting the relative abundances of Liolaemus multimaculatus and Liolaemus wiegmannii at different scales (landscape units of 100, 300 and 500 m radius) in the Pampas coastal dunes of Argentina. Parameter estimates are weighted averages (using AICc weights) from all models, and standard errors (SE) were calculated from all candidate models from unconditional variances. References: AFD, total area of forested dune; ASD, total area of semi-fixed dune; dRO, distance to road; ABE, total area of beach; SAD, geometric complexity of active dune patches; AAD, total area of active dune; TAE, total length of active dune edge; SSD, geometric complexity of semi-fixed dune patches; NPS, number of semi-fixed dune; PAS, percentage of like adjacencies of semi-fixed dune; TSE, total length of semi-fixed dune edge. Explanatory variables with a confidence interval (CI) excluding zero are in bold.
| Scale | L. multimaculatus | L. wiegmannii | ||||
|---|---|---|---|---|---|---|
| Explanatory variables | Coefficient ± SE | CI | Explanatory variables | Coefficient ± SE | CI | |
| 100 m | Intercept | −1.61 ± 0.20 | [−2.00, −1.22] | Intercept | −1.10 ± 0.08 | [−1.26, −0.94] |
| AFD | −0.36 ± 0.20 | [−0.75, 0.04] | SSD | 0.25 ± 0.11 | [0.03, 0.47] | |
| ASD | −0.3 ± 0.13 | [−0.54, −0.05] | NPS | 0.34 ± 0.10 | [0.13, 0.55] | |
| dRO | 0.28 ± 0.17 | [−0.05, 0.62] | AAD | 0.11 ± 0.08 | [−0.05, 0.28] | |
| ABE | −0.39 ± 0.15 | [−0.68, −0.09] | PAS | 0.17 ± 0.11 | [−0.05, 0.48] | |
| TAE | 0.05 ± 0.12 | [−0.20, 0.30] | AFD | −0.05 ± 0.08 | [−0.20, 0.10] | |
| 300 m | Intercept | −1.77 ± 0.22 | [−2.20, −1.34] | Intercept | −1.49 ± 0.11 | [−1.73, −1.26] |
| AFD | −0.83 ± 0.29 | [−1.40, −0.26] | TSE | 0.09 ± 0.08 | [−0.06, 0.25] | |
| dRO | 0.30 ± 0.17 | [−0.04, 0.64] | AFD | −0.08 ± 0.08 | [−0.24, 0.08] | |
| TAE | −0.14 ± 0.12 | [-0.37, 0.10] | ABE | 0.03 ± 0.07 | [−0.11, 0.18] | |
| SAD | 0.13 ± 0.11 | [-0.08, 0.35] | AAD | 0.01 ± 0.08 | [−0.14, 0.17] | |
| 500 m | Intercept | −1.79 ± 0.15 | [−2.08, −1.50] | Intercept | −1.50 ± 0.11 | [−1.73, −1.26] |
| AFD | −0.83 ± 0.23 | [−1.28, −0.38] | SSD | 0.14 ± 0.07 | [−0.01, 0.29] | |
| ASD | −0.26 ± 0.12 | [−0.50, −0.03] | TSE | 0.11 ± 0.08 | [−0.04, 0.27] | |
| ABE | −0.31 ± 0.11 | [−0.52, −0.10] | ABE | 0.09 ± 0.08 | [−0.05, 0.25 | |
| AAD | 0.21 ± 0.11 | [0.01, 0.42] | AFD | −0.06 ± 0.08 | [−0.23, 0.10] | |
Associations between the relative abundance of Liolaemus multimaculatus and Liolaemus wiegmannii and significant landscape metrics from generalized linear mixed models (GLMMs). (A) Associations between the relative abundance of L. multimaculatus and the total area (ha) of forested dunes (AFD), semi-fixed dunes (ASD), beach (ABE), and active dunes (AAD) within 500 m radius landscape units (LU500). (B) Associations between the relative abundance of L. multimaculatus and the total area (ha) of forested dunes (AFD) within 300 m radius landscape units (LU300). (C) Associations between the relative abundance of L. multimaculatus and the total area (ha) of semi-fixed dunes (ASD) and beach (ABE) within 100 m radius landscape units (LU100). (D) Associations between the relative abundance of L. wiegmannii and the number of patches (NPS) and patch shape complexity of semi-fixed dunes (SSD) within 100 m radius landscape units (LU100).
For L. wiegmannii, the mean value of ψ was 1.00, with a cumulative detection probability of 0.99 and a P value of 0.27 (95% CI: 0.24–0.30). The mean relative abundance of L. wiegmannii also varied spatially, ranging from 0.20 (95% CI: 0.13–0.30) individuals/transect at Mar de las Pampas to a maximum of 0.33 (95% CI: 0.22–0.50) individuals/transect at Punta Médanos (Fig. 2). This variation in abundance was explained exclusively by metrics that described habitat configuration, only at the smallest scale (100 m radius), where seven models were significant (Supplementary Table S3). These models showed a positive association with both semi-fixed dune patch complexity (SSD) and number of patches (NPS) (Table 1; Fig. 3). At intermediate and large scales, four and nine models, respectively, explained the relative abundance of L. wiegmannii (Supplementary Table S3), but none of the variables within these models were statistically significant (Table 1).
DiscussionOur findings using the co-occurrence of the habitat specialist L. multimaculatus and the habitat generalist L. wiegmannii in the coastal Pampas dunes strongly support the hypothesis that ecological specialization critically mediates species' vulnerability and response to landscape change in anthropogenically modified landscapes. The occupancy analysis revealed a clear difference in estimated occupancy probability between these two sympatric sand-dwelling lizards. The specialist, L. multimaculatus, exhibited a lower mean occupancy probability (ψ = 0.84) compared to its generalist congener (ψ = 1.00), consistent with its high association with the specific, restricted, and threatened active dune habitat (Block et al., 2023). Although the single-visit detection probability (P = 0.25) of L. multimaculatus was slightly lower than that of L. wiegmannii (P = 0.27), possibly due to behavioral traits like cryptic coloration or a tendency to bury (Stellatelli et al., 2015a), the high cumulative detection probability (0.98) confirms that the lower occupancy estimate accurately reflects its less widespread distribution, rather than a sampling artifact. Conversely, the habitat generalist, L. wiegmannii, showed a near-certain mean occupancy (ψ = 1.00), indicating its ubiquitous presence in virtually all surveyed suitable patches, which aligns with its description as a species with a relatively wide distribution across sandy substrates of multiple biomes (Abdala et al., 2021). This clear difference in occupancy estimates underscores the value of occupancy modeling, which accounts for imperfect detection, a common issue for cryptic or mobile species like these lizards (Sewell et al., 2012; Turner et al., 2023), suggesting contrasting responses to landscape attributes due to natural or anthropogenic causes.
Consistent with our first hypothesis, landscape composition metrics accounted for all observed variation in L. multimaculatus abundance across all three analyzed scales. The species exhibited a strong positive association with the total area of active dunes—the specialized habitat type—at the largest scale (500-m radius) and a significant negative association with the total area of other habitat categories (semi-fixed dunes, beach, and exotic forestations) at all three scales. This scale-dependent pattern is in alignment with regional findings in the southern coastal Pampas, where patches retaining natural dune structure support local species abundance (Stellatelli et al., 2015b) and lizard occurrence positively correlates with large areas of active dunes (Block et al., 2023). These results confirm our prediction that the abundance of the endemic specialist L. multimaculatus would be significantly associated with the amount of suitable active dune habitat, and that changes in habitat amount would explain a greater proportion of its abundance variation than configuration metrics, providing strong support for the Habitat Amount Hypothesis (Fahrig, 2013). Our findings echo broader ecological knowledge, indicating that specialist reptiles are negatively impacted by native vegetation loss and encroachment of exotic plants, which disrupt food availability and thermal habitat quality (Valentine et al., 2007; Ghosh and Basu, 2020). The strong negative association of this species with land-use change resulting compositionally altered patches emphasizes that these adjacent land covers—especially exotic forestation—function as a dispersal barrier, strongly confining the specialist to its core suitable habitat (Young et al., 2018; Wenner et al., 2022), thereby highlighting the critical importance of natural habitat structure at multiple scales (Mulhall et al., 2022).
In stark contrast, the abundance of the generalist L. wiegmannii was explained exclusively by configuration metrics, and only at the small scale (100-m radius). Specifically, L. wiegmannii abundance was positively associated with patch shape complexity and the total number of semi-fixed dune patches. This partially supports our second hypothesis regarding the generalist, which predicted a positive association with fragmentation metrics, reflecting its greater tolerance of disturbance and its ability to exploit structurally complex, heterogeneous landscapes (Ramiadantsoa et al., 2018). This finding is congruent with Block et al. (2016), who reported that the probability of L. wiegmannii occurrence increased in disaggregated patches of semi-fixed dunes. This positive association with local-scale fragmentation is consistent with the generalist's inherent flexibility, allowing it to capitalize on complex patch edges, which likely facilitate dispersal and resource exploitation between nearby patches, a form of landscape complementation (Fahrig, 2017). This finding aligns with observations of generalist herpetofauna benefiting from more geometrically complex patches (Clements et al., 2024), suggesting that, for this species, edge habitats and the surrounding matrix are utilized rather than acting as a strong barrier (Devictor et al., 2008).
However, these results should be interpreted with caution, as the absence of a significant response in L. wiegmannii abundance to landscape composition or configuration at intermediate or broad scales (300-m and 500-m radii) warrants discussion. Landscape structure and function are scale-dependent, and, therefore, landscape metrics also depend on the scale at which they are measured (Turner, 1989). The positive relationship between lizard abundance and total number of semi-fixed dune patches at small spatial extents may be associated with a bias between the size of the selected landscape unit and landscape variability (Block et al., 2016). For instance, some metrics, such as patch density, tend to decrease when the extent of the landscape unit increases because the number of patches is higher near the edges than farther away (Saura and Martínez-Millán, 2001). This could partially explain the lack of a broader-scale effect of landscape fragmentation on L. wiegmannii abundance. While this result contrasts with our second prediction—which anticipated a positive response to fragmentation at broader scales—it supports the hypothesis that the generalist would be less sensitive to the total coverage of its preferred semi-fixed dune habitat. The finding suggests that the factors driving L. wiegmannii abundance are highly localized, a pattern that differs from previous studies which reported small- and intermediate-scale associations between the species' presence and landscape composition and configuration (Block et al., 2016). This discrepancy highlights that different quantification measures (abundance vs. presence) and ecological processes can manifest at different spatial scales (Dibner et al., 2017; Rotem et al., 2020), further emphasizing the importance of fine-scale factors that determine habitat quality, such as vegetation structure and microhabitat features, for generalist reptiles (Holland and Bennett, 2009). Ultimately, the spatial segregation previously reported—L. multimaculatus in active dunes and L. wiegmannii in semi-fixed dunes (Block et al., 2016; Vega et al., 2019; Block et al., 2023)—along with the tolerance of L. wiegmannii to slight changes in patch structure and the thermal environment (Stellatelli et al., 2013, 2015b), likely underlie this observed differential, scale-dependent response, with small-scale heterogeneity being key for the generalist.
These contrasting results—showing composition as the driver for the specialist across all scales and configuration influencing the generalist only at the small scale—underscore that species-specific habitat specialization is the critical modulator of responses to landscape variables, a finding consistent with broader reptile ecology studies (Ghosh and Basu, 2020). Consequently, the conservation priority for sand-specialist lizards, such as L. multimaculatus, must be the protection and expansion of core active dune habitat (Leavitt and Fitzgerald, 2013; Block et al., 2023). Specifically, this necessitates the systematic control of exotic plant expansion, a primary stressor that negatively affects landscape composition and functional connectivity (Schlesinger et al., 2020). To ground the conservation implications in evidence, conservation efforts for the endangered L. multimaculatus must focus on meeting demographic viability requirements. Population viability analysis models of Kacoliris et al. (2019) have determined that populations of >2,400 individuals are required for long-term viability (i.e., the minimum viable population), and thus, active dune restoration and protection should be quantified to ensure the resulting habitat patches are large enough to support this population size, meeting the species' minimum area requirement. Conversely, management for the generalist L. wiegmannii requires a focus on maintaining small-scale structural heterogeneity and the complexity of the natural semi-fixed dune matrix. Therefore, the management of this anthropized coastal ecosystem must adopt a multi-faceted, species-specific approach, which is essential for mitigating the loss of the specialist and preserving the overall biodiversity of the Pampean coastal dunes.
ConclusionIn summary, our multi-scale analysis demonstrates that a species' degree of habitat specialization critically dictates its response to landscape change, providing an essential framework for conservation in the Southern Pampean Coastal Dunes. The habitat specialist L. multimaculatus responded entirely to scale-dependent habitat composition (amount of active dunes), confirming its reliance on core habitat availability and its extreme vulnerability to habitat loss. Conversely, the generalist L. wiegmannii was influenced exclusively by small-scale configuration metrics related to semi-fixed dune patchiness and shape complexity, reflecting its adaptive ability to exploit structural heterogeneity across the fragmented landscape. While this clear dichotomy strongly supports theoretical predictions on how specialization dictates vulnerability (Chandler and Hepinstall-Cymerman, 2016), the observed abundance patterns are also highly susceptible to finer-scale ecological processes and confounding factors. The lizards' low abundance and sparse distribution are likely exacerbated by unmeasured threats, including predation pressure (Stellatelli et al., 2015a), direct anthropogenic disturbance from off-road vehicles (Vega et al., 2000), and intrinsic microclimatic vulnerability due to habitat alterations (Block et al., 2013; Stellatelli et al., 2013, 2015b). Furthermore, although current spatial segregation suggests ecological partitioning (Vega et al., 2019), persistent anthropogenically-driven habitat loss risks intensifying interspecific competition by forcing these sympatric species into shared habitats, potentially leading to competitive exclusion by the more resilient generalist (Romano et al., 2022). This integrated perspective necessitates a multi-faceted management strategy. Prioritizing land protection and habitat restoration for the specialist with evidence-based targets aimed at achieving a minimum viable population, while simultaneously managing small-scale structural heterogeneity for the generalist. The generalist’s positive response to fragmentation per se aligns with global ecological literature, suggesting that configuration, independent of habitat loss, can enhance functional connectivity and landscape complementation. Ultimately, future research employing controlled experiments or fine-scale spatiotemporal modeling is warranted to rigorously disentangle the roles of local factors and species interactions from the dominant landscape-scale drivers.
Declaration of Generative AI and AI-assisted technologies in the writing processDuring the preparation of this work the author used GEMINI in order to check the grammar. After using this tool, the author reviewed and edited the content as needed and takes full responsibility for the content of the publication.
FundingThis work was supported by the Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT)- Fondo para la Investigación Científica y Tecnológica [PICT 2016-0677, PICT 2016-0266]; the Universidad Nacional de Mar del Plata [15/E1093, EXA1133/23]; the CONICET doctoral fellowship; and the Neotropical Grassland Conservancy (NGC) nonprofit corporation.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
We thank the Dirección de Áreas Protegidas PBA (Ministerio de Ambiente, Buenos Aires, Argentina) for granting research permit (N° 003/18).
The following is Supplementary data to this article:
