Research Paper
Research Paper
Modeling the potential distribution of floristic assemblages of high Andean wetlands dominated by Juncaceae and Cyperaceae in the Argentine Puna
expand article infoElvira Casagranda, Andrea E. Izquierdo§
‡ Universidad Nacional de Tucumán, Tucumán, Argentina
§ Universidad Nacional de Córdoba, Córdoba, Argentina
Open Access


Aims: The aim of this work was to model the distribution of suitable environmental conditions of vegas with specific floristic characteristics. Vegas are high Andean wetlands that represent the main sequestered carbon stocks, biodiversity hotspots, and water regulating systems in the region. In these wetlands, plant communities are the main biological factor that determines functional processes, and plant species assemblages are associated with different ecogeographic features. Study area: Argentine Central Andean Puna ecoregion. Methods: For two different floristic assemblages of vegas, we develop ecological niche models of n-dimensional minimum volume ellipsoids through NicheToolBox, then obtain potential distribution maps. One floristic assemblage was dominated by the cushion-structured plant Oxychloe andina (Juncaceae) and the other by plants of the Cyperaceae family. Results: Elevation and precipitation were the main environmental factors determining the distribution of the two floristic assemblages. Juncaceae dominated vegas tend to be located in high, humid, and cold places, while Cyperaceae vegas are found at a lower elevation, with less humidity, and higher temperatures. According to the dominant climatic gradient in the region, potential distribution maps show that vegas of Juncaceae are commonly found towards the Northeast of the Puna while Cyperaceae vegas are more frequent at lower elevations to the South of the region. Conclusions: This study represents the first approach to niche modeling based on plant communities in vegas of the Argentine Puna, providing knowledge on the environmental factors that limit their distribution. This information could serve as a planning tool in a region exposed to growing perturbations such as mining and climate change.

Taxonomic reference: Zuloaga et al. (2019).

Abbreviations: AUC = Area Under the ROC Curve; NDVI = Normalized Difference Vegetation Index; ROC = Receiver Operating Characteristic.


Central Puna ecoregion, Cyperaceae, floristic assemblages, Juncaceae, niche models, plant communities, vegas, vegetation distribution


Vegas (also known as bofedales; see White-Nockleby et al. 2021) are characteristic high-altitude wetlands along the Andes that develop under extreme environmental conditions such as high altitudes, daily frosts, and hypoxic soils (Squeo et al. 2006; Otto et al. 2011). These wetlands are one of the most productive ecosystems above 3,000 m a.s.l. (Squeo et al. 2006; Domic et al. 2021), with multiple and important ecosystem functions such as carbon sequestration, habitat provision for a large number of species, forage for native and domestic herbivores, as well as fresh water for human populations (Squeo et al. 2006; Maldonado Fonkén 2014; Izquierdo et al. 2018a). Vegas originate from groundwater aquifer outcrops in depressed areas and valley bottoms or slopes, constituting permanently irrigated plant communities, or at least for a large part of the year (Izquierdo et al. 2018a). In vegas, the vegetation is capable of modifying hydrology at the local level by retaining and regulating water flows, which combined with the lack of oxygen promotes the accumulation of organic matter below the surface (Ruthsatz 2012; Carilla et al. 2018; Izquierdo et al. 2018a). Various disturbances can alter the stability of vegas, including overgrazing (Prieto et al. 2003; Domic et al. 2018), mining (Izquierdo et al. 2015a, 2018b), tourism (Izquierdo et al. 2018a; Troncoso, 2018), and climate change (Carilla et al. 2013; Morales et al. 2015, 2018); factors that present complex relationships and different spatial-temporal trends (Izquierdo et al. 2018b; Navarro et al. 2020).

In these wetlands, plant species dominance gradually changes north-south throughout their distribution, and is determined by ecological factors and the different histories of colonization of the species (Ruthsatz 2012). Some cushion and cespitose species such as Distichia muscoides (Juncaceae) are widely distributed in the tropical and subtropical Andes from Colombia to northern Chile and northwestern Argentina (Carilla et al. 2018; Ruthsatz et al. 2020; Benfield et al. 2021). While other species such as Oxychloe andina (Juncaceae), Zameioscirpus and Carex (Cyperaceae), are more restricted (POWO 2022). Cushion plants are considered to be vega forming species (i.e. ecosystem engineers) since they are fundamental in determining soil properties, and thus their diversity and composition (Badano and Cavieres 2006). In addition to the dominant ones, there are also species of accompanying flora, which colonize typical micro-environments such as ponds and watercourses, sand and gravel sediments (Navarro et al. 2011; Ruthsatz 2012; Ruthsatz et al. 2020). The accompanying species are able to settle in these extreme environments because the cushion plants provide them water and protection against frost, herbivory, and erosion, among others (Ruthsatz 2012). Therefore, vegas represent ecosystems of great taxonomic plant richness, including endemic species and habitat specialists (Ruthsatz 2012; Polk et al. 2019; Izquierdo et al. 2020; Ruthsatz et al. 2020) which characterize these environments (Ruthsatz 2012; Maldonado Fonkén 2014; Izquierdo et al. 2020).

Several studies have shown that the distribution of the vega forming species seems to respond to changes in environmental conditions, indicating that they present a high environmental specificity, as is the case of the Andean Oxychloe andina cushion plant, which is commonly found in high, wet and cold vegas (Ruthsatz 2012; Izquierdo et al. 2020; Ruthsatz et al. 2020; Domic et al. 2021). In contrast, species of the accompanying flora present greater ecological amplitude and therefore a wider distribution, such as species of the genus Eleocharis (Izquierdo et al. 2020; Ruthsatz et al. 2020). Other plant species, such as halophytes, predominate in vegas with soils with higher salinity, adapted to lower altitude, drier and warmer conditions (Montesinos 2012; Ruthsatz 2012; Izquierdo et al. 2020). Soil moisture and organic matter content have been found to play an important role in determining the plant community formation in high Andean wetlands, showing a strong differentiation in the species composition between flooded and dry wetland communities (Domic et al. 2021). Another important variable that should be considered is pH, which shows a latitudinal variation with more alkaline wetlands found in the tropical regions of the Central Andes and of lower pH towards the southern end of its distribution (Ruthsatz 1993; Ruthsatz et al. 2020). Also, the depth and stability of the groundwater, which in turn determines the accumulation of salts, are important in shaping the composition and richness of these communities (Navarro 2020).

