Environmental, land cover and land use constraints on the distributional patterns of anurans: Leptodacylus species (Anura, Leptodactylidae) from Dry Chaco

Selected species

We selected ten species of the genus Leptodactylus, whose main distribution area is in the Dry Chaco ecoregion: Leptodactylus bufonius, L. chaquensis, L. elenae, L. fuscus, L. gracilis, L. laticeps, L. latinasus, L. latrans, L. mystacinus and L. podicipinus. Among these species, Leptodactylus laticeps is the only one considered as near threatened by IUCN assessments (Cortez et al., 2004); and vulnerable by the categorization of Argentinean amphibians (Schaefer & Céspedez, 2012).

Occurrence data

Presence points databases were constructed based on revised specimens of collections, databases of herpetological collections, bibliography and free databases. The material revised proceeded from the following herpetological collections: Argentina: Fundación Miguel Lillo (FML), Museo Argentino de Ciencias Naturales (MACN), Museo de la Plata (MLP), private collections of Dr. María Laura Ponssa (L) and Dr. Julián Lescano (JL); Brazil: Museu Nacional de Rio de Janeiro (MNRJ), Museu de Zoologia de São Paulo (MZUSP), and Paraguay: Instituto de Investigaciones Biológicas del Paraguay (IIBP), Museo Nacional de Historia Natural de Paraguay (MNHNP). The database of Laboratorio de Genética Evolutiva (LGE, Instituto de Biología Subtropical, Misiones, Argentina) was also considered in the analysis. One hundred thirty five scientific articles were analyzed to extract presence points. When presence points were scarce, we completed species records with data obtained from Global Biodiversity Facility Information GBIF (accessed in March 2016). We cleaned this information to avoid mistakes from outdated taxonomic arrangements. All presence records are available in Supplemental Information 1.

We confirmed the species identifications when it was possible to avoid primary source of error in niche modeling. We identified the species recognizing external morphology features proposed in literature as descriptions of species, taxonomic and phylogenetic revisions (e.g. Heyer, 1970; Heyer, 1978; Heyer, 2014; de Sá et al., 2014). We used a binocular stereoscope when it was necessary.

We used the geographic coordinates of each specimen to georeference species occurrences, using Google Earth whenever these coordinates were missing. We estimated the uncertainties of these locations by considering the extent of the locality: each locality was defined as a circle, with a point marking the most likely position (the geographic center of the named place), and a radius representing the maximum distance from the point within which the locality is expected to occur (Wieczorek, Guo & Hijmans, 2004). The database was cleaned of records which were outside the known altitudinal range of each species and whose uncertainty was higher than 10,000 m. Sampling bias in geographic space (i.e. localities which are overrepresented in the collections) may lead to biases in the inferences about the environmental requirements of the species. To avoid this spatial bias and its consequent autocorrelation, we disaggregated occurrences of each species at 10,000 m with ‘ecospat’ package (Broenniman et al., 2014; Boria et al., 2014) implemented in R software. The conserved occurrence points were plotted in maps to check the consistency with the known ranges of the species. We analyzed a total of 2,306 localities, and the number of occurrences per species varied between 33 for L. laticeps and 467 for L. latrans (Table 1).

Ecological niche model and environmental characterization

To model the presence of a species, we used the theoretical Biotic-Abiotic-Mobility approach (BAM) of Soberon & Peterson (2005), which captures and links the geographic and environmental dimensions of species distributions (Peterson & Soberón, 2012). The BAM model is composed of three components: 1) “B” components are the biotic conditions, i.e. the appropriate suite of both present species (e.g. food) and certain absent species (e.g. strong competitors, diseases). In general, this component is not included in the modeling process due to the fact that it is very difficult to make accurate spatial quantification of this kind of data (Barve et al., 2011). According to this and considering that Leptodactylus species exhibit a broad food spectrum, for which considering the component B becomes much more difficult, we do not included it in the present work. 2) The “A” components represent the abiotic conditions, e.g. bioclimatic variables; and 3) The “M” component is the region of the world which has been accessible to the species via dispersal over relevant periods of time. This heuristic scheme assumes that stable populations of a species will be found only in the intersection of the B, A, and M components, (B∩A∩M) (Soberon & Peterson, 2005).

