Improving the knowledge of plant potential biodiversity-ecosystem services links using maps at the regional level in Southern Patagonia

Biodiversity supports multiple ecosystem services, whereas species loss endangers the provision of many services and affects ecosystem resilience and resistance capacity. The increase of remote sensing techniques allows to estimate biodiversity and ecosystem services supply at the landscape level in areas with low available data (e.g. Southern Patagonia). This paper evaluates the potential biodiversity and how it links with ecosystem services, based on vascular plant species across eight ecological areas. We also evaluated the habitat plant requirements and their relation with natural gradients. A total of 977 plots were used to develop habitat suitability maps based on an environmental niche factor analysis of 15 more important indicator species for each ecological area (n = 53 species) using 40 explanatory variables. Finally, these maps were combined into a single potential biodiversity map, which was linked with environmental variables and ecosystem services supply. For comparisons, data were extracted and compared through analyses of variance. The plant habitat requirements varied greatly among the different ecological areas, and it was possible to define groups according to its specialization and marginality indexes. The potential biodiversity map allowed us to detect coldspots in the western mountains and hotspots in southern and eastern areas. Higher biodiversity was associated to higher temperatures and normalized difference vegetation index, while lower biodiversity was related to elevation and rainfall. Potential biodiversity was closely associated with supporting and provisioning ecosystem services in shrublands and grasslands in the humid steppe, while the lowest values were related to cultural ecosystem services in Nothofagus forests. The present study showed that plant species present remarkable differences in spatial distributions and ecological requirements, being a useful proxy for potential biodiversity modelling. Potential biodiversity values change across ecological areas allowing to identify hotspots and coldspots, a useful tool for landscape management and conservation strategies. In addition, links with ecosystem services detect potential synergies and trade-offs, where areas with the lowest potential biodiversity are related to cultural ecosystem services (e.g. aesthetic values) and areas with the greatest potential biodiversity showed threats related to productive activities (e.g. livestock).