Ecological niche modeling is useful in estimating the environmental requirements of the species, which can be projected geographically to help identify areas suitable for their survival (Barve et al. 2011). Niche models link the occurrence of species in known locations with environmental and spatial information, making it possible to generate spatial models of habitat suitability (Hirzel and Le Lay 2008), thus providing information on the potential distribution of species (Lobo et al. 2015). Although most of the studies have traditionally focused on modeling the individual species niche, biodiversity loss and global climate change have highlighted the relevance of expanding the niche concept to groups of functionally closely associated species (Hirzel and Le Lay 2008). The modeling of individual species takes into account their particular needs and acquires special importance when applied to species of conservation interest (i.e. threatened or emblematic species), but it does not directly address the pattern of biological diversity as a whole, particularly when dealing with diverse taxa in sparsely studied regions: there will be many species with few records that will not allow them to be effectively modeled individually (Ferrier et al. 2002). One way to solve this is by integrating spatial models with numerical classification techniques that analyze patterns in data sets, where the data matrix (presences, absences, or relative abundances) is classified into groups of species found at similar sites (or groups of sites with similar species), which can then be modeled and extrapolated across the region of interest (McKenzie et al. 1989; Ferrier et al. 2002). By relating the distribution of multiple species to sets of environmental variables, these models analyze the collective properties of species distribution (Ferrier and Guisan 2006).

The present contribution provides an initial model that predicts the environmental variables that determine the ecological niche of different floristic assemblages of vegas at regional scale. For this, we: a) model the potential niches of two floristic assemblages of vegas, and b) analyze their potential spatial and ecological distribution into the study area. This approach provides information about environmental requirements of ecological engineer species and allows us to identify potential sites throughout the region with the environmental conditions that these vegas need to exist.

Study area

The study area includes the Argentine sector from the Central Andean Dry Puna, Central Andean Puna and part of Southern Andean Steppe of the terrestrial ecoregions of the world (Olson et al. 2001) covering an area of approximately 12,000,000 ha (Fig. 1). It is delimited to the north by part of Jujuy province, to the west by the international borders between Argentina and Bolivia, and Argentina and Chile, the southern limit is defined by the San Guillermo Biosphere Reserve in the province of San Juan; while the eastern limit is given by the altitudinal level of 3,200 m a.s.l. (Izquierdo et al. 2015b). The study area contains 8,593 vegas, which occupy an area of 54,880 ha, representing 0.46% of the total study area. Vegas are located in an altitudinal average range between 2,500 and 5,389 m a.s.l. Almost all vegas (99.6%) are small, with sizes between 1 and 100 ha, 0.39% have between 100 and 1,000 ha, and only one exceeds 2,000 ha. The region’s water resources are conditioned by an orographic control of rainfall, as the north-south mountain ranges act as a barrier to the humid Atlantic winds. This results in a climate characterized by aridity, with an average annual rainfall of less than 400 mm that occurs mainly during the summer (Cabrera 1976; Morales et al. 2018), and decreasing in a Northeast-Southwest direction, with sectors with less than 100 mm per year (Paoli et al. 2002; Reboratti 2005). In addition to scarce rainfall, the combination of extreme climatic variables such as evapotranspiration associated with strong winds and high solar radiation results in a negative water balance throughout the year (Aceituno 1993; Izquierdo et al. 2018a).

Figure 1. 

Study area modified from Izquierdo et al. (2015b). Ecoregions are based on Olson et al. (2001).


Niche model data

For this study we used the classification of vegas that was obtained by Izquierdo et al. (2022) which was based on the vegetation sampling of 50 vegas distributed throughout the study area, located between 3,277 and 4,827 m a.s.l. In the sampling, the composition of vascular plant species and their percentage cover were recorded using a quadrat of 1m2 divided into 100 sub-quadrats of 10 cm2, equivalent to 1% of the total. The quadrats were distributed to capture the greatest diversity of vegetation, ranging from six to 41 quadrats in each sampled vega, depending on their size and environmental heterogeneity. The identification of the plants at the species level was verified by specialists from the National University of Tucumán based on Flora Argentina (Zuloaga et al. 2019). A total of 111 plants were identified at the species level (only one at the genus level) belonging to 67 genera and 28 families. The highest diversity was found in a vega with 35 species, and the lowest in a vega with seven species, being 17 the average value of species per vega (Izquierdo et al. 2022). In that study, through a Correspondence Analysis (Benzecri 1992) and k-means analysis (Likas et al. 2003) five floristic assemblages of vegas were grouped. These floristic assemblages were characterized by spatial and spectral variables that represent the ecological and geographical context in which these communities exist. The grouping was explained in two dimensions, one driven by altitude and stability of vegetation productivity, represented through the amplitude and maximum value of the NDVI, and the other by longitude and geographical latitude, salinity and humidity.

For niche modeling we built two groups using data belonging to three of these five floristic assemblages, covering a total of 41 vegas (Fig. 2). The two remaining floristic assemblages are represented by few vegas (seven and two vegas respectively), resulting in too few occurrence points to build niche models, and were therefore discarded in this analysis. One of the groups we used for modeling (formed by floristic assemblages 1 and 2 of Izquierdo et al. 2022) is dominated by cushion plants of the Juncaceae family (Oxychloe andina) (Table 1) associated mainly with Poaceae, totaling 17 vegas (hereafter “Juncaceae vegas”, Fig. 2). Juncaceae vegas are typically found at stream headwaters, forming dense and continuous vegetation tapestries (Carilla et al. 2018). The second group modeled (floristic assemblage 3 of Izquierdo et al. 2022) contains 24 vegas dominated by species of the families Cyperaceae and Campanulaceae, among others (Table 1), with Eleocharis pseudoalbibracteata and Zameioscirpus atacamensis (both Cyperaceae) being the dominant species of the group (hereafter “Cyperaceae vegas”, Fig. 2). Cyperaceae vegas are generally found in lower altitude sectors, where there is increased salinity or are associated with salt flats (Carilla et al. 2018).