To delineate the “A” component, we selected environmental variables (i.e. climate and soil variables) matching the temporal and geographic resolution of the studied species occurrences. Climate variables refer to temperature and precipitation. Since anurans are ectotherms, they are more vulnerable to environmental changes than other tetrapods (Duellman & Trueb, 1994). Temperature is the most pervasive factor affecting the rate of biological reactions and physiological processes of amphibians (Rome, Don Stevens & John Alder, 1992); while precipitation is the most important extrinsic factor controlling the seasonality of reproduction in anurans (Duellman & Trueb, 1994). As a matter fact, in most of the species of temperate latitudes, reproductive activity depends on temperature and precipitation (Duellman & Trueb, 1994). Soil was hypothesized to be an important predictor of the distribution of arid-adapted anurans, since it can affect water loss rates during aestivation (Schalk, Montaña & Springer, 2015). Soil pH influences hatching success in Pseudophryne bibronii (Chambers, Wilson & Williamson, 2006), thus, soil characteristics could affect the success of the survival rate of terrestrial egg-laying species, such as most Leptodactylus species. Climate information was obtained from WorldClim database (Hijmans et al., 2005), which includes 19 variables of temperature and precipitation. Soil information was obtained from SoilGrids1km database (Hengl et al., 2014, through ISRIC–WDC Soils), which includes eight global soil variables that summarize different aspects of soil at six depths. Most of these variables are not independent; consequently, Pearson correlations were used to detect and exclude highly correlated variables (r ≥ 0.8). For soil variables, correlation analyses were made within each depth and between depths (2.5, 22.5 and 45 cm depth). When high correlations were detected, we kept the variables that most likely affect the ecology of anurans. Eight bioclimatic variables were retained: BIO1 = Annual Mean Temperature, BIO2 = Mean Diurnal Range (Mean of monthly (max temp − min temp)), BIO4 = Temperature Seasonality (standard deviation of mean month temperature × 100), BIO5 = Maximum Temperature of Warmest Month, BIO6 = Minimum Temperature of Coldest Month, BIO12 = Annual Precipitation, BIO13 = Precipitation of Wettest Month, BIO14 = Precipitation of Driest Month. Because the correlations of each soil variable between depths were high (r > 0.9), only the variables of the most superficial depth (2.5 cm) were kept. Six soils variables were retained: 1) Sand content, 2) Coarse Fragments, 3) Soil organic carbon, 4) pH index, 5) Bulk density, and 6) Cation-exchange capacity. The environmental databases used presented a spatial resolution of 2.5′ (approx. 5 × 5 km per pixel at the equator).

A critical step for calibration, validation and comparison of ENMs is the definition of the “M” component; i.e. the extent of the region where the model calibration will be performed. The “M” component must be customized for each species and must be related to an explicit a priori hypothesis (Barve et al., 2011). It is reasonable to take the set of biotic regions within which a species is known to occur as the hypothesis of the M region (Barve et al., 2011). Here, the M region was defined for each species by identifying the ecoregions (sensu Olson et al., 2001) with known presences and by buffering the resulting region to generate a soft edge. The buffer size was between 0.1° and 0.5°, depending on each species data. A buffer of 0.5° was applied when the points were located marginally in one ecoregion, and a smaller buffer (less than 0.5°), when the points were located more uniformly in the ecoregions. This prevents from including areas where the species had no access. We pruned the selected areas when they included some regions where species obviously did not occur. For example, since Leptodactylus species are mainly distributed in the lowlands, we excluded higher altitudes from the “M” regions according to the altitudinal range known for each species.