Introduction
Biodiversity is critical to support multiple ecosystem services (ES) (MEA 2005;Cardinale 2012;Mace et al. 2012). Several studies agree that plant biodiversity strongly affects supporting and regulating ES, e.g. soil nutrients cycling, productivity, and erosion control (Chen et al. 2018;Peri et al. 2019a;Ma et al. 2020). Furthermore, plant communities provide essential habitats for a wide range of species (Wen et al. 2020). Traditional management only focused on few provisioning ES, which is one of the main threats to biodiversity (Kohler et al. 2017). As a consequence, biodiversity loss endangers the provision of many other ES (MAES 2013; Shih 2020), affecting ecosystem resilience and resistance capacity to face climate change or species invasions (Kohler et al. 2017;Pires et al. 2018). Nevertheless, understanding these complex interactions represent a challenge, due to the multiple roles that biodiversity plays in the delivery of ES (MAES 2013; Mori et al. 2017).
Different factors affect species diversity distribution, where environmental heterogeneity (e.g. climate, soil, and topography) is the main driver for landscape variation in species composition across different taxa and biomes (Stein et al. 2014;Li et al. 2020). In consequence, plant biodiversity changes according to the different ecosystem types, e.g. temperate forests present unique understory species associated with tree canopy structure and composition (Barbier et al. 2008;Lencinas et al. 2008), while grasslands present characteristic species associated with temperature and rainfall gradients Li et al. 2020). In addition, soil characteristics as stock carbon and water balance, and topography as elevation and slope, influence plant distribution, richness, and diversity (e.g. Dubuis et al. 2013;Pottier et al. 2013;Peri et al. 2019a;Jordan et al. 2020). However, it is not totally understood how plant species distribution is linked with regional environmental variables (Kreft and Jetz 2007;Xu et al. 2016;Martínez Pastur et al. 2016a) and how biodiversity is related to ES supply (Schneiders et al. 2012;Martínez Pastur et al. 2017). Remote sensing techniques became a powerful tool for qualitative and quantitative modelling of plant species distribution and habitat prediction at the landscape level (Wiersma et al. 2011). One of the most employed methods is based on the combination of field plots with satellite data Pottier et al. 2013), where the Biomapper software showed its advantage to the development of potential habitability suitability (PHS) models (Hirzel et al. 2001(Hirzel et al. , 2006 in areas with scarce data as Southern Patagonia (Rosas et al. 2017). These models relate regional environmental characteristics (e.g. climatic, topographic, landscape) with the occurrence of species in one particular area (Guisan and Zimmermann 2000;Hirzel et al. 2002) determining the PHS of each species. During the last years, studies in plant species distribution were significantly increased (Bradley et al. 2012), mostly related to woody species (Silva et al. 2017), and invasive (Chai et al. 2016;Wan and Wang 2018) or endangered plant species (Abdelaal et al. 2020). Nevertheless, few studies combined multiple PHS to create a unique potential biodiversity map (PBM) that synthetizes the information for multiple species at the landscape level (e.g. Martínez Rosas et al. 2019a). PHS maps allow to understand the ecology of ecosystem communities (Martínez Pastur et al. 2016a), recognize environmental predictors of species distribution Pottier et al. 2013), identify potential hotspots of biodiversity (Wulff et al. 2013), and analyse the natural protected area network effectiveness and representativeness (Rosas et al. 2019a). However, few studies have linked the biodiversity and the ES supply to understand their relationship and support landscape management and conservation planning Shih 2020).
Southern Patagonia, specifically Santa Cruz Province, has a well-conserved wilderness landscape with a unique ecosystem variety, where temperature and rainfall patterns strongly defined the natural vegetation types. The region can be classified into eight main ecological areas, where dry steppes dominate in the northeast (central plateau, shrub-steppe of San Jorge Gulf, and mountains and plateaus), and shrublands and humid steppes (humid and dry Magellanic grass steppe) prevail in the southeast, while sub-Andean grasslands, Nothofagus forests, and alpine vegetation occupy a narrow strip in the west (Oliva et al. 2004) (Fig. 1). The vegetation types influence the provisioning ES delivery, generating several trade-offs between biodiversity and other ES supplies . Traditionally, steppe and shrubland areas present provisioning ES related to livestock production, where different management proposals have been implemented . Peri et al. (2013) also showed that different grazing strategies influence plant communities at the landscape level, where intensive management using mechanical shredding to remove shrubs significantly modified the arthropod community (Sola et al. 2016). Besides, extreme climatic conditions and overgrazing in dry steppes determined the highest values of desertification (Del Valle et al. 1998;Peri et al. 2016b;Gaitán et al. 2019). Recently, it was possible to map cultural ES , where higher values were related to cities close to protected areas and touristic routes . Forestry was conducted in the forests close to Los Andes; N. antarctica forests are related to livestock grazing in private lands (Peri et al. 2016c), while N. pumilio and evergreen forests are mainly under protected areas ). These studies also determined the impact of forest management alternatives on N. antarctica forests, where some multi-objective strategies (e.g. silvopastoral systems) allowed to maintain biodiversity values and increase ES delivery on time.
Regional plant diversity studies stand out for defining dominant vegetation types (Oliva et al. 2004;Oyarzabal et al. 2018) at the coarse resolution level and were used to identify potential negative trade-offs with different management activities at the local scales (Peri et al. , 2016c. The efforts to increase landscape ecology studies represent a challenge in remote areas such as Santa Cruz Province, where available data is scarce. Biodiversity mapping was studied through endangered species (Rosas et al. 2017), as well as lizards and darkling beetles , which represent good proxies for specific vegetation types. Martínez  proposed understory plants as proxies of potential biodiversity in Tierra del Fuego forests, because they have a wider distribution than other taxonomic groups, and Rosas et al. (2019a) also apply this methodology for Santa Cruz forests. However, these Fig. 1 Characterization of the study area. a Location of Argentina (dark grey) and Santa Cruz Province (black). b Mean annual temperature (8.6 to 13.5°C) where blue is low and red is high. c Mean annual rainfall (136 to 1681 mm year −1 ) where dark blue is higher values. d normalized difference vegetation index (NDVI) (0-1) where red is low and green is high. e Main ecological areas, where light grey is dry steppes (central plateau, shrub-steppes of San Jorge Gulf, and mountains and plateaus), grey is humid steppes (humid and dry Magellanic grass steppes), medium grey is shrublands, dark grey is sub-Andean grasslands, and black is forests and alpine vegetation (modified from Oliva et al. 2004) efforts for mapping biodiversity are not enough to compare the different vegetation types at the regional level, e.g. lizards and darkling beetles occurred in arid environments and not in forests. To understand the role that biodiversity plays in the ES delivery (MAES 2013), it is necessary to link biodiversity and ES supply at the landscape level, using the same proxy taxonomic group. With one unique map, it will be possible to analyse the potential trade-offs with other ES deliveries at the regional level (e.g. Martínez . In this context, the objective was to evaluate the potential biodiversity of the entire Santa Cruz Province (Argentina) and their linking with ES delivery, based on vascular plant species distribution as a proxy along the different ecological areas of the region. We also define the following specific objectives: (i) determine how much variation presents the habitat for indicator plants across the different ecological areas, (ii) relate PBM and environmental variables, and (iii) link PBM and ES supply. We hypothesized that (i) the distribution of the species is related to environmental gradients (Riesch et al. 2018), which together with other factors (e.g. biotic factors such as predation risk) define the potential biodiversity of one particular ecosystem; (ii) the provision of the different ES is not independent, interaction among them according to the species assemblage and biophysical characteristics of the ecosystems (Martínez Martínez Pastur et al. 2017); and (iii) human pressures on natural ecosystems due to economic activities (e.g. livestock use, silvopastoral management, forestry) generates synergies and trade-offs between provisioning ES and the other ES and biodiversity conservation (e.g. Peri et al. 2016b;Martínez Pastur et al. 2017). We expect to bring some evidences about these links between ES and biodiversity and provide tools for managers and decision-makers to propose more effective conservation alternatives at the regional level.

Study area
The study area includes Santa Cruz Province (Argentina) between 46°00′ and 52°30′ S, and between 66°00′ and 73°00′ W (Fig. 1a). The annual temperature increases from south to north (6.1 to 9.6°C, Fig. 1b) with a mean annual range from − 2.6 to 19.5°C (coldest and warmest months), while precipitation increases from east to west (204.1 to 899.7 mm year −1 , Fig. 1c) with a mean annual range from 13.7 to 30.3 mm month −1 (driest and wettest months). This environmental pattern strongly influences the normalized difference vegetation index (NDVI), where the highest values were detected in the west, medium values in the south, and lowest values in the northeast (Fig. 1d). As was described before, the province was classified in different ecological areas considering soil qualities, regional climate, and the dominant plant communities (modified from Oliva et al. 2004) (Fig. 1e): (i) Dry steppes (DS) are located in the north including central plateau steppes (CP), shrubsteppe San Jorge Gulf (SSJG), and steppes in mountains and plateaus (MP), covering an area of 156,764 km 2 . This area presents the lowest vegetation cover, mainly associated with small shrubs (e.g. Nassauvia glomerulosa) and grasses (e.g. Stipa spp.). (ii) Shrublands (SL) cover 28,373 km 2 and are dominated by Mulguraea tridens. (iii) Humid steppes (HS) are divided into humid and dry Magellanic grass steppes (HMGS and DMGS, respectively), covering an area of 17,831 km 2 , and are composed mainly by grasses (e.g. Festuca spp. and Carex spp.) in the south. Finally, other two ecological areas are located near to the Andes Mountains in the west, (iv) sub-Andean grasslands (SAG) that cover 8473 km 2 and are dominated by Festuca pallescens and (v) Nothofagus forests and alpine vegetation (FAV) that cover 6658 km 2 .

Database
A total of 53 plant species were selected for PHS modelling (Appendix 1). We used 977 plot surveys, where 235 plots belong to the PEBANPA network (Parcelas de Ecología y Biodiversidad de Ambientes Naturales en Patagonia Austral) , and 742 plots belong to FAMA laboratory (Forestal, Agrícola y Manejo del Agua of INTA EEA-Santa Cruz). This database includes vascular plant species, bare soil, and inferior plants (e.g. bryophytes and lichens) based on cover (%) and frequency of occurrence (%). We calculated a coveroccurrence index selecting the 15 most important plant species for each ecological area. Each species must present at least 20 presence points following the methodology described in Rosas et al. (2019a). Then, the database was complemented with the available data for the selected species in the Sistema Nacional de Datos Biológicos of Ministerio de Ciencia, Tecnología e Innovación Productiva (www.datosbiologicos.mincyt.gob.ar) (final presence points = 5914). Additionally, we explored 40 potential explanatory variables, which were rasterized at 90 × 90 m resolution using the nearest resampling technique on the ArcMap 10.0 software (ESRI 2011), where all variables (climate, topography, landscape) used in this study were previously described in Rosas et al. (2017Rosas et al. ( , 2018. Potential biodiversity was analysed considering the occurrence across the landscape and their links with the ES type supply through the comparison of the selected environmental variables in the PHS modelling (Appendix 2) and the main ecological areas. In addition, to link PBM values with ES supply (CICES V5.1, Haines-Young and Potschin 2018) we used different proxies to discuss potential trade-offs, according to available maps (grids) for the study area: (i) Supporting and regulating ES were related to soil quality and carbon sequestration through soil organic carbon and nitrogen , annual net primary productivity (Zhao and Running 2010), and desertification index (Del Valle et al. 1998) as a reverse proxy of soil control erosion. Soil organic carbon and nitrogen, as well desertification index, were modelled by related field data with climate and topographic variables. Annual net primary productivity was modelled using global productivity data of different ecosystems and remote sensing. (ii) Provisioning ES were related to food production through the probability of sheep occurrence (Pedrana et al. 2010), which were modelled according to climate and other bio-physical variables. Finally, (iii) cultural ES were related to spiritual, symbolic, and other interactions with the natural environment including aesthetic, existence, local identity, and recreational values Rosas et al. 2019b). The cultural ES (aesthetic value, existence value, recreation, and local identity) were mapped using geo-tagged digital images that local people and visitors posted on web platforms, which were related to social and biophysical landscape features.

Models and statistical analyses
For modelling, we follow previous regional proposals for forests in Tierra del Fuego (Martínez  and Santa Cruz (Rosas et al. 2019a) and methodological analyses of trade-offs to compare biodiversity and ES ) (see the chart flow presented in Fig. 2). PHS models for each plant species were based on the environmental niche factor analysis (ENFA) (Hirzel et al. 2002) using the Biomapper 4.0 software (Hirzel et al. 2004). ENFA employed the concept of ecological niche (Hutchinson 1957), comparing the ecogeographical variable distribution of presence dataset (locations of plant species) with the predictor distribution for the whole study area (Hirzel et al. 2001). ENFA summarizes all predictor variables and calculates two specific ecological relevance indexes: (i) Marginality (0 to 1) describes how far the optimum species habitat is from the mean environmental conditions of the study area, where higher values indicate that species grow in extreme environmental conditions. And (ii) global tolerance or specialization (tolerance −1 ) varied from 0 to infinite and describes the narrowness of the species niche, where higher values indicate that species tend to live in a narrow range of environmental conditions (Hirzel et al. 2002;Martínez Pastur et al. 2016a). We used a distance of geometric-mean algorithm to perform each PHS, which provides a good generalization of the niche (Hirzel and Arlettaz 2003). Each ENFA model results in a continuous PHS map with a range from 0 (minimum) to 100 (maximum habitat suitability). We used only the presence data in order to assess the statistical fit of each PHS model. ENFA used a cross-validation (Hirzel et al. 2006), evaluating the robustness and the predictive power of the models through Boyce index (B), continuous Boyce index (Bcont), proportion of validation points (P), absolute validation index (AVI), and contrast validation index (CVI) (Boyce et al. 2002;Hirzel and Arlettaz 2003;Hirzel et al. 2004Hirzel et al. , 2006. All these indices were described in Rosas et al. (2017Rosas et al. ( , 2018. We visualized PHS maps for each species into a geographical information system (GIS) with a resolution of 90 × 90 m, and we combined them with a mask based on NDVI < 0.05 to detect bare soil, ice fields, and water bodies (Lillesand and Kiefer 2000). Then, the 53 PHS maps were combined (average values for each pixel) to obtain the final PBM. This map had scores that varied from 1 to 46% (average values of PHS for all the studied species), and it was re-scaled by a lineal method from 0 to 100%. PBM was analysed using one-way ANOVAs. For this, we extracted the data from PBM using the hexagonal binning processes (each hexagonal = 250,000 ha) and calculated the average of potential biodiversity (PB) for each hexagon. After that, we classified each hexagon according to its PB values: (i) low was < 51% (32 hexagons), (ii) medium was 52-62% (44 hexagons), and (iii) high was > 63% (41 hexagons). The thresholds for each category were defined considering an equal number of pixels in the PBM. Finally, we analysed the extracted PB values considering the variables employed in the modelling, the main ecological areas, and ES proxies for each hexagon.
The species with the highest presence points was Osmorhiza chilensis (n = 357), while the lowest was Poa ligularis (n = 21) (Table 1). Three species presented a wide distribution and were related to six different ecological areas, e.g. Festuca pallescens was the most important species in SAG (frequency = 60%), P. spiciformis was the second important species in SL (frequency = 72%), and Pappostipa chrysophylla was the third important species in CP and MP (frequency = 61 and 60%, respectively). Another three species were related to five different ecological areas, e.g. Berberis microphylla was the second important species in SAG (frequency = 50%), while Carex andina was the second (frequency = 88%) and Nardophyllum bryoides was the third important species (frequency = 63%) in DMGS. Four species were related to four ecological areas, e.g. Festuca gracillima was the most important species in DMGS (frequency = 81%), as well as Nassauvia glomerulosa in CP (frequency = 72%) and N. ulicina in SSJG (frequency = 62%), and Senecio filaginoides was the third important species in SAG (frequency = 54%). Eight species (Acaena poeppigiana, Avenella flexuosa, Bromus setifolius, Carex argentina, Empetrum rubrum, Mulguraea tridens, Pappostipa ibarii, P. sorianoi) were related to three ecological areas, e.g. four presented high-ranking values (up to third): Avenella flexuosa was the most important species in HMGS (frequency = 59%), as well as Mulguraea tridens in SL (frequency = 52%) and Pappostipa sorianoi in MP (frequency = 70%), while Empetrum rubrum was the second important species in FAV (frequency = 33%). Twelve species (Acaena magellanica, Azorella prolifera, Chiliotrichum diffusum, Chuquiraga aurea, Clinopodium darwinii, Festuca magellanica, Galium aparine, Gaultheria mucronata, Hordeum comosum, Osmorhiza chilensis, Pappostipa chubutensis, Rytidosperma virescens) were related to two ecological areas, e.g. O. chilensis presented high-ranking value in FAV (frequency = 48%) and Ch. diffusum the third important species in  Table 1 Code for vascular plant species, number of presence point (N), mean cover (%), and frequency of occurrence (%) between brackets, followed by the rank (1 to 15) for the most important species of each ecological area (FAV forests and alpine vegetation, SAG sub-Andean grasslands, HMGS humid Magellanic grass steppes, DMGS = dry Magellanic grass steppes, SL shrublands, CP central plateau steppes, SSJG shrub-steppes of San Jorge Gulf, MP mountains and plateaus) (24) (24) 10 In the modelling of PHS maps, nine variables were the most appropriate to include in the analyses due to lower correlation values among all variables, and where only four variables presented higher values (> 0.80) based on Pearson's correlation index (Appendixes 3 and 4). The most useful variables in the 53 PHS-adjusted maps were annual mean temperature (AMT), annual precipitation (AP), and normalized difference vegetation index (NDVI). The correlation index among the nine selected variables varied between 0.03 and 1.00. The lowest correlation was detected between the minimum temperature of the coldest month (MINCM) and the distance to rivers (DR) (0.03), and the highest was between the maximum temperature of the warmest month (MAXWM) and the global potential evapotranspiration (EVTP) (1.00), AMT and MAXWM (0.98), AMT and EVTP (0.96), and MINCM and elevation (ELE) (− 0.85).
PHS maps presented remarkable differences in the spatial distribution (Appendix 7) and the habitat requirements (e.g. marginality and specialization) for each species (Fig. 3). Some species presented higher values of PHS in the steppes of the eastern areas of the province, showing the lowest specialization (1.60 to 7.54) and marginality (0.52 to 0.97) values. Among these species, nine (ADVO, CHAV, COSU, EPCH, FEAR, LYCH, MIRG, POLA, POLI) were related to a single ecological area, while four species (CLDA, PACH, CHAU, HOCO) were associated to two ecological areas, and ten species (ACPO, CAAR, MUTR, NABR, NAGL, NAUL, PACHR, PAIB, PASO, POSP) were associated to more ecological areas. Moreover, another species presented the highest PHS values in the western steppes near the ecotone, close to the mountain forests, showing low specialization values (2.06 to 6.76) and middle marginality values (1.35 to 2.16). Among these species, seven (AGPE, ARMA, BAMA, CAMA CAUN, HOPU, JUBA, PERE) were related to only one ecological area, while three species (AZPR, FEMA, RYVI) were associated to two ecological areas and six species (BEMI, BRSE, CAAN, FEGR, FEPA, SEFI) to more ecological areas. Finally, thirteen species showed the highest PHS values in a narrow area close to were related to only one ecological area, while five species (ACMA, CHDI, GAAP, GAMU, OSCH) were associated to two ecological areas and two species (AVFL, EMRU) to more ecological areas.

Potential biodiversity map of vascular plant species
The fitted PHS maps were combined into a single map, integrating all the selected vascular plant species (Fig. 4).
The PBM presented the lowest values (1-51%) near the mountain forests and ice fields in the west, and medium values (52-62%) in the dry steppes in the east. The highest PBM values (63-100%) occurred in three different steppe areas: (i) one important area of sixteen hexagons in the south (4.00 million ha), (ii) eleven hexagons in the north area of the province (2.75 million ha), and (iii) the smaller area with five hexagons located in the centraleast (1.25 million ha). The analyses of the modelled environmental variables allowed us to understand the occurrence of different PB levels in the landscape (Fig. 5). Temperature variables greatly influenced PB (AMT: F = 14.25, p < 0.001; MAXWM: F = 14.58, p < 0.001; MINCM: F = 8.69, p < 0.001), e.g. medium and high PB qualities occurred in areas where the mean annual temperature was higher (8.0°C) than the mean value of the whole province.
Annual precipitation negatively influenced PB (F = 14.24, p < 0.001), e.g. low PB quality occurred in areas where rainfall was above the mean value of the province (246 mm year −1 ). The higher rainfall was related to the upper mountains above the tree-line and closeness to ice fields. Elevation (F = 21.64, p < 0.001) followed the same pattern, while other landscape variables (distance to lakes and rivers) did not significantly influence PB qualities. Finally, PB qualities increased with the global potential evapotranspiration (F = 13.20, p < 0.001), while extreme values (low and high qualities) were presented in areas with higher NDVI values (F = 4.57, p = 0.012).
PB greatly changed through the main ecological areas (F = 19.93, p < 0.001), where shrublands (SL) presented the highest values followed by humid (HS) and dry steppes (DS). The lowest PB values were found in the sub-Andean grasslands (SAG) and in the mountain forests and alpine vegetation (FAV) (Fig. 6). Considering the supporting and regulating ES proxies, soil organic carbon stock (F = 2.73, p = 0.069), soil total nitrogen stock (F = 3.16, p = 0.046), annual net primary productivity (F = 4.00, p = 0.021), and desertification index (F = 9.61, p < 0.001) greatly varied according to PB qualities. The highest ES values were related to soil nutrient stock and carbon sequestration, which were mainly associated with low and high PB qualities. The desertification index (reverse proxy of soil control erosion) increased with PB qualities. The provisioning ES proxy, sheep presence probability, significantly changed across PB qualities (F = Fig. 3 Specialization (low species' variance compared to the global variance of all sites) vs. marginality (large difference of species' mean compared to the mean of all sites) of indicator vascular plant species, classified according to the number of times that the species were selected for each ecological area (blue = one ecological area, green = two ecological areas, orange = more than three ecological areas) 18.97, p < 0.001), increasing with PB quality. Finally, cultural ES proxies presented different connections with PB qualities: aesthetic (F = 6.84, p = 0.002), existence (F = 0.71, p = 0.492), and local identity (F = 3.14, p = 0.047) values decreased with PB qualities, and higher recreational values (F = 3.47; p = 0.035) were associated with higher PB qualities.

Potential habitat suitability of vascular plant species
The interest to understand species distribution at different scales (e.g. regional level) had been increased during the last years (Chen et al. 2018;Peri et al. 2019a). Nevertheless, the need of big data is a challenge in remote areas as Southern Patagonia (Martínez . In this context, PHS modelling based on ENFA is one feasible alternative (Hirzel et al. 2002) that links the presence data with environmental variables (Hirzel and Le Lay 2008), which was used for understory plant species in forests of Patagonia (Martínez Pastur et al. 2016a;Rosas et al. 2019a). In our study, data presence belonged from different database networks (e.g. local and national) covering all the landscape and the different ecological areas of the Santa Cruz Province and include the most relevant species of the local flora (e.g. 15 indicator plants for each ecological area).
The PHS maps support previous vegetation studies (e.g. Oliva et al. 2004;Oyarzabal et al. 2018), where different climate and environmental drivers strongly influenced plant species distribution (Fig. 1). Besides this, these studies also determine that some species are representative of multiple ecological areas (Table 1), e.g. Festuca pallescens is one of the most conspicuous and wide-dispersed species in Patagonia, whereas some autoecological studies (e.g. seed germination) associate this species with cold environments (López et al. 2019) and under 400 mm year −1 (Mancini et al. 2012). These findings support our results, where high PHS values were related to sub-Andean grasslands in the west. In addition, Mancini et al. (2012) indicated that shrub species (e.g. Empetrum rubrum and Gaultheria mucronata) grow close to Nothofagus forests in the west, while Gargaglione et al. (2014) showed that Avenella flexuosa were associated to N. antarctica forests in open environments (e.g. ecotone between grasslands and mountain forests). Despite the species associated to different ecological In the left graph, higher intensity colour represents greater potential biodiversity, while in the right graph, hexagons are classified in three levels of mean potential biodiversity: low potential = pale green (1-51%), medium potential = green (52-62%), and high potential = dark green (63-100%) areas, our ENFA results showed that their high marginality and specialization values indicated that these species tend to live in a narrow range of environmental conditions and where their optimum habitat differed from the mean environmental values of the province (Fig. 3). Rainfall gradient, which decreases from west to east, also influences vegetation distribution, e.g. Bertiller et al. (1995) identified in the extra-Andean steppes of Northern Patagonia that Festuca pallescens decreases while shrubs increases (e.g. Senecio filaginoides and Nassauvia glomerulosa) their cover along the rainfall gradient (from 600 to 170 mm year −1 ). In humid Magellanic grass steppes in the south, F. pallescens is also associated with Poa spiciformis and Carex andina (Peri and Lasagno 2010). However, in dry Magellanic grass steppes, F. gracillima is the dominant species accompanied by Nardophyllum bryoides and Nassauvia ulicina (Mancini et al. 2012). These distinct plant communities of steppes were mainly due to long-term differences in water availability (Peri et al. 2013). These studies support our results, where most of the species presented similar ecological requirements (e.g. low specialization and middle marginality values); however, N. bryoides, Nassauvia glomerulosa, and N. ulicina presented lower marginality values indicating that their optimum habitat is near to the mean environmental conditions of the province.
Other species presented a narrow distribution (e.g. two ecological areas) (Table 1), e.g. Acaena magellanica, Chiliotrichum diffusum, Galium aparine, and Osmorhiza chilensis. These species were identified by Huertas Herrera et al. (2018) growing in deciduous forests at lower elevations in Tierra del Fuego, Argentina. This is coincident with our results, where similar ecological requirements (e.g. high specialization and marginality values) and high PHS values were shown in forested areas and near humid steppes in the west. Other species presented F-value and probability (between brackets) are presented for each analysis, where different letters showed significant differences among categories of potential biodiversity using the Tukey test at p < 0.05. The numbers in the boxes indicate how many times the variable was used for modelling each plant species. AMT, annual mean temperature (°C); MAXWM, max temperature of the warmest month (°C); MINCM, min temperature of the coldest month (°C); AP, annual precipitation (mm year −1 ). In brown colour: ELE, elevation (m.a.s.l.); DLK, distance to lakes (km); DR, distance to rivers (km); EVTP, global potential evapotranspiration (mm year −1 ); NDVI, normalized difference vegetation index a more restricted habitat (e.g. Baccharis magellanica), showing a middle specialization and high marginality values. Gargaglione et al. (2014) associated this species to open environments, e.g. gap grasslands surrounding N. antarctica forests. Other species presented the highest marginality values (e.g. Blechnum penna-marina and Viola magellanica) and were related to evergreen forest in humid areas, close overstory canopies (Martínez Pastur et al. 2012), and other forested areas where some shrub species (e.g. Escallonia rubra) occurred (Mancini et al. 2012).
Some species with high specialization values (e.g. Berberis empetrifolia) were related to extreme xeric conditions near the Andean Mountains (Bottini et al. 2000), while others (e.g. Agrostis capillaris) were associated to ecotone areas between grasslands and N. antarctica forests (Alonso et al. 2019). López et al. (2019) found that Poa ligularis was associated to steppes in Northern Patagonia due to some of their ecological requirements (e.g. seed germination was related to a higher temperature than Festuca pallescens). In our study, this species was related to the steppes in the north-west of the province (mountains and plateaus), where high PHS values were related to the dry steppe ecological area. In addition, Chuquiraga avellanedae, Pappostipa chrysophylla, P. sorianoi, Nassauvia ulicina, Lycium chilense, and Festuca argentina have been also described in these extreme natural environments (Oliva et al. 2004;Peri et al. 2013;Oyarzabal et al. 2018), where ENFA analyses showed low specialization and marginality values, and PHS values associated to the dry steppes.

Potential biodiversity map of vascular plant species
Usually, PHS maps of plants species have been used to link species occurrence with environmental variables Pottier et al. 2013). Recently, new studies combine multiple maps to synthetizes the information of multiple species (Martínez Pastur et al. 2016a;Rosas et al. 2018), and particularly, Rosas et al. (2019a) create one PBM of understory species of forests in Santa Cruz. However, those studies focalized the analyses in particular ecosystems (e.g. forests), and not allowed the comparisons among distinct ecological areas. For example, most of the conservation efforts were focused on forests ), due to the lack of information about the importance of grasslands and shrublands for conservation in the Patagonia (Peri et al. 2021). This study considered the entire region Fig. 6 One-way ANOVAs comparing the following: a ecological areas-considering DS, central plateau, shrub-steppes of San Jorge Gulf, and mountains and plateaus; HS, humid and dry Magellanic grass steppes; SL, shrublands; SAG, sub-Andean grasslands; FAV, forests and alpine vegetation for the potential biodiversity as the main factor; b potential biodiversity (L < 51, M = 51-62, H > 62) for ES proxies as the main factors. The dotted line indicates the mean value for Santa Cruz Province; F-value and probability between brackets are presented for each analysis, where different letters showed significant differences among categories of potential biodiversity using the Tukey test at p < 0.05. Proxies for supporting ES were as follows: SOC, soil organic carbon (kg m −2 ); STN, soil total nitrogen (kg m −2 ); ANPP, annual net primary productivity (gr C m 2 year −1 ); and DES, desertification index (dimensionless). Proxy for provisioning ES was as follows: SPP, sheep presence probability (sheep presence probability km −2 ). Proxies for cultural ES were as follows: AES, aesthetic values (dimensionless); EXS, existence values (dimensionless); LI, local identity values (dimensionless); and REC, recreational values (dimensionless) and included the different ecosystem types across all the ecological areas (n = 8, from dry steppes to alpine vegetation) and a greater number of plant species (fifteen indicator species per ecological area). Different studies in Patagonia have focused on floristic heterogeneity description (e.g. physiognomy, dominant species, structural heterogeneity) (Paruelo et al. 2004), e.g. Oyarzabal et al. (2018) have re-described the main vegetation types of Patagonia, and Oliva et al. (2004) made an update for Santa Cruz Province. However, the available information is still scarce about plant biodiversity and their relationships with the regional environmental variables. One useful study characterized Patagonian ecosystem based on different approaches using climate as the main proxy, e.g. life zones or functional variables (e.g. NDVI) (Paruelo et al. 2001;Derguy et al. 2019). These approaches refer to vegetation physiognomy but not to the observed vegetation cover in the field, where more particular drivers (e.g. environmental heterogeneity) can influenced local biodiversity (Li et al. 2020). Our PHS maps were coincident with different vegetation cover studies (Oliva et al. 2004;Oyarzabal et al. 2018), where our PBM allowed to improve ecosystem characterization by capturing potential biodiversity changes according to the different indicator species. In a general trend, PBM values increased with temperature (south to north) and evapotranspiration values, where shrublands and steppes were predominant, while decreased with precipitation and elevation (west to east), where forests and alpine vegetation prevail ( Fig. 1). At the global scale, rainfall is strongly correlated with species richness along latitudinal gradients (Kreft and Jetz 2007), e.g. more richness was observed in rainy areas at lower tropical latitudes. In our study, low PB values were associated with higher precipitation values, while medium and high values of PB were related to low precipitation values. In our study area, the lowest biodiversity was related to dry steppes with the lower rainfall values located in the inner land. PB values increased across the rainfall gradient and closeness to seashores. However, below the close canopy of the mountain forest, the richness decreased despite the rainfall amounts, and in the upper mountains, the temperature becomes the limiting factor, influencing the number of plant species that can survive. Cleland et al. (2013) found a positive relationship between species richness and mean annual precipitation in grasslands of the USA. However, in Patagonia (46°00′-52°30′ S), biodiversity is more controlled by the availability of environmental temperature, where Kreft and Jetz (2007) found a significant positive effect for vascular plant species richness, supporting the described results.
Another challenge of this study was to link PB values with ES supply (e.g. Martínez ) that can help to understand the role that plant biodiversity plays in the delivery of these services (MAES 2013) at the regional scale. Biodiversity is assumed to be critical for ES supply (MEA 2005), where services can be obtained only if ecosystems include the species richness that guarantees the functional processes (Cardinale et al. 2011;Thompson et al. 2011;Mori et al. 2017). In fact, some authors indicated that biodiversity is itself an ES, because it is the basis for nature-based tourism or regulation of diseases (Mace et al. 2012). Supporting ES also had been demonstrated to present a strong relationship with plant biodiversity (Cardinale et al. 2011;Chen et al. 2018;Peri et al. 2019a). In this study, the proxies used for each ES type were soil organic carbon stock , soil total nitrogen stock , annual net primary productivity (Zhao and Running 2010), and desertification index (Del Valle et al. 1998). These proxies were previously proposed by several authors for similar studies (e.g. Stephens et al. 2015;Martínez Pastur et al. 2017;Chen et al. 2018;Wen et al. 2020) and were related with ecosystem functions and plant biodiversity considering different management scenarios (Wang et al. 2017), e.g. Peri et al. (2019a) found that soil carbon was a powerful predictor for threatened vascular plant biodiversity in Southern Patagonia. In our study, the higher values of ES proxies occurred at both extremes of PB qualities (low and high), e.g. low qualities were related to Nothofagus forests and high qualities to humid steppes and shrublands. Experimental studies showed that plant diversity increased their ES by improving soil microbial diversity and activity of grasslands (Fornara and Tilman 2008). Native Patagonian forests presented lower PB values, where high ES values were mainly related to the environmental conditions that improve organic residues and soil microbial biomass (Peri et al. 2016c. In this context, Peri et al. (2017) found that N. antarctica forests presented more soil organic carbon stock than other temperate forests. Provisioning ES were mainly related to regional economic activities (MEA 2005) and were represented here by sheep presence probability (Pedrana et al. 2010). This proxy increased with PB values, where humid steppes and shrublands prevail. Both areas support sheep and cattle livestock grazing , where some plant species (e.g. Festuca gracillima) constitutes about one-third of forage intake during winter (Oliva et al. 2005). Livestock grazing reduces plant species diversity, productivity, and vegetation cover and induces changes in soil structure and nutrient contents (Jiang et al. 2011;Dlamini et al. 2016). Peri et al. (2016b) indicate a negative effect on plant biodiversity due to continuous overgrazing in heterogeneous large paddocks in Patagonia. The steppe vulnerability is widely recognized (Gaitán et al. 2019), where intensive grazing decreases carbon and nitrogen stocks and increases desertification, generating losses in the ES supply (Del Valle et al. 1998;Peri et al. 2016b). The supply of cultural ES presented complex relationships with biodiversity, where particular vegetation types (e.g. forests) are closely related with high provision levels (Martínez  contributing to different dimensions of human well-being (Vilardy et al. 2011;Russell et al. 2013). In our study, aesthetic, existence and local identity were lower than recreational values at high PB qualities. Martínez  identified that most of the cultural ES in Santa Cruz were influenced by the closeness to water bodies and mountains (low PB values), while recreational values were influenced by the accessibility to cities and main routes in the steppes (low PB values) ). According to these results, provisioning ES can generate tradeoffs with the biodiversity conservation at the regional scale, while cultural ES can produce synergies with the biodiversity, allowing the creation of natural reserves (Rosas et al. 2019a).
Several human-related activities in Southern Patagonia (e.g. livestock) negatively influenced the original ecosystems by modifying the plant biodiversity, soil properties, and structure (Peri et al. 2013(Peri et al. , 2016c. At the regional level, most of the natural reserves are located close to the Andes Mountains, protecting mainly native forests (N. pumilio and N. betuloides) and specific ES related to tourism (e.g. aesthetic ES) (Martínez . However, N. antarctica forests presented the highest PB values (Rosas et al. 2019a), and less than 20% of these forests are under formal protection (e.g. National Parks or Provincial Reserves). This highlights the importance of conservation strategies on private lands (e.g. land-sharing strategy), where most of the economic activities influence plant biodiversity ). There is an increasing interest to integrate ES arguments within management plans, where new strategies (e.g. silvopastoral systems) have been developed to reconciles the productive needs and requirements of local people with biodiversity conservation (Peri et al. 2016c). In the other extreme, Patagonian steppes presented lower natural reserves than any other ecological region (e.g. Bosques Petrificados Natural Monument and Monte León National Park), and where most of them are located near the sea coast. Steppe areas presented more complex interactions with the environment and ES; however, a general trend showed that high PB values are coincident with high values of different ES supplies (provisioning, regulating, and supporting). For more than 100 years, livestock grazing influenced and modelled the region (Peri 2011) in large paddocks that varied from 1000 to 20,000 ha (Ormaechea and Peri 2015), increasing the desertification processes (Del Valle et al. 1998). Bjerring et al. (2020) experimentally demonstrated that the adjustment of the stocking rate or rotational grazing had a more positive effect on plant biodiversity and biomass growth than continuous grazing in the steppes. Some grazing managements, as warm-season grazing, improve ES delivery and biodiversity conservation in temperate steppes of northern China (Wang et al. 2017). This sustainable management reduces grazing intensity using specific areas during summer and autumn, and it is considered as an appropriate management strategy for maintaining plant species diversity and soil texture (Cingolani et al. 2005;Bjerring et al. 2020). In addition, this type of management is more efficient when some paddocks were excluded (throughout the year by fencing) due to degradation of ES delivery and biodiversity in the temperate steppes (Wang et al. 2016). To date, there are no management strategies that protect the biodiversity (e.g. land sharing) (Fischer et al. 2014), where shrublands are identified as an obstacle for livestock, e.g. some management proposals promote the removal of the shrubs to increase grass cover (Sola et al. 2016). It is necessary to reconcile these management proposals and the maintenance of the biodiversity in the hotspot areas identified in this work. The regional approach of PBM of the most relevant plants for each regional vegetation type allows us to improve the knowledge about plant potential habitat suitability and supports other studies about vegetation description at the regional level. Moreover, PBM allowed us to relate plant biodiversity and ES supply and help us to understand their relationships to identify potential negative impacts and trade-offs with economic activities.

Conclusions
Potential habitat suitability models using the Biomapper software allowed us to develop maps of relevant indicator vascular plant species at the landscape level using the presence data only and environmental variables (e.g. climate, topography) free availability on the public internet. ENFA indexes (marginality and specialization) allowed us to characterize plant habitat requirements and link them to environmental gradients of the landscape, helping in categorizing the plants according to its distribution in the landscape. The final map of potential biodiversity synthetized the information of the different plant species and assisted to identify hotspot areas across the province (e.g. mainly in humid steppes and shrublands). These areas were outside of the formal reserve network, and they are not included in the conservation strategies implemented in the region. In addition, the link between potential biodiversity and ES allowed us to understand potential trade-offs and synergies across the landscape and for the different ecological areas (e.g. from dry steppes to alpine vegetation). Areas with the greatest biodiversity values also showed high provisioning ES, and this spatial congruence determinates potential trade-off with conservation and other ES supply. The proposals explored in this paper is an effective tool for areas with low availability of databases to detect biodiversity hotspots and potential trade-offs that can be rapidly incorporated into decision-making to raise better spatial planning for plant biodiversity conservation.    Hijmans et al. (2005) (2) Farr et al. (2007) ( 3) SIT -Santa Cruz (http://www.sitsantacruz.gob.ar) Consortium for Spatial Information (CSI) (Zomer et al. 2008) ORNL DAAC (2008) Zhao and Running (2010) VIMA X X X X X Appendix 4