Figure 2. 

On the left side the map of the studied area is shown with the location of the vegas classified by Izquierdo et al. (2022), grouped into Juncaceae and Cyperaceae. The right panel shows panoramic views of the vegas a. Tocomar (VG001); b. Chorrillos (VG002) and c. Colorada (VG026), with Tocomar and Chorrillos representing Juncaceae vegas and Colorada Cyperaceae vegas. Some typical species of these floristic assemblages are shown in d. Oxychloe andina (Juncaceae vegas), e. Lobelia oligophylla (Cyperaceae vegas), f. Zameioscirpus atacamensis (Cyperaceae vegas) and g. Triglochin concinna (Juncaceae and Cyperaceae vegas). Photos: Julieta Carilla.

Table 1.

Vega groups classified by Izquierdo et al. (2022), here grouped into Juncaceae vegas (floristic assemblages 1 and 2), and Cyperaceae vegas (floristic assemblage 3). The ten species with the highest mean percent cover within each group are shown. The family to which each species belongs is indicated in parentheses.

Juncaceae vegas Cyperaceae vegas
Group 1 species Mean cover (%) Group 2 species Mean cover (%) Group 3 species Mean cover (%)
Oxychloe andina (Juncaceae) 27.52 Oxychloe andina (Juncaceae) 34.48 Eleocharis pseudoalbibracteata (Cyperaceae) 14.19
Festuca nardifolia (Poaceae) 9.30 Zameioscirpus atacamensis (Cyperaceae) 8.16 Zameioscirpus atacamensis (Cyperaceae) 9.62
Deyeuxia hackelli (Poaceae) 6.58 Triglochin concinna (Juncaginaceae) 3.99 Juncus balticus (Juncaceae) 7.43
Distichia muscoides (Juncaceae) 5.21 Deyeuxia eminens (Poaceae) 3.73 Eleocharis atacamensis (Cyperaceae) 4.90
Trichophorum rigidum (Cyperaceae) 5.07 Eleocharis pseudoalbibracteata (Cyperaceae) 3.40 Triglochin concinna (Juncaginaceae) 4.38
Zameioscirpus muticus (Cyperaceae) 3.44 Eleocharis atacamensis (Cyperaceae) 3.10 Oxychloe andina (Juncaceae) 4.30
Deyeuxia curvula (Poaceae) 2.67 Deyeuxia curvula (Poaceae) 2.13 Festuca argentinensis (Poaceae) 4.17
Phylloscirpus deserticola (Cyperaceae) 2.36 Festuca argentinensis (Poaceae) 1.83 Lobelia oligophylla (Campanulaceae) 3.96
Zameioscirpus atacamensis (Cyperaceae) 1.79 Juncus stipulatus (Juncaceae) 1.62 Phylloscirpus acaulis (Cyperaceae) 2.36
Rockausenia pygmaea (Asteraceae) 1.51 Calandrinia acaulis (Montiaceae) 1.39 Carex macrorrhiza (Cyperaceae) 2.16

Ecological niche modeling and environmental characterization

We modeled the ecological niche of the floristic assemblages using the biotic-abiotic mobility (BAM) theoretical approach (Soberón and Peterson 2005), in which environmental and geographic dimensions of species distribution are linked (Peterson and Soberón 2012). There are three components in the BAM model. Component B refers to the biotic conditions necessary for the maintenance of populations (i.e. food availability, presence and influence of competitors and predators). Component A represents the abiotic conditions necessary for species survival and growth, and is considered independent of species abundance or presence (i.e. bioclimatic variables). Finally, the M component refers to the geographic region that has been accessible to the species for relevant periods of time.

Despite the importance of biotic interactions in determining the distribution of species and species assemblages at regional, continental and global scales (Wisz et al. 2013), these are generally not included in niche modeling of individual species due to their complex nature, which makes it difficult to spatially quantify this type of data (Barve et al. 2011). However, when modeling groups of species it is possible to consider that biotic interactions between them are implicitly considered (Baselga and Araújo 2009). To determine component A, we used 19 bioclimatic variables provided by WorldClim 2.0, corresponding to the average of the period 1970–2000 (Fick and Hijmans 2017), with spatial resolution of 30” (~1 km2) (Table 2), and the variable altitude with the same spatial resolution. Climatic variables, especially over large areas, are among the most important factors modeling species distribution (Grinnell 1917; Guisan and Theurillat 2000), as they directly influence the physiology of organisms. In the case of plants, these variables are of particular importance since plants cannot evade adverse climatic conditions by moving or sheltering (Hirzel and Le Lay 2008). Particularly, WorldClim bioclimatic variables are widely used in ecological modeling as they offer a better fit than monthly or annual averages (Hirzel and Le Lay 2008). On the other hand, topography affects species indirectly through its correlation with temperature and precipitation, thus topographic variables are crucial for plants (Guisan et al. 1998). The model calibration area or “M” was the study area previously described, taking into account that the species that compose the floristic groups have known distribution within the ecoregions of the study area and within their altitudinal limits.

Table 2.

Variables used for niche modeling. WorldClim 2.0 bioclimatic variables are derived from temperature and precipitation data for the period 1970–2000.

BIO1: mean annual temperature BIO11: average temperature of the coldest quarter
BIO2: daytime temperature range BIO12: annual precipitation
BIO3: isothermality (BIO 2/ BIO 7)*100 BIO13: precipitation of the rainiest month
BIO4: seasonality of temperature (σ*100) BIO14: precipitation of the driest month
BIO5: maximum temperature of the warmest month BIO15: seasonality of precipitation
BIO6: minimum temperature of the coldest month BIO16: precipitation of the rainiest quarter
BIO7: annual temperature range (BIO 5 - BIO 6) BIO17: precipitation of the driest quarter
BIO8: average temperature of the rainiest quarter BIO18: precipitation of the warmest quarter
BIO9: average temperature of the driest quarter BIO19: precipitation of the coldest quarter
BIO10: average temperature of the warmest quarter Altitude