ENMs were calibrated using a maximum entropy method, MAXENT v3.3.3K (Phillips, Anderson & Schapire, 2006) in batch mode implemented in R-package “ENMGadgets” (Barve & Barve, 2016). Maxent fits models based on the probability distributions that show maximum entropy (i.e. closest to uniform), subject to the correlations of known presences of species with the environmental conditions across the study area (Phillips, Anderson & Schapire, 2006). Most default settings were kept but we defined five bootstrap replications and raw outputs (i.e. continuous probability values). We randomly partitioned occurrence data in two subsets: 50% of the occurrences were used to calibrate the model, and the remaining 50% were used for model evaluation. Given the observed limitations of the Receiver Operating Characteristics Curve (ROC) approach we used a partial ROC area under the curve (AUC) to evaluate the performance of the model (Lobo, Jiménez-Valverde & Real, 2008). Partial ROC allows for differential weighting of omission and commission errors and focuses on meaningful predictions for model evaluation (Peterson, Papeş & Soberón, 2008). We restricted ROC space to predictions corresponding to omission error ≤ 5%, and we then randomly sub-sampled 50% of the available evaluation data 1,000 times and estimated the ratios of simulated AUCs and null expectations of AUC (Peterson, Papeş & Soberón, 2008) using Partial-ROC software (Barve & Barve, 2016). Probabilities were determined by direct count of null replicate frequencies of AUC ratios falling below the observed value.

To identify the variables which contributed the most to the model, we used a jackknife test on the original environmental variables. This test removes one environmental variable from the full list at a time and recalculates the model to quantify the contribution of each variable to the overall model performance.

Spatial analyses

The resulting ENM can be projected in the geographic space, thereby depicting the potential distribution area of each species (Araújo & Peterson, 2012). We assumed that the species analyzed are in ecological equilibrium (i.e. their populations are not migrating) and that the potential distribution obtained corresponds to their historical distribution. The assumption of equilibrium of a population within its ecological niche (or single environment) implies ignoring mutation pressures from the environment and migration processes (Levene, 1953). State of stable equilibrium should not be viewed as a primary property of ecological systems, but it is a feature that can emerge from extrapolation to large spatial scales (DeAngelis & Waterhouse, 1987). The estimates of the geographic range extent obtained by ENM techniques have proved to be more representative than those obtained by traditional methods (e.g. minimum convex polygon (MCP)), particularly when data are scarce or when species are rare (Marcer et al., 2013; Pena et al., 2014; Syfert et al., 2014); and also avoid potential subjective bias of experts (Fourcade et al., 2013). We used a minimum threshold (Anderson, Lew & Peterson, 2003) to convert raw model outputs (i.e. continuous probability values) to distribution estimates (i.e. a binary representation). We modified this threshold–which, by default, includes 100% of the training data–to consider the potential error in the occurrence data, by choosing the highest threshold which included (100–E)% of the training data, assuming E = 5% (Table 2). E is an estimate of the proportion of the occurrence data that is likely to include the georeferenciation error (Peterson, Papeş & Soberón, 2008). This is a conservative method that minimizes the commission rate. This threshold defines that 95% of the observed presences are predicted as such, representing, in other words, the sensitivity of the model (Allouche, Tsoar & Kadmon, 2006). We then projected the distributional estimates of the ENM in space to obtain a quantitative potential geographic distribution. The spatial reference used for the analyses was Lambert Azimuthal Equal Area (LAEA) projection at a spatial resolution of approximately 5 × 5 km per pixel at the equator. The analyses and maps were developed with QGis software v. 2.16.2 (QGIS Development Team, 2016).

Ecological niche model and environmental characterization

The partial ROC tests indicated significant predictive ability of the models (p < 0.001). For all species, bioclimatic variables (temperature and precipitation) were the most important predictors. The temperature variables (temperature seasonality, the maximum temperature of warmest month and the minimum temperature of coldest month) were the most important variables for seven out of 10 species. Only for one species’ potential distribution (L. podicipinus) the precipitation (precipitation of wettest month) was important. For two species’ potential distribution a soil variable (PH index) was the most important (Table 2).

Spatial analysis

The potential distribution areas ranged from less than 40,000 km² for L. laticeps to 7.5 million km² for L. fuscus, with values between 1 and 3.5 million km² for most species (Fig. 1). Leptodactylus fuscus and L. latrans were the species with the largest potential distribution ranges (above 5 million km²) (Fig. 1).

Land cover and ecoregion characterization

The level of overlaps of the potential distribution areas of Leptodactylus species ranged from 30 to 72% with “forest habitats” and from 29 to 71% with “open habitats” (Table 3). Among forest habitats categories, the highest overlaps were found in Subtropical broadleaf deciduous forests, Tropical broadleaf evergreen forests and Subtropical broadleaf evergreen forests. Among open habitat categories, the highest overlap occurred with croplands.