To perform the niche models we used the NicheToolBox package (Osorio-Olvera et al. 2020a) implemented through scripts in the R program (R Core Team 2020). NicheToolBox allows the estimation of ecological niches using different algorithms such as BIOCLIM, Maxent and Minimum Volume Ellipsoids. In this study we generated the niches using n-dimensional minimum volume ellipsoids (MVEs) (Van Aelst and Rousseeuw 2009), a function that uses Mahalanobis distances to the centroid of the ellipsoid with the idea that the maximum environmental suitability occurs at this centroid (Osorio-Olvera et al. 2019). The structure of fitness in niche space has been hypothesized to be approximately ellipsoidal, with some empirical support for this idea (Maguire 1973; Osorio-Olvera et al. 2020b). Ellipsoids are simple models that require three parameters to be defined: (a) a niche-centroid, which is the point in ecological space where fitness has maximum value (Martínez-Meyer et al. 2013); (b) a shape matrix, which measures how dependent two niche axes are and how they change together; and (c) the proportion of observations to be included in the ellipsoid (Van Aelst and Rousseeuw 2009). NicheToolBox allows users to calibrate models and perform model selection based on statistical significance (partial ROC and AUC) and model predictive performance (omission rates). The partial ROC (Receiver Operating Characteristic) test has been modified to improve on the classical ROC test (Peterson et al. 2008), being a statistical significance test more appropriate for modeling algorithms using only presence data, as in this study, and giving more weight to omission errors than to commission errors (Peterson et al. 2008). The AUC (Area Under the ROC Curve) measures the area under the ROC curve and its values vary from 0 to 1, indicating that the model is statistically better than chance the closer its value is to 1.

The presence points we used correspond to the coordinates of the center of each vega in each of the groups, obtaining two sets of occurrence data, one for the Juncaceae (17 points) and another for the Cyperaceae group (24 points). Each occurrence dataset was randomly divided into training and test data in a 70:30 ratio, and environmental information was extracted for both training and test points. For each set the correlations between the environmental variables were estimated in order to identify the variables that were least correlated to each other to avoid redundant information. With these variables the models were fitted, specifying the number of them to be used (i.e. the number of dimensions in which the ellipsoid models will be constructed). Other important parameters that were defined before running the models are the proportion of training points that will be used to fit the model (Van Aelst and Rousseeuw 2009), in our study defined as 0.99; the collection of information from the environmental layers to calculate the statistical significance of the models (AUC and partial ROC), for which we used 50,000 points. Finally, we specify the omission rate below which the best models will be selected; in this case we defined an omission rate of 6% (accepted below 10%). Models were built with all the parameters previously described, and then the best model was selected according to the omission rate for training and test data below 6%, partial ROC test with p-value ≤ 0.5, and maximum AUC value. A scheme of the workflow is presented in Fig. 3.

Figure 3. 

Steps performed for ecological niche modeling and obtaining potential distribution maps using NicheToolBox.

Once the best model was chosen for each vegas group, minimum volume ellipsoid and the map of potential habitat suitability in the geographic space (or potential distribution map) were obtained. Habitat suitability varies between 0 and 1, with 1 indicating that the cell contains environmental conditions more similar to the sites where the vegas are present. The suitability map was clipped with the vega layer of the studied area, and then a cut-off threshold was applied to conserve all those vegas that presented suitability ≥ 0.1, considering that beyond this threshold, there are already suitable conditions for a vega to be present. From the maps obtained after applying the threshold, for each group, the number of vegas with suitability ≥ 0.1, surface area, average, minimum and maximum altitude were obtained.

Data availability

Data on the species that compose the plant communities of the studied vegas at this study are openly available in GBIF API at, reference number


Niche modeling and vegas potential distribution

The environmental variables used to fit the models (i.e. the least correlated with each other) in Juncaceae vegas were: altitude, mean annual temperature (BIO1), diurnal temperature range (BIO2), isothermality (BIO3), minimum temperature of the coldest month (BIO6), annual precipitation (BIO12), precipitation of the driest month (BIO14) and seasonality of precipitation (BIO15). In Cyperaceae vegas, the variables used were: altitude, mean annual temperature (BIO1), isothermality (BIO3), seasonality of temperature (BIO4), annual temperature range (BIO7), annual precipitation (BIO12), precipitation of the rainiest month (BIO13) and seasonality of precipitation (BIO15). For both groups, ellipsoid models were constructed in three, five and seven dimensions using all possible combinations of the mentioned least correlated variables, generating 120 models in both cases. The best niche model for Juncaceae vegas obtained was constructed with three variables: mean annual temperature (BIO1), diurnal temperature range (BIO2) and annual precipitation (BIO12), whose values at the centroid were 5.3°C (BIO1), 18.5°C (BIO2) and 108 mm (BIO12). For Cyperaceae vegas, the best model was built with five variables: altitude, mean annual temperature (BIO1), seasonality of temperature (BIO4), annual temperature range (BIO7) and seasonality of precipitation (BIO15); the centroid of this model was located at 3,762 m a.s.l. (altitude), 7.2°C (BIO1), 3.2°C (BIO4/100), 23.4°C (BIO7) and 97.4 mm (BIO15). Of the 120 models obtained in each group, two of the models generated for Juncaceae vegas exceeded the omission rate criterion of less than 6% for the training data, while none of the models exceeded this criterion for the test data, so the selection criterion was to choose the one with the highest AUC among the two models that exceeded the omission rate for training data. In the case of Cyperaceae vegas, 51 of the 120 models passed the omission criterion for both training and test data. The parameters of the best selected models are shown in Table 3.

Table 3.

Best minimum volume ellipsoid models for Juncaceae and Cyperaceae groups using the model selection and calibration protocol in NicheToolBox. For each group, the model that had omission rates ≤ 0.06, significant partial ROC value (p < 0.01) and the highest AUC value are shown.

Vegas group Best model variables (BIO) Omission rate (training) Omission rate (test) Partial ROC p-value AUC Total number of calibrated models
Juncaceae 1, 2, 12 0.00 0.28* 0.00 0.71 120
Cyperaceae Altitude, 1, 4, 7, 15 0.05 0.00 0.00 0.85 120

Potential distribution maps for both vegas groups are shown in Fig. 4. The best model for Juncaceae has a total of 1,439 vegas with suitability values ≥ 0.1 (Fig. 4A), which occupy an area of 15,418 ha, representing 28.0% of the total vega area in the study area (54,880 ha). These vegas are located at an average mean, maximum and minimum altitude of 4,423, 5,160 and 3,873 m a.s.l respectively. Of the 1,439 vegas, 1,051 have suitability values between 0.1–0.4, while 388 between 0.41–1.0. For Cyperaceae vegas, the best model encompasses 1,397 vegas with suitability values ≥ 0.1 (Fig. 4B), occupying 19,016 ha, representing 34.6% of the vegas surface of the study area. The average mean altitude of these vegas is 3,922 m a.s.l; while maximum and minimum is 4,389 and 3,252 m a.s.l respectively. Regarding their distribution in suitability ranges, 1,243 present values between 0.1–0.4, and 154 between 0.41–1.0. In addition, we found that between the models of both groups there is an overlap of 170 vegas, i.e. the same vegas have suitability values ≥ 0.1 in both models. Of these vegas, 134 have higher suitability values in the Juncaceae model than in the Cyperaceae one, while in the remaining 36 vegas the values are higher for the Cyperaceae group.

Fig. 5 shows the distribution in relation to mean altitude and mean annual precipitation values of all vegas in both modeled groups (Juncaceae and Cyperaceae) with suitability values ≥ 0.1.

Figure 4. 

Potential distribution maps for Juncaceae (A) and Cyperaceae (B) vegas according to habitat suitability threshold ≥ or ≤ 0.1.

Figure 5. 

Distribution of vegas predicted by the models performed for both groups (habitat suitability ≥ 0.1), according to their average altitude and annual precipitation values.


In this study we modeled the ecological niche and potential distribution of two floristic assemblages of the Argentine Central Andean Puna vegas according to the environmental variables that determine their distribution ranges. Niche modeling of functionally close species groups allows us to study the niches and potential distribution of these plant communities even with relatively few data (Ferrier et al. 2002), which is especially important in isolated or difficult to access study areas such as high mountain deserts, and with much lower density of data on the location of species compared to other ecoregions. Models obtained by using this methodological approach concur with previous studies (Ruthsatz 2012; Carilla et al. 2018; Izquierdo et al. 2020) and confirm the pattern shown in Izquierdo et al. (2022), which was the grouping base of the floristic assemblages used for this analysis. In short, sites with suitable habitat conditions for Juncaceae vegas are found at sites of higher elevation and generally higher rainfall than for Cyperaceae, according to different eco-geographic conditions related to altitude and humidity which are determinant factors in their spatial distribution (Fig. 5).

On the other hand, NicheToolBox has performed well in predicting sites with ideal environmental conditions for the existence of plant communities in both studied floristic assemblages. The reliability of the models is indicated by the relatively high AUC values reported, 0.71 and 0.85 for Juncaceae and Cyperaceae vegas respectively, with which it is possible to consider that they identify the sites where these vegas have been reported quite efficiently (Phillips et al. 2006; Ortíz-Yusty et al. 2014). In addition, and despite the fact that the Juncaceae model did not pass the selected criterion for omission rate in the test data, the partial ROC test showed a statistically significant value (p < 0.5) (Table 3). Although the AUC values are not particularly high, it is important to consider that this statistic is known to be method dependent (Peterson et al. 2008). It is important to consider not only the number of presence points used for modeling, but also that they correspond to the centroid of each vega, so that in vegas with wide altitudinal range there may be an important variation in environmental conditions that might not be contemplated in a single central point, even when the vegetation sampling was designed considering the heterogeneity of the vegetation to obtain a good representation of the floristic composition. In addition, with respect to elevation as an important factor in determining the distribution of vegas, it is possible that, in the case of small vegas, the spatial resolution of the data (1 km2) does not always take into account the range of altitudinal variation of a single pixel. Regarding other possible limitations for niche modeling, the use of WorldClim bioclimatic variables, although they are widely used in this field of ecology due to their biological importance (Hirzel and Le Lay 2008), in these high mountain regions climate models and their derivatives (such as bioclimatic variables) are poor due to the lack of instrumental climate records (Carilla et al. 2013; Casagranda et al. 2019). However, despite the limitations of the climatic data, it is interesting to highlight how these systems that are considered azonal, where vegetation depends mainly on local edaphic conditions (Ahumada and Faúndez 2009), are shown to be related to climatic conditions that are not strictly local.

The environmental suitability map for Juncaceae vegas places them in sites mainly located in the NE region of the study area, while that of Cyperaceae shows a potential distribution that covers areas of lower altitudes and reaches sites further south in the region (Fig. 4). These potential distribution maps for both groups coincide with those reported in Izquierdo et al. (2022), where the different groups show an association with spatial and spectral variables in an ecological gradient from wetter regions in the NE to more arid zones in the SW, and also related to a topographic gradient of elevation and latitude/longitude. That study found that the Juncaceae group is characterized by being in high zones at lower latitudes and longitudes, while the Cyperaceae group is found in less elevated zones and with higher minimum NDVI values.