Table 3:

Percentages of overlap of the distributions models of chacoan species of Leptodactylus with land cover categories from Blanco et al. (2013).

	Land cover	L. bufonius	L. chaquensis	L. elenae	L. fuscus	L. gracilis	L. laticeps	L. latinasus	L. latrans	L. mystacinus	L. podicipinus
Forest habitats	Tropical broadleaf evergreen forest	7.46	18.06*	28.21*	36.18*	0.36	1.93	0.30	21.12*	11.22*	34.70*
	Tropical broadleaf deciduous forest	0.55	0.27	0.47	0.15	0.00	0.06	0.00	0.08	0.09	0.29
	Sub-tropical broadleaf evergreen forest	5.87	8.56	9.79	6.37	23.38	0.68	6.65	9.34	11.28*	6.62
	Sub-tropical broadleaf deciduous forest	41.73*	14.95*	22.70*	5.03	5.00	69.12**	22.36*	6.73	10.99*	9.14
	Temperate broadleaf evergreen forest	0.00	0.00	0.00	0.01	0.03	0.00	0.01	0.03	0.02	0.00
	Temperate broadleaf deciduous forest	0.00	0.03	0.00	0.01	0.09	0.00	0.01	0.02	0.03	0.00
	Total	55.62**	41.87*	61.18**	47.74**	28.86**	71.80**	29.33**	37.32**	33.63**	50.75**
Open habitats	Tropical shrubland	0.09	6.70	6.86	15.09**	0.01	0.04	0.00	12.46**	7.70	9.55
	Tropical grassland	0.00	0.00	0.01	7.60	0.00	0.00	0.00	0.18	0.00	0.19
	Sub-tropical shrubland	17.32**	8.38	5.57	1.52	11.61**	10.07**	14.69**	4.69	7.28	4.61
	Sub-tropical grassland	0.30	1.80	0.01	0.36	10.29	0.00	10.61	2.60	3.80	0.36
	Temperate shrubland	0.00	0.11	0.00	0.00	0.83	0.00	0.03	1.90	1.17	0.00
	Temperate grassland	0.00	0.00	0.00	0.03	0.24	0.00	0.22	0.26	0.03	0.00
	Cropland	22.87**	36.25**	21.24**	24.99**	41.08**	15.17**	38.69**	36.27**	42.98**	29.24**
	Others	3.80	4.89	5.13	2.66	7.08	2.93	6.43	4.32	3.42	5.30
	Total	44.38**	58.13**	38.82**	52.26**	71.14**	28.20**	70.67**	62.68**	66.37**	49.25**

DOI: 10.7717/peerj.2605/table-3

Notes:

* Values show overlap percentages higher than 10%.

** Values show overlap percentages higher than 50%.

All the species considered in this study distributed beyond the Dry Chaco ecoregion (Fig. 2). Leptodactylus laticeps, L. bufonius, L. latinasus, L. elenae and L. chaquensis are mostly Dry chacoan species, with at least 25% of their geographic range occurring in the Dry Chaco. These species potentially inhabit less ecoregions than the other species (Supplemental Information 2). On the other hand, the remaining five species (L. latrans, L. fuscus, L. podicipinus, L. mystacinus, L. gracilis) potentially inhabit in a greater number of ecoregions (Supplemental Information 2).

Land transformation in Dry Chaco

By the end of 2013, 1.8 million km² of Dry Chaco (19% of the natural area in the entire ecoregion) were transformed to agriculture (Vallejos et al., 2014). The percentage of remaining natural areas inhabited by Leptodactylus chacoan species in the Dry Chaco decreased every year from 1976 to 2013 due to agricultural expansion, although the greatest transformation occurred after 1996 (Fig. 3). The largest percentage of potential chacoan area loss was 16% (L. mystacinus and L. podicipinus), while for seven species the loss of chacoan area were between 15% (L. bufonius, L. gracilis, L. latinasus and L. latrans) and 14% (L. chaquensis, L. fuscus and L. laticeps). The minimum loss of potential chacoan area was 13% (L. elenae) (Fig. 3).