In general, organisms respond to complex interactions between environmental variables (Rydgren et al. 2003), and for plants in particular, niches are affected to a greater extent by interactions between climatic variables (Huntley et al. 1989; Prentice et al. 1991). As proposed by Maguire (1973), the internal structure of ecological niches has a centroid where the optimal environmental conditions for species survival are found. The niche models obtained here showed that the centroid was determined by climatic variables in both groups of vegas, and in Cyperaceae also by altitude, a topographic variable with a strong influence on precipitation and temperature (Guisan et al. 1998) and an important factor in determining the plant composition of these wetlands (Ruthsatz 2012). In the case of the selected model for Juncaceae vegas, niche centroid is determined by mean annual temperature (5.3°C), diurnal temperature range (18.5°C) and annual precipitation (108 mm), indicating that Juncaceae-dominated communities establish well in sites with wet, cool conditions and considerable diurnal thermal amplitude. This is in agreement with work that has reported that Oxychloe andina (the dominant cushion plant of this group) is highly likely to be found in vegas located on higher, wetter and cooler sites compared to other vegas (Ruthsatz 2012; Izquierdo et al. 2020; Ruthsatz et al. 2020), and is also an obligate species of flooded areas (Domic et al. 2021). It has also been observed that O. andina grows in places where the water table is shallow (20–25 cm) and with great seasonal variability (> 40 cm) (Navarro 2020), being a species relatively resistant to salinity (Ruthsatz 2012). It would be expected to find suitable conditions for its establishment and that of its associated flora, in sites with environmental conditions that may differ from those considered by the model in the determination of the ecological niche. For the Cyperaceae vegas, on the other hand, the centroid of the best niche model was formed by five variables: altitude (3,762 m.a.s.l), mean annual temperature (7.2°C), seasonality of temperature (3.2°C), annual temperature range (23.4°C) and seasonality of precipitation (97.4 mm); indicating that the ideal habitat conditions are less cold and of lower altitude than for the Juncaceae vegas, with little difference in temperature but great difference in precipitation between seasons. These conditions reflect in part what has been previously reported for the two species with the highest percent cover in the Cyperaceae group: Eleocharis pseudoalbibracteata and Zameioscirpus atacamensis. E. pseudoalbibracteata is a species that, although it has a wide distribution, is generally found in vegas located at lower altitudes and in drier and warmer conditions (Izquierdo et al. 2020); while Z. atacamensis is a cushion plant that predominates in sites with a high saline content and lower altitudes (almost always below 4,200 m.a.s.l., Ruthsatz 2012) than O. andina; however O. andina can also be found in this group of Cyperaceae, although less frequently (Ruthsatz 2012; Ruthsatz et al. 2020). The Cyperaceae group was the largest, most diverse and species-rich of the groups classified by Izquierdo et al. (2022), so it is expected that the environmental conditions that determine the niche centroid are also more variable and diverse.

The overlap of the niche models in predicting vegas with suitability values greater than 0.1 for both groups could be due to the fact that these groups share both dominant cushion plant species of the communities (e.g. Oxychloe andina and Zameioscirpus atacamensis) and accompanying flora (e.g. Eleocharis pseudoalbibracteata, Triglochin concinna). In addition, these vegas classified to both floristic assemblages contain at least 20 vegas with an altitudinal variation of more than 100 to 200 m, which could be partly the reason why both floristic assemblages can be found in these vegas. This overlap can also be understood taking into account the gradual replacement of dominant species (and therefore associated flora) along the distribution gradient of these wetlands in the tropical and subtropical Andes (Ruthsatz 2012), and which occurs along the previously mentioned gradient of environmental conditions of decreasing humidity in a NE-SW direction as described by other authors (Izquierdo et al. 2015a; Casagranda et al. 2019; Ruthsatz et al. 2020; Izquierdo et al 2022). Species replacement between vegas groups is high, from vegas dominated by O. andina at the wetter end of the gradient, with a high percentage cover of species growing in extreme conditions, towards floristic assemblages of vegas that were not modeled in this study but have low productivity and include species adapted to saline conditions such as the cushion plant Amphiscirpus nevadensis (Izquierdo et al. 2022), towards the more arid end. In this gradient, the Cyperaceae vegas, where the dominant cushion plant is Z. atacamensis, are located in the middle of the two extremes, with the largest number of species found only in this group (21 sp.) and also sharing species with the groups of vegas located at the extremes of the environmental gradient (e.g. O. andina in wetter and colder conditions and T. concinna, a species more associated with halophyte communities in the arid extreme) (Izquierdo et al. 2022).

Despite the importance of climatic variables in determining the distribution ranges of the different types of plant communities in vegas, in the future it would be interesting to have other variables that contribute to the understanding of their habitat requirements. It has been shown, for example, that in similar wetlands in the northern hemisphere, the distribution and ecological niche of Cyperaceae species are not only determined by the climatic gradient and the chemical components of the water (especially salinity levels), but also by shade and the level of the water table (Gignac et al. 2004). Regarding the biotic interactions that occur within communities and play an important role in shaping species distributions, it is possible that interactions such as herbivory and competition (Ruthsatz 2012) play a fundamental role in the spatial configuration of communities, and although it has been postulated that multi-species distribution models can capture these types of interactions (Baselga and Araújo 2009), the analysis of these factors is complex and beyond the scope of this study.

This contribution represents an initial step towards ecological niche modeling based on plant communities for vegas in the Argentinean Central Andean Puna. The results obtained contribute to the knowledge of the environmental factors that limit the distribution of vegas dominated by different species assemblages, which give them particular characteristics in terms of productivity, organic matter accumulation and carbon storage capacity, among others. Knowledge of the current and potential distribution of these wetlands is a valuable tool, useful for planning in a region exposed to growing anthropogenic disturbances such as lithium mining and climate change.

Author contributions

E.C. performed the niche models and led the writing. A.E.I. conducted the field sampling and project administration. Both authors have contributed to the development of the ideas discussed in this article, the writing, and critically revised the manuscript.


This study was possible due to financial support from CONICET, Grant from PICT2018-04228. We thank all the colleagues who collaborated with field sampling and discussion of ideas during the last few years.


  • Aceituno P (1993) Elementos del clima en el altiplano sudamericano. Revista Geofísica 44: 37–55.
  • Ahumada M, Faúndez L (2009) Guía descriptiva de los sistemas vegetacionales azonales hídricos terrestres de la ecorregión altiplánica (SVAHT). Ministerio de Agricultura de Chile, Santiago de Chile, CL, 60 pp.
  • Badano EI, Cavieres LA (2006) Impacts of ecosystem engineers on community attributes: Effects of cushion plants at different elevations of the Chilean Andes. Diversity and Distributions 12: 388–396.
  • Barve N, Barve V, Jiménez-Valverde A, Lira-Noriega A, Maher SP, Peterson AT, Soberón J, Villalobos F (2011) The crucial role of the accessible area in ecological niche modeling and species distribution modeling. Ecological Modelling 222: 1810–1819.
  • Benfield AJ, Yu Z, Benavides JC (2021) Environmental controls over Holocene carbon accumulation in Distichia muscoides-dominated peatlands in the eastern Andes of Colombia. Quaternary Science Reviews 251: e106687.
  • Cabrera AL (1976) Regiones fitogeográficas argentinas. ACME, Buenos Aires, AR, 85 pp.
  • Carilla J, Grau HR, Paolini L, Morales M (2013) Lake fluctuations, plant productivity, and long-term variability in high-elevation tropical Andean ecosystems. Arctic, Antarctic, and Alpine Research 45: 179–189.
  • Carilla J, Grau A, Cuello AS (2018) Vegetación de la Puna argentina. In: Grau HR, Babot MJ, Izquierdo AE, Grau A (Eds) La Puna argentina: naturaleza y cultura [Serie de Conservación de la Naturaleza 24]. Fundación Miguel Lillo, Tucumán, AR, 143–156.
  • Casagranda E, Navarro CJ, Grau HR, Izquierdo AE (2019) Interannual lake fluctuations in the Argentine Puna: relationships with its associated peatlands and climate change. Regional Environmental Change 19: 1737–1750.
  • Domic AI, Capriles JM, Escobar-Torrez K, Santoro C, Maldonado A (2018) Two thousand years of land-use and vegetation evolution in the Andean highlands of northern Chile inferred from pollen and charcoal analyses. Quaternary 1: e32.
  • Domic AI, Capriles JM, Meneses RI, Pacheco P (2021) Plant community assembly is predicted by an environmental gradient in high-altitude wetlands in the semiarid western Bolivian Andes. Mires and Peat 27: 1–12.
  • Ferrier S, Drielsma M, Manion G, Watson G (2002) Extended statistical approaches to modelling spatial pattern in biodiversity in northeast New South Wales. II. Community-level modelling. Biodiversity & Conservation 11: 2309–2338.
  • Fick SE, Hijmans RJ (2017) WorldClim 2: new 1‐km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37: 4302–4315.
  • Gignac LD, Gauthier R, Rochefort L, Bubier J (2004) Distribution and habitat niches of 37 peatland Cyperaceae species across a broad geographic range in Canada. Canadian Journal of Botany 82: 1292–1313.
  • Guisan A, Theurillat JP, Kienast F (1998) Predicting the potential distribution of plant species in an alpine environment. Journal of Vegetation Science 9: 65–74.
  • Huntley B, Bartlein PJ, Prentice IC (1989) Climatic control of the distribution and abundance of beech (Fagus L.) in Europe and North America. Journal of Biogeography: 551–560.
  • Izquierdo AE, Grau HR, Carilla J, Casagranda E (2015a) Side effects of green technologies: the potential environmental costs of Lithium mining on high elevation Andean wetlands in the context of climate change. GLP news Newsletter of the Global Land Project 12: 53–56.
  • Izquierdo AE, Aragón R, Navarro CJ, Casagranda E (2018a) Humedales de la Puna: principales proveedores de servicios ecosistémicos de la región. In: Grau HR, Babot MJ, Izquierdo AE, Grau A (Eds) La Puna argentina: naturaleza y cultura [Serie de Conservación de la Naturaleza 24]. Fundación Miguel Lillo, Tucumán, AR, 96–111.
  • Izquierdo AE, Grau HR, Navarro CJ, Casagranda E, Castilla MC, Grau A (2018b) Highlands in transition: Urbanization, pastoralism, mining, tourism, and wildlife in the Argentinian Puna. Mountain Research and Development 38: 390–400.
  • Izquierdo AE, Carilla J, Nieto C, Osinaga-Acosta O, Martin E, Grau HR, Reynaga MC (2020) Multi-taxon patterns from high Andean peatlands: assessing climatic and landscape variables. Community Ecology 21: 317–332.
  • Izquierdo AE, Blundo C, Carilla J, Foguet J, Navarro CJ, Casagranda E, Chiappero MF, Vaieretti MV (2022) Floristic types of high‐Andean wetlands from northwest Argentina and their remote‐sensed characterization at a regional scale. Applied Vegetation Science 25: e12658.
  • Lobo JM, Herrero A, Zavala MA (2015) ¿Debemos fiarnos de los modelos de distribución de especies? In: Herrero A, Zavala MA (Eds) Los bosques y la biodiversidad frente al cambio climático: Impactos, vulnerabilidad y adaptación en España. Informe de Evaluación. Ministerio de Agricultura, Alimentación y Medio Ambiente, Madrid, ES, 407–417.
  • Maguire Jr B (1973) Niche response structure and the analytical potentials of its relationship to the habitat. The American Naturalist 107: 213–246.
  • Maldonado Fonkén M (2014) An introduction to the bofedales of the Peruvian High Andes. Mires and Peat 15: 1–13.
  • Martínez-Meyer E, Díaz-Porras D, Peterson AT, Yáñez-Arenas C (2013) Ecological niche structure and rangewide abundance patterns of species. Biology Letters 9: e20120637.
  • McKenzie NL, Belbin L, Margules CR, Keighery GJ (1989) Selecting representative reserve systems in remote areas: a case study in the Nullarbor region, Australia. Biological Conservation 50: 239–261.
  • Montesinos DB (2012) Vegetación halófila de tres localidades andinas en la vertiente pacífica del sur de Perú. Chloris chilensis: revista chilena de flora y vegetación 15(2).
  • Morales MS, Carilla J, Grau HR, Villalba R (2015) Multi-century lake area changes in the Southern Altiplano: a tree-ring-based reconstruction. Climate of the Past 11: 1821–1855.
  • Morales MS, Christie DA, Neukom R, Rojas F, Villalba R (2018) Variabilidad hidroclimática en el sur del Altiplano: pasado, presente y futuro. In: Grau HR, Babot MJ, Izquierdo AE, Grau A (Eds) La Puna argentina: naturaleza y cultura [Serie de Conservación de la Naturaleza 24]. Fundación Miguel Lillo, Tucumán, AR, 95–91.
  • Navarro CJ (2020) Respuesta funcional de las vegas de la Puna argentina a la interacción entre cambios climáticos y cambios de uso del suelo. PhD thesis, Universidad Nacional de Tucumán, Tucumán, AR.
  • Navarro CJ, Izquierdo AE, Aráoz E, Foguet J, Grau HR (2020) Rewilding of large herbivore communities in high elevation Puna: geographic segregation and no evidence of positive effects on peatland productivity. Regional Environmental Change 20: 1–11.
  • Navarro G, De la Barra N, Goitia E, Maldonado M (2011) Propuesta metodológica para la clasificación de los humedales altoandinos en Bolivia. Revista Boliviana de Ecología y Conservación Ambiental 29: 1–22.
  • Olson DM, Dinerstein E, Wikramanayake ED, Burgess ND, Powell GV, Underwood EC, D’amico JA, Itoua I, Strand HE, … Kassem KR (2001) Terrestrial ecoregions of the world: A new map of life on earth: A new global map of terrestrial ecoregions provides an innovative tool for conserving biodiversity. BioScience 51: 933–938.[0933:TEOTWA]2.0.CO;2
  • Ortíz-Yusty C, Restrepo A, Páez VP (2014) Distribución potencial de Podocnemis lewyana (Reptilia: Podocnemididae) y su posible fluctuación bajo escenarios de cambio climático global. Acta Biológica Colombiana 19: 471–481.
  • Osorio‐Olvera L, Lira‐Noriega A, Soberón J, Peterson AT, Falconi M, Contreras‐Díaz RG, Martínez-Meyer E, Barve V, Barve N (2020a) NTBOX: An R package with graphical user interface for modelling and evaluating multidimensional ecological niches. Methods in Ecology and Evolution 11: 1199–1206.
  • Osorio‐Olvera L, Yañez‐Arenas C, Martínez‐Meyer E, Peterson AT (2020b) Relationships between population densities and niche‐centroid distances in North American birds. Ecology Letters 23: 555–564.
  • Otto M, Scherer D, Richters J (2011) Hydrological differentiation and spatial distribution of high altitude wetlands in a semi-arid Andean region derived from satellite data. Hydrology and Earth System Sciences 15: 1713–1727.
  • Paoli H, Bianchi AR, Yañez CE, Volante JN, Fernández DR, Mattalía MC, Noé YE (2002) Recursos Hídricos de la Puna, Valles y Bolsones áridos del Noroeste Argentino. Convenio INTA EEA CIED, Salta, AR, 274 pp.
  • Peterson AT, Soberón J (2012) Species distribution modeling and ecological niche modeling: getting the concepts right. Natureza & Conservação 10: 102–107.
  • Polk MH, Young KR, Cano A, León B (2019) Vegetation of Andean wetlands (bofedales) in Huascarán National Park, Peru. Mires and Peat 24: 1–26.
  • Prentice IC, Bartlein PJ, Webb III T (1991) Vegetation and climate change in eastern North America since the last glacial maximum. Ecology 72: 2038–2056.
  • Prieto G, Alzérreca H, Laura J, Luna D, Laguna S (2003) Características y distribución de los bofedales en el ámbito boliviano del sistema T.D.P.S. In: Rocha OO, Saéz C (Eds) Uso pastoril en humedales altoandinos: Talleres de capacitación para el manejo integrado de los humedales de Argentina, Bolivia, Chile y Perú. Convención RAMSAR, WCS/ Bolivia, La Paz, BO, 13–40.
  • R Core Team (2020) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, AT. [accesed 28 Nov 2021]
  • Reboratti C (2005) Situación ambiental en las ecorregiones Puna y Altos Andes. In: Brown A, Martinez Ortiz U, Acerbi M, Corcuera J (Eds) La situación ambiental Argentina 2005. Fundación Vida Silvestre Argentina, Buenos Aires, AR, 33–39.
  • Ruthsatz B (1993) Flora und ökologische Bedingungen hochandiner Moore Chiles zwischen 18°00’ (Arica) und 40°30’ (Osorno) südl. Br. Phytocoenologia 23: 157–199.
  • Ruthsatz B, Schittek K, Backes B (2020) The vegetation of cushion peatlands in the Argentine Andes and changes in their floristic composition across a latitudinal gradient from 39°S to 22°S. Phytocoenologia 50: 249–278.
  • Soberón J, Peterson AT (2005) Interpretation of models of fundamental ecological niches and species’ distributional areas. Biodiversity Informatics 2: 1–10.
  • Troncoso C (2018) Valoración turística: tendencias recientes. In: Grau HR, Babot MJ, Izquierdo AE, Grau A (Eds) La Puna argentina: naturaleza y cultura [Serie de Conservación de la Naturaleza 24]. Fundación Miguel Lillo, Tucumán, AR, 426–440.
  • Van Aelst S, Rousseeuw P (2009) Minimum volume ellipsoid. Wiley Interdisciplinary Reviews: Computational Statistics 1: 71–82.
  • Wisz MS, Pottier J, Kissling WD, Pellissier L, Lenoir J, Damgaard CF, Dormann CF, Forchhammer MC, Grytnes JA, … Svenning JC (2013) The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling. Biological reviews of the Cambridge Philosophical Society 88: 15–30.

E-mail and ORCID

login to comment