Protected areas

The highest proportion of the potential distribution area in protected areas was found in Leptodactylus fuscus, L. latrans and L. mystacinus (between 6 and 11%), while the remaining species (L. bufonius, L. chaquensis, L. elenae, L. gracilis, L. latinasus, L. laticeps and L. podicipinus) exhibited between 3.5 and 5% of their potential distribution in protected areas (Fig. 4). The near threatened (Cortez et al., 2004) and vulnerable (Vaira et al., 2012) species Leptodactylus laticeps showed only 3.6% of its potential distribution in protected areas (Fig. 5).

Ecological niche model and environmental characterization

Our results show that the most important bioclimatic variable explaining species occurrence (in 70% of the studied species) was temperature. Other studies have also shown that climate (especially temperature and rainfall) is a determinant factor for amphibians (e.g. Parris, 2004; Urbina-Cardona, Olivares-Pérez & Reynoso, 2006). Since amphibians are ectothermic organisms, environmental temperature is probably a major constraint to their ecophysiological traits; for example, their mobility and their energy balance may be affected by temperature (Bennett, 1990). Physiological traits turn amphibians strongly sensitive to environmental conditions (Ferder & Burrgren, 1992), and consequently to climatic change (Pounds, Fogden & Campbell, 1999; Pounds et al., 2006). Projected climate changes over the 21st century include rises in surface temperature, and in the frequency and duration of heat waves and the occurrence of extreme precipitation events (IPCC, 2014). Also, seasonality will be affected by more frequent hot temperature extremes and fewer cold temperature extremes over most land areas (IPCC, 2014). Particularly, in south-eastern South America, which includes Dry Chaco ecoregion, changes in precipitation produce wetter conditions (Christensen et al., 2007). Under these scenarios, it is expected that populations change their geographic ranges to coldest areas, i.e. to higher altitudes or latitudes (Harsch et al., 2009; Thomas, 2010; Chen et al., 2011), although in the case of amphibians, this is unlikely due to their low dispersal ability. Measurements of the home range size of amphibians show a mean home range size of 40 m² (Wells, 2007). Moreover, the capacity for amphibian behavioral thermoregulation is limited because the cooling effect of evaporative water loss from skin counteracts heat gain by basking (Hutchinson & Dupe, 1992). Studies revealed that the effects of climate change over amphibians involve changes in the timing of breeding of some species (Blaustein et al., 2010; Carey & Alexander, 2003), instead of changes in their distributions (Corn, 2005). For example, no colonization of higher altitude areas associated with declines of Anaxyrus boreas has been found (Livo & Yeakley, 1997). In subtropical areas, the period of activity of amphibians and their population dynamics are strongly influenced by climate seasonality, e.g. reproduction behavior of most species occurs in the warm season of the year (Conte & Machado, 2005). For many species of Leptodactylus genus, the reproductive season matches with wet and/or warm season, e.g. L. latinasus, (Gallardo, 1964; Canavero et al., 2008), L. bufonius (Cei, 1980), L. gracilis, (Canavero et al., 2008; Ximenez & Tozetti, 2015), L. mystacinus, (Gallardo, 1964), L. fuscus, (Martins, 1988; Prado, Uetanabaro & Haddad, 2005; Lucas et al., 2008), L. elenae (Prado, Uetanabaro & Haddad, 2005), L. chaquensis (Gallardo, 1987; Prado, Uetanabaro & Haddad, 2005; Canavero et al., 2008) and L. podicipinus (Cei, 1980; Vizotto, 1967; Rossa-Feres & Jim, 1994); although it was reported that populations of L. podicipinus from the Brazilian pantanal exhibit a continuous reproductive cycle throughout the year (de Almeida Prado, Uetanabaro & Lopes, 2000; Prado, Uetanabaro & Haddad, 2005). The Dry Chaco presents a monsoonal climate, in which rains are concentrated in the warm season. Since the reproduction of most amphibians is associated to water, modifications in water balance in the area could have significant effects on their breeding timing and behavior. Thus, our results suggest that climate change could affect Leptodactylus chacoan species in different and unpredictable ways.

Spatial analyses