Ecological niche modeling, niche overlap, and good old Rabinowitz’s rarities applied to the conservation of gymnosperms in a global biodiversity hotspot

Biodiversity hotspots harbor 77% of endemic plant species. Patagonian Temperate Forest (PTF) is a part of a biodiversity hotspot, but over the past centuries, has been over-exploited, fragmented and replaced with exotic species plantations, lately also threatened by climate change. Our aim is to better understand patterns of habitat suitability and niche overlap of nine endemic gymnosperm species, key elements of the PTF, complementing traditional approaches of biodiversity conservation. Using R packages and 3016 occurrence data, we deployed ecological niche models (ENM) in MaxEnt via kuenm, and classified species according to Rabinowitz’s types of rarity. We then overlapped their niches calculating Schoener's D index, and considered types of rarity in a spatial ecological context. Finally, we overlay high species’ suitability and protected areas and detected conservation priorities using GapAnalysis. We generated simplified ENMs for nine Patagonian gymnosperms and found that most niches overlap, and only one species displayed a unique niche. Surprisingly, we found that three species have divergent suitability of habitats across the landscape and not related with previously published geographic structure of neutral genetic variation. We showed that the rarer a species is the smaller niche volume tend to have, that six out of nine studied species have high conservation priority, and that there are conservation gaps in the PTF. Our approach showed that there are unprotected suitable areas for native key species at high risk in PTF. Suggesting that integrating habitat-suitability models of multiple species, types of rarity, and niche overlap, can be a handy tool to identify potential conservation areas in global biodiversity hotspots.


3
Vol:. (1234567890) species via high speciation rates and/or preservation of species over time, via low extinction rates (Jetz et al. 2004), usually hosting ancient and young species (Harrison and Noss 2017). Patagonian Temperate Forest (PTF) runs through the spine of the Southern Andes in Argentina and Chile (Armesto et al. 2001;Smith-Ramirez 2004, Fajardo et al 2021. This region includes part of a global biodiversity hotspot in Central Chile (Myers et al. 2000), and is also of conservation concern as part of the World Wildlife Fund's, denominated Valdivian Ecoregion (www. world wildl ife. org/ ecore gions/ nt0404). However, over the past centuries, has been over-exploited, fragmented and replaced with exotic species plantations (Miranda et al. 2016), and is lately also threatened by climate change. This forest is frequently dominated by an ancient group of phylogenetically related gymnosperm species, whose gene pool has shown its ability to keep pace with past climatic changes. Although some of them have been focus of ecological and genetic studies, either for their cultural or economic value, others are relatively understudied, due to their remote distribution and inaccessibility (Premoli et al. 2000(Premoli et al. , 2002Allnutt et al. 2001;Donoso 2006;Quiroga and Premoli 2010;Souto et al. 2015a). Since many of its native species are considered rarities, included in the IUCN red list, urgent conservation actions are in need (Smith-Ramirez et al. 2004;Alaniz et al. 2016;Miranda et al. 2016), moreover as sometimes, rare species are the most prone to extinction risk (Gaston 1994;Gaston et al. 2000;Harnik et al. 2012). Given limited resources, conservation efforts should arguably be directed toward taxa and/ or areas that represent the greatest amount of unique evolutionary history (Scherson et al. 2014). Currently new aims in conservation contemplate large geographic areas and/or long-time scales, to preserve evolutionary potential. Highly diverse groups that share a set of ecological adaptations along an environmental gradient, are exemplar opportunities to test and develop specific actions and target areas aimed to protect a biome. Such studies might shed light to understand ecological niche differences and handy to develop conservation and mitigation actions against global change impacts (Aguirre- Gutiérrez et al. 2015).
Current approaches to the conservation of biodiversity hotspots include the application of biogeographical principles, theories, and analyses, concerned with range shifting of taxa (Whittaker et al. 2005). Species distribution models (SDMs), which describe both, a species' niche and habitat suitability that supports it (Richardson and Whittaker 2010), and ecological niche models (ENM) that identify potential areas of distribution (Soberón and Nakamura 2009), constitute essential tools in conservation biogeography (Franklin 2010;Scoble and Lowe 2010). In these studies, the niche conservatism concept allows inference of species responses to climate changes in the past and how they might shift species' distribution, especially regarding changes in their habitat specificity (Wiens et al. 2010). Also, similarities and differences in niche characteristics among populations within a species can be statistically analyzed using ENM (Wiens and Graham 2005). Advances in studies on how global climate change affects species have highlighted the importance of knowing the differences between species in terms of their environmental requirements in a geographical context and in an extension comparable to that of the species' ranges (Warren et al. 2008). ENMs are correlative models unrelated to the biological processes underlying genetic differentiation and adaptation of a species to its environment (Morente-López et al. 2020), so they can identify areas with different levels of environmental suitability for a species within its range. In this way, ENMs allow a user not only to estimate the distribution of a species and identify putative important environmental variables that shape it (Searcy and Shaffer 2016), but also to postulate hypotheses on divergent selection within a species' range (Blanquart et al. 2013). ENMs are a valid tool to identify potential intraspecific differences, which may involve historical events of reproductive isolation, and or divergence in ecological tolerances that could be evidence of potential speciation events (Silva et al. 2014;Zhao et al. 2019Zhao et al. , 2021. The niche overlap framework has been shown as a robust statistical method to test differences between-species niches', as well as subspecies or populations of the same species which are therefore likely to experience different climatic conditions (Broennimann et al. 2012).
Species conservation depends on the preservation of their habitats (Dengler et al. 2008), and its goals have evolved over time from a more speciescentric approach to one that considers keystone species and their role in the ecosystem (Lacher 2018), including rare ones. Types of rarity theory, states that range size, habitat specificity, and local population size, are environmental constrains that generate many ways in which a species can be rare (Rabinowitz 1981). There are eight forms of rarity, defined according to three species' characteristics: 1. Range size, so species can be widespread or narrow; 2. Habitat specificity, classifying species in broad or restricted, and 3. Local population size, then species have, either large, somewhere dominant, or small, non-dominant populations. From the combination of these three two-state variables, the eight possible types of rarity emerge. Current, online abundant occurrence records could potentially be used to define plant species types of rarity and their level of endemism, adding a new tool to the design of conservation plans (Fleishman et al. 2006;Choe et al. 2019). Geographic ranges of closely related species can vary dramatically, yet we do not fully grasp the mechanisms underlying such variation. The niche volume hypothesis posits that species that have evolved broad environmental tolerances can achieve larger geographic ranges than species with narrow environmental tolerances (Slatyer et al. 2013). This assumption has been used in multiple studies of conservation and distribution of dissimilar lineages (Yu and Dobson 2000;Choe et al. 2019;Ramírez-Bautista et al. 2020), giving a conceptual framework to make multiple comparisons.
In this study we combine the background presented above to deploy patterns of habitat suitability and niche overlap of nine gymnosperm species, key elements of the PTF biodiversity hotspot, complementing traditional approaches of biodiversity conservation. Particularly, we aim to: (1) Define ENM for focal species; (2) Test pairwise niche overlap to detect shared environmental constrains among focal gymnosperms, and within a species' range if indication of divergent climatic conditions were observed; (3) Classify each species according to Rabinowitz's types of rarity and test niche volume hypothesis that widespread species with low habitat specificity, and continuous populations have larger niche volume compared to narrow species, with high habitat specificity and scattered populations; finally (4) Combining ENM, niche overlap and forms of rarity, to pinpoint and prioritize conservation gaps in a global biodiversity hotspot.

Materials and methods
The Patagonian Temperate Forest have a high number of endemism, making it a global biodiversity hotspot (Smith-Ramírez 2004;Scherson et al. 2014), with sectors of high genetic diversity evidencing a complex evolutionary history (Souto et al. 2015b;Scherson et al. 2017). Product of these events, PTF is considered a biogeographic island, arisen along several geological/climatic events (Armesto et al. 1998), with a dynamic history, undergoing almost continuous changes since the Paleocene, stablishing the current species' associations (Hinojosa and Villagrán 1997;Barreda and Palazzesi 2007).

Studied species and presence records
We focus our study in an endemic phylogenetically related group of species of biogeographical interest, the Coniferophyta II (here after gymnosperms) from PTF. This group includes nine species from three families: Araucariaceae, with one species (Araucaria araucana (Molina) K. Koch); Cupressaceae, with three monotypic genera (Austrocedrus chilensis (D. Don) Pic.Serm. & Bizzarri, Fitzroya cupressoides (Molina) I.M. Johnst. and Pilgerodendron uviferum (D. Don) Florin); and Podocarpaceae with five species (Lepidothamnus fonkii Phil., Podocarpus nubigenus Lindl., Podocarpus salignus D. Don, Prumnopitys andina (Poepp. ex Endl.) de Laub. and Saxegothaea conspicua Lindl.). These are all ancient taxa that have long successfully navigated climate and geology, been confined to a narrow strip of land in Southwest South America. Araucariaceae and Podocarpaceae diverged between 176 and 230 Ma ago (Leslie et al. 2012;Quiroga et al. 2016), while Austrocedrus, Pilgerodendron, and Fitzroya diverged circa 120 Ma ago (Yang et al. 2012). All these species are emblematic either for being taxonomic rarities, relict lineages, long-lived or dominant species. Other authors have measured evolutionary distinctiveness and scarcity to inform conservation at biogeographical scales, and suggested that lineages with deep time of origin in these lineages and their phylogenetic uniqueness are elements to prioritize the preservation of evolutionary diversity, which is one way to maximize genotypic and functional diversity (Cadotte and Davies 2010). The Cupressaceae and Saxegothaea are mono-specific and unique in the world. The other species have many related taxa within the continent, as is the case of Araucaria, Prumnopitys and Podocarpus, as well as in other former Gondwanan continents.

Climatic variables and ecological niche models
We completed our personal data collection of focal species' presence records, with available occurrence sources: GBIF repository, Concepcion (Chile) regional herbarium (CONC), and CONAF territorial information system (http:// sit. conaf. cl/). Uninformative or ambiguous localities of herbarium labels for each species, were thoroughly reviewed by comparing the presence according to CONAF maps (http:// sit. conaf. cl/). Also our field knowledge was used as a criterion to evaluate some uninformative or ambiguous records. We checked each datum before including it in our data base and cleaned all species' data sets with a minimum distance of 30 s among records using NicheToolBox (NTB) R-script (Osorio-Olvera et al. 2020). After removing duplicates and neighbors, we obtained a total of 3016 presence records (Table 1).
In order to define the fundamental niche of the nine gymnosperms of the PTF, we clipped 19 bioclimatic variables (BIO) from the present-day World-Clim database at 30 s resolution   (Fick and Hijmans 2017; http:// www. world clim. org), representing summaries of means and variation in temperature and precipitation. Since, each species has its ecological tolerances, we defined species-specific M zones according to their potential distribution range as suggested by Soberón and Nakamura (2009). Each M zone includes regions close to distribution records of each species and, according with current environmental conditions, that could be possible areas of distribution under an assumption of natural migration from current localities within the boundaries of PTF. Even among some of them there is no report of sympatry, as a result of independent evolution. Each M zone includes regions close to distribution records of each species and, according with current environmental conditions, could be possible areas of distribution under an assumption of natural migration from current localities within the boundaries of PTF. For each species, in its respective M zone, we performed suitability models using MaxEnt v3.3.3 (Phillips et al. 2017), that iteratively maximizes the likelihood of the occurrence dataset using one or a combination of: linear, quadratic, and product with or without threshold and hinge features, recording in all cases the logistic output (Elith et al. 2011). We used MaxEnt because is an algorithm that can be use with background pseudoabsences, since true absences were not available for all species (Elith and Graham 2009). To estimate bioclimatic variable's relevance in each model, we run a jackknife procedure. We cross validate each model with 10 replicates, 500 iterations, 10 percentile training presence, and regularization multiplier of 1. We analyzed logistic distribution of bioclimatic variables at species level in order identified no Gaussian distributions. Species with this anomalous variable distribution, were partitioned in two groups and analyzed using the same M zone, that include total distribution. We used the same M zone to evaluate divergent bioclimatic conditions, in order to fully evaluate the entire species, and to determine potential areas of suitability in the other sector. Exceptionally, we doubled the number of iterations for L. fonkii, an extremely rare species, with only 37 available records. Once we obtained each species' general model in MaxEnt, we applied dimension-reduction techniques to minimize strong covariance among environmental variables and to avoid over fitted models that reduce accuracy (Synes and Bsborne 2011;Boria et al. 2014), using kuenm R-script (Cobos et al. 2019), that automates calibration and evaluation steps in ENMs. We performed Pearson's correlations among the 19 bioclimatic variables for each species, since each one has its own M zone. For each species we compared alternative models and sets of variables considering (1) statistical significance, (2) predictive power, and (3) model complexity. We used jackknifing procedures in MaxEnt and bioclimatic variables correlation matrixes to build alternative models for each species selecting different sets of variables that contributed the most to each model. To allow broader number of sets of variables to compete in kuenm, we set a higher than usual threshold (|0.85|). Each comparison included different sets of variables, with multipliers (0.5, 1, 1.5, 2, 4, 6, 8) and features (lineal, quadratic, product, hinges and threshold). For each species we test all possible models until obtaining a statistically significant simplified ecological niche model (S3ENM). This model calibration is a Bayesian procedure that finds the combination of variables and parameters that best represents the phenomenon of interest fitting with the data.

Niche overlap
We estimated niche overlap, to test shared environmental constrains between pairs of species or between two different sets of populations within a species' range when they showed different niches across their range as an indication of divergent selection. For each pairwise comparison we run scripts (downloaded in 2019) provided in Broennimann et al. (2012), including only bioclimatic variables selected in S3ENM for each species, and both species occurrences were pooled and randomly split into two datasets, maintaining the number of occurrences as in the original dataset. The analysis of niche overlap involves three steps: (1) calculation of the density of occurrences and of environmental factors along the environmental axes of a multivariate analysis (2) measurement of niche overlap along the gradients of this multivariate analysis and (3) statistical tests of niche equivalency and similarity (Warren et al. 2008). Then the niche overlap statistic Schoener's D was calculated, which shows values that range from 0 (no overlap) to 1 (complete overlap) (Schoener 1968). This method uses cell-by cell comparisons of geographical predictions of occurrences and randomization tests to quantify niche differences and assess their statistical significance (Warren et al. 2008), with 100 repetitions, to ensure that the null hypothesis can be rejected with high confidence, building a histogram of simulated values. If the observed D value falls within the 95% density of simulated values, the null hypothesis of niche equivalency cannot be rejected. We also performed a PCA, applying to this multivariate approach a smoothers kernel to densities of species occurrence in gridded environmental space, in order to calculate metrics of niche overlap (Broennimann et al. 2012).

Rabinowitz' rarities and niche volume
In this study we followed Rabinowitz's (1981) classification of rarities which is constructed with qualitative and flexible categories according to geographic range, habitat specificity, and local populations size (Supplementary information SI4). A species was considered "widespread" if its range spans more than 8 degrees of latitude, while a species was "narrow" when its range is 8 degrees of latitude or less. Regarding habitat specificity, we assigned each species either to "broad", when it occurs in several habitats, or to "restricted" if it shows high fidelity to its habitat. Finally, we classified species according to its local population size and continuity, in "large-dominant" or "small-isolated". These classifications were made following species' descriptions available in the literature along with our personal observations in the field (Donoso 1995(Donoso , 2006Luebert and Pliscoff 2006).
We also estimate each species' Euclidian niche volume in NTB R-package using bioclimatic variables selected in S3ENM, except for L. fonkii north and south populations, where we included the three variables with higher percent of contribution in Max-Ent models, performed without hinges or threshold features. For each species, we ran 10 replicas of Euclidian niche volume calculation. Then, we tested differences in niche volume between rarity types via lineal models, with range width, habitat specificity, and population size as explanatory variables. We also assigned an increasing rarity score to each of the eight cells of the Rabinowitz's classification of species: 1→ common (wide, broad, dominant), 2→ predictable, 3→ sparse, 4→ discontinuous, 5→ unlike, 6→ endemic, 7→ rare, 8→ maximum rarity (narrow, restricted, small/isolated). With these values, we tested if increasing rarity is related with niche volume, deploying lineal models using niche volume as response variable performed in lmer package R 3.4.4 (R Development Core Team 2017).

Conservation priorities and gap areas
We produce metrics that estimate the degree of representation of each species in protected natural areas and identified their conservation priorities applying ecogeographic methods and tools using GapAnalysis R package (Carver et al. 2021). We performed the option in situ using our occurrence records data set for each species, each species' S3ENM, a moderate suitability value of 0.5, and WDPA protected areas layer (UNEP-WCMC 2021 retrieved from www. prote ctedp lanet. net). This package calculates four scores for each species: a sampling representativeness score (SRS), a geographical representativeness score (GRS), an ecological representativeness score (ERS), and, a final conservation score (FCS), averaging all previous scores. All scores are bound between 0 and 100, with 0 representing an extremely poor state of conservation, and 100 representing comprehensive protection. Also, to graphically identify potential conservation gap areas, we overlaid suitability areas with values above 0.5 for all species together and regional protected areas layer in QGIS Development Team (2019). Considering this conservation gap map, the degree of representation of each species in protected natural areas, ENMs, niche overlap results and forms of rarity, we discuss an integrated prioritization scheme.

Niche modelling
Complete models for the nine species including the 19 bioclimatic variables were run, including or not Hinges and Threshold features, with ten replicates and a multiplier of 1, to obtain contribution and permutation values, logistic distribution and correlation matrices. After a detailed analysis of these results, we proceeded to assemble different sets comprising a reduced number of low correlated bioclimatic variables, with high contribution and permutation in the complete model, considering also if the logistic distribution of each variable was adequate (each model is in Supplementary Information SI1, including the AUC values, correlations tables are in Supplementary Information SI2). Significant reduce models included variables with correlation values of |0.75| or less, with the maximum value of r = 0.714 (In supplementary information SI2 find Pearson's correlation values for the bioclimatic variables along the M of each species). We obtained S3ENM including 3 bioclimatic variables (Fig. 1, Supplementary Information SI1 for detailed model results) that represent a possible set of low correlated variables explaining in a simplified way the potential distribution of each species. Precipitation seasonality (BIO15) was a recurrent environmental constrain, contributing between 21 and 41% in four out of the six models deployed ( Supplementary Information SI1, SI3). The three remaining species delivered low AUC values, when using the full set of bioclimatic variables in Max-Ent. In these species we observed bimodal logistic distributions: for P. uviferum in BIO15; L. fonkii in BIO1, BIO2, BIO17 and BIO18; S. conspicua in BIO3 and BIO15; Supplementary Information SI3). Pilgerodendron uviferum and L. fonkii showed ENMs with relatively low AUC values (0.848 and 0.721 respectively) and S. conspicua, although displayed an AUC value of 0.917, the percentage contribution and permutation importance of climatic variables exhibited discrepant response curves. In order to overstep observed inconsistencies, and based in bimodal gaussian distribution, we split the records of these three species into two geographic sectors each. For P. uviferum and L. fonkii, we split the records in a latitudinal fashion, separating north and south of − 44.5°L and − 47°L, respectively. And for S. conspicua we split the records longitudinally, distinguishing between coastal (west) and Andean (east) localities (Fig. 1). For each sector, we obtained ecological niche models with higher AUC values than those obtained with complete sets of records, and also more reliable suitability maps ( Supplementary Information SI3). We identified different climatic drivers within each sector, as well as accurate response curves, and jackknifing values ( Supplementary Information SI3). We assume these results as two divergent ecological sectors for these three species, and tested specific S3ENMs for each sector, using the same background zone to model both sectors, and the projection of one sector did not predict the species' presence in the other (Supplementary Information SI1). Once we partitioned L. fonkii data, the reduced model analysis test lack convergence, due to low number of records, so for further analysis we retain the three bioclimatic variables with higher percent of contribution in MaxEnt models performed without hinges and threshold features ( Supplementary Information SI3). Precipitation seasonality (BIO15) was also a recurrent environmental constrain, contributing between 24.4 and 51.4% in five out of these six models deployed (two

Niche overlap
We analyzed 65 pairs of combinations (Table 2), representing the niche of the species along the two first axes of the PCA for each pair [only shown four possible cases Fig. 2(1)]. We estimate the contribution of the climatic variables on the two axes of the PCA and the percentage of inertia explained by the two axes (data not shown). We pairwise compared the overlap of S3ENMs for all focal species, and found that 17% of the niches overlapped, with eleven statistically significant pairs (Table 2, Fig. 2). Araucaria araucana niche showed no overlap with any other studied species (Table 2), even when it is sympatric with three of them (A. chilensis, P. andina, and P. salignus). We found that the following pairs occupy the same portion of the niche space: P. andina and A. chilensis (Fig. 2a), F. cupressoides and P. andina (Fig. 2b), F. cupressoides and S. conspicua east, and also the northern sectors of P. uviferum and L. fonkii (Table 2, Fig. 2a-a′, c-c′). Many of the studied species are partially sympatric, however their niches do not overlap, this is the case of A. araucana and P. andina (Fig. 2c). In addition, P. uviferum north overlaps niche with P. nubigenus, P. salignus and S. conspicua west; while P. uviferum south overlaps niche with L. fonkii south and P. nubigenus, and L. fonkii south niche is overlapped with P. nubigenus. Finally, P. andina and S. conspicua east have overlapped niches. For the three species showing two divergent ecological sectors across their range as S. conspicua east and west, P. uviferum north and south and L. fonkii north and south (Fig. 2d-d′); pairwise niche equivalence test yielded no niche overlap between sectors ( Table 2).

Comparison of species considering Rabinowitz' rarities
We classified the nine focal species of this study following Rabinowitz's types of rarity (Table 3), and tested the difference in niche volume between rarity categories (Supplementary Information SI4). We found statistically significant differences in niche volume between species with wide and narrow range, with widespread species showing bigger niche volume (F = 16.55; p < 0.001) than narrow species. We also found statistical differences in niche volume between species with broad and restricted habitat specificity (F = 4.54; p = 0.03), showing that species with restricted habitat specificity present smaller niche volume. Also, we observed significantly larger niche volume in species with large dominant populations than species with small isolated populations (F = 18.38; p < 0.001). Finally, we tested the relationship between niche volume and types of rarity, finding that niche volume significantly decreases with increasing rarity (t = − 7.48, p < 0.01, r 2 = 0.32).

Conservation priorities and gap areas
Using metrics that estimate the degree of representation of each species in protected natural areas we found that species with high conservation priority (HP) are F. cupressoides, P. uviferum north, L. fonkii north, P. nubigenus, P. salignus and S. conspicua west, with low final conservation scores (FCS) ranging from 0.46 to 24.25. We obtained intermediate scores of conservation priority for P. andina (FCS = 41.20), and A. araucana, A. chilensis and P. uviferum south presented low conservation priority (50 < FCS < 75). Meanwhile, L. fonkii south and S. conspicua east are sufficiently conserved (SC) with FCS > 75 (Table 4). We identified areas of high suitability among nine native gymnosperm species of PTF (Fig. 3). One conservation gap area is located between − 35° to − 40°S latitude in the Costal Range, where the P. salignus has the higher suitability values (Fig. 3a). A second conservation gap area is located southward, between -40° to − 43°S latitude, including a suitability hotspot in central Chiloé island, where S. conspicua west and F. cupressoides showed concordant high suitability along the coastal range (Fig. 3b). When overlapping protected areas layer, we detected few, very small natural sanctuaries scattered within this high suitability hotspot (Fig. 3a, b). The southernmost region of the PTF, between − 43° to − 55°S latitude is well conserved in vast protected areas, where P. uviferum south, L. fonkii south, and P. nubigenus showed concordant high suitability values (Fig. 3c).

Discussion
Studying ENMs of nine species of gymnosperms and integrating niche overlap and estimation of levels of  rarity, allowed us to prioritize species and identify conservation gap areas in a global biodiversity hotspot. We generated reliable ENMs for nine Patagonian gymnosperms and found that some of their niches overlap, and only one species displayed a unique niche. Unexpectedly, we found that three species have divergent suitability of habitats along their range. We showed that the rarer a species is the smaller niche volume tend to have, and that there are unprotected suitable areas for key species of the PTF. We suggest that urgent conservation actions are needed, if we aim to preserve gymnosperms' niches in this global biodiversity hotspot.

ENM of PTF's gymnosperms
We obtained reduced ENMs that suggested both temperature and precipitation seasonality, as constrain factors depicting the environmental envelope of vegetation in PTF. Although, precipitation seasonality (BIO15) is present in nine out of twelve models deployed and in the region gymnosperms are usually considered to define the forest type (i.e. Araucaria forest) (Donoso 2006;Luebert and Pliscoff 2006), here we showed that each species can also harbor its particular set of environmental constrains. This species-specific set of constrains might be relevant to guide climate change mitigation projects through assisted migration programs for focal species in this hotspot (Casazza et al. 2021). This strategy can be especially relevant for endemic species that might not be able to follow rapid changes in climate, as they may have little dispersal capacity (Essl et al. 2011), and/or a reduced climatic niche, so negative effects of warming and climate change deepen their survival (Ozinga et al. 2009), particularly for cold tolerant species. Species that already show a narrow range may increase their vulnerability, if range shrinking due to climate change couples with other threats such as habitat fragmentation and degradation (Orsenigo et al. 2018), or land-use change (Miranda et al. 2016). Regardless of current and future climate change, the introduction of new actions in conservation sensu IUCN (2013) can compensate for species range loss due to threats other than climate change. Future studies, including climate projections would be of great value to elucidate potential responses of these gymnosperm species in this hotspot. Surprisingly, three focal species have divergent suitability of habitats across the landscape. Pilgerodendron uviferum and Lepidothamnus fonkii displayed two ecological groups separating their populations into north and south. This result is less surprising when we consider that these species have the largest latitudinal range across the PTF (encompassing more than 15° of latitude), occurring in heterogeneous environments. Both P. uviferum and L. fonkii, have small isolated populations, and these two species niches overlapped within sectors, although, we found within species divergent suitability of  (Table 2). Moreover, both south sectors were classified as discontinuous species with wide range and an intermediate rarity score of 4, while both north sectors were classified as improbable species and have narrow range with a high rarity score of 7. Suggesting that overlapped north sectors of P. uviferum and L. fonkii are more threatened than southern sectors. Instead, Saxegothaea conspicua showed longitudinal divergence, between populations occurring along the Andes Cordillera (in the east), and populations occurring along the Coastal range (in the west). The divergent suitability of habitats observed in three out of nine studied species is not related to previously observed geographic genetic structure measured with neutral molecular markers. These markers indicated local persistence in glacial refugia along their ranges at different latitudes for cold-tolerant species such as F. cupressoides (Premoli et al. 2000), P. uviferum (Premoli et al. 2002), P. nubigenus (Quiroga and Premoli 2010) and A, araucana (Bekessy et al. 2003;Martín et al. 2014). And for mesothermic species, such as A. chilensis, showed range retraction to milder climates during ice ages (Souto et al. 2015a). Extrapolated ENM to the past for A. chilensis are in concordance with this (Souto et al. 2015a). Modern molecular techniques would be more appropriate for further testing adaptive divergence particularly for species showing divergent suitability of habitats. Likewise, studies detecting characters under divergent natural selection could be useful to design conservation plans under predicted climate change conditions, to allow in situ adaptive divergence (Blanquart et al. 2013;Wang and Bradburd 2014). Studies with molecular markers, would be essential to corroborate environmentally triggered divergence hypotheses.
Within-species niche differences might suggest new hypotheses on ecological, genetic and evolutionary divergences and support conservation decisions.
High species diversity in these hotspots is often attributed to an accumulation of narrow-ranged species via high speciation rates and/or preservation of species over time, via low extinction rates (Jetz et al. 2004), usually hosting ancient and young species (Harrison and Noss 2017). In this way, niche divergences may favor genetic isolations that could lead to potential speciation. Therefore, the comparison between different geographic areas within the same species may be of interest to detect areas with divergent evolutionary potential. As ENMs are independent of genetic or phenotypic based adaptive responses, research of potential speciation mechanisms and the development of regional-scale conservation strategies might be of interest (Morente-López et al. 2020). From the evolutionary point of view, the observed disjunction between populations could indicate gene flow barriers, affecting sectors connectivity, that in turn could encourage reproductive isolation through genetic divergence (Blanquart et al. 2013).

Niche Overlap
We tested pairwise niche overlap to detect shared environmental constrains among gymnosperms, and compared between different geographic areas within the same species to detect evolutionary potential. Although, the number of statistically significant equivalent niche pairs represents only 17% of the total pairs compared, we found that eight out of nine of the studied species showed equivalent niches with at least one of the other gymnosperms of the PTF, and only A. araucana held a unique niche and occupies a set of environmental conditions not exploited by any of the other studied gymnosperms. It has been suggested that phylogenetically distant species will also be environmentally distant (Martínez-Méndez et al. 2016). Surprisingly, F. cupressoides and P. andina are parapatric and showed niche equivalency (Table 2). Conversely, we found that P. uviferum has multiple niche equivalences with phylogenetically distant species. We had considered here ancient lineages that had successfully surfed at least Pleistocene/Holocene climatic oscillations, so convergent and divergent niches can be explained by phylogenetic relationships, and other factors such as climatic conditions at the time of lineage diversification (Herrera 2020). Particularly in gymnosperms, climate adaptation has been shown to be genetically constrained, and plastic responses to Fig. 3 Areas of moderate to high overlapping suitability for Patagonian Temperate Forest gymnosperm species (yellowred) and preserved areas (in blue). a: Detailed overlapping suitability for Podocarpus salignus, Araucaria araucana and Prumnopitys andina in the Costal Range from -35° to -40° latitude; b: Details overlapping suitability for Austrocedrus chilensis, Fitzroya cupressoides and Saxegothaea conspicua west and east from − 39° to − 43° latitude; C: Moderate to high overlapping suitability for nine gymnosperm species and preserved areas of PTF ◂ climate appear to be relatively conserved and highly polygenic (Yeaman et al. 2016).

Rabinowitz's types of rarity
It has been suggested that magnitude of local abundance in combination with the extent of spatial distribution of a species may increase the risk of global extinction, with rare species being the most prone to it (Gaston 1994;Gaston et al. 2000). Thus, rare, endemic and improbable species, such as P. andina, A. araucana, F. cupressoides, L. fonkii north and P. uviferum north, behold higher risk, compared with widespread species, with large populations. We classified our focal species according to Rabinowitz's types of rarity and tested niche volume hypothesis that widespread species with low habitat specificity, and continuous populations have larger niche volume compared to narrow species, with high habitat specificity and scattered populations. We showed that there is a significant relationship between geographic range, population size and niche volume. Niche volume tend to be smaller when the level of rarity is higher, as previously proposed (Slatyer et al. 2013). Traditionally, much effort has been devoted to conserve rare species (Broennimann et al. 2005) being the Rabinowitz's approach broadly accepted for assessing forms of rarity, but sometimes the data needed to perform this classification is not available. We calculated niche volume for nine emblematic species of the PTF, from open sources of occurrence data. Our results and the recent broad and public availability of georeferenced data has made niche volume relatively easy to estimate for any species, suggesting it as a handy tool as a proxy of rarity for conservation purposes.

Conservation gaps in a global biodiversity hotspot
We combined ENM, and forms of rarity, to detect conservation gaps in a global biodiversity hotspot, deploying patterns of habitat suitability, niche overlap and calculating prioritization scores for nine gymnosperm species, key elements of the PTF biodiversity hotspot, complementing traditional approaches of biodiversity conservation. We observed that there are conservation gaps in this biodiversity hotspot by overlaying areas of high suitability for gymnosperm species of the PTF and protected areas in Chile and Argentina. Particularly, we identify two areas with moderate to high suitability or hotspots shared by two or more gymnosperm species, one on the northerncentral Coastal Range, and the other on the Andes Cordillera. These areas are in agreement with previous reports based on genetic diversity patterns (Souto et al. 2015b;Premoli et al. 2018). Thus, the areas that are currently suitable for the distribution of gymnosperm species are also areas of high genetic diversity and biodiversity, for these and other plant species (Durán et al. 2013;Moreira-Muñoz and Troncoso 2014;Scherson et al. 2014;Martínez-Tilleria et al. 2017;Fuentes-Castillo et al. 2020). These represent relatively stable areas since Holocene and should be considered priority if we aim to preserve gymnosperms' niches in the PTF. These species are also threatened by land use change, so we plea for concrete actions in the design and management of conservation areas for their environments. Currently, both land use change studies (Miranda et al. 2016) and protected area surveys (Moreira-Muñoz and Troncoso 2014) indicate a strong deterioration and absence of conservation areas, particularly in the Coastal Range as in Nahuelbuta (Otavo and Echeverria 2017;Premoli et al. 2018).
Protected areas have a core role in sustaining biodiversity, and native forests remnants are the key to preserve the entire ecosystem (Moreira-Muñoz and Troncoso 2014). Recently, in Patagonia, especially in Chile, a new strategy of public/private actions has been developed to create small biodiversity conservation areas, although with partial representation (Schutz 2018). Since such small areas, make insufficient representation of local biodiversity, and are biasedly distributed, and particularly missing in the Chilean central valley (Durán et al. 2013). Conservation planning requires setting priorities at the same spatial scale at which decision-making processes are undertaken considering all levels of biodiversity (Souto et al. 2015b). Many studies highlight the importance and need to carry out concrete conservation actions in central Chile (Scherson et al. 2014;Martínez-Tilleria et al. 2017;Fuentes-Castillo et al. 2020) given its high number of endemism, a strong impact of land use change into exotic species plantations, a high frequency of fires, added to the forecasted impacts of climate change for the next 100 years.
The current availability of spatially explicit occurrence data from multiple species can help to identify biodiversity hotspots and aid the prioritization of conservation areas to establish science-based protected areas that might preserve the evolutionary potential of key habitats and species. In our study, we found large areas with high overlapping suitability, and our results allow the identification of species that require concrete conservation actions. Particularly, the area between -35º and -38.5º south latitude, both in Andes and Costal ranges where Araucaria, P. salignus, and P. andina which are endemic, improbable and rare respectively have high overlapping suitability. The complete range of these species is within this area, with all their populations being highly threatened by anthropogenic disturbances. In addition, these species have high rarity scores, which increases their potential vulnerability, along with being paleoendemisms, making this area of high conservation value. Moreover, at this latitude P. salignus has low conservation scores, while P. andina has moderate scores, suggesting urgent conservation actions are needed. In Chilean central valley about -39.5° south latitude, P. salignus and L. fonkii north shared high overlapping suitability, both species have high rarity values and low conservation scores, again calling for conservation attention. Moving southward between -40 and -43° south latitude, along the coastal range, and including Chiloe Island, we found high overlapping suitability for P. uviferum north along with F. cupressoides and S. conspicua. These species are improbable, endemic or unlikely, respectively and all showed high conservation priority scores. In this region few, little and scattered conservations areas are settled, yet again claiming for conservation actions.
Our study supports the idea that conservation must contemplate large geographic areas and/or long-time scales so that evolutionary effects and events can progress (Lacher 2018). In such a process, the natural environment and ecosystem are protected and maintained so that all constituent species (known or unknown) are conserved and benefit (Jaisankar et al. 2018). Biodiversity protection and management through in situ conservation involves the creation of certain specific areas, be they national parks, sanctuaries, and biosphere reserves. Current conservation targets and goals have evolved over time from a more species-centric approach to one that considers a species' role in ecosystem functioning. We suggest that conservation efforts for the PTF, should consider such approaches, as they may be much more beneficial than specific actions for one or a few species in a small area. Although, Patagonian Temperate Forest is home of umbrella species, is a biodiversity and genetic hotspot, and represents a freshwater reservoir for the region, is also under many threats like highly invasive species' plantations, dump constructions and other anthropogenic disturbances. Thus, we appeal to make it a priority in the need for active conservation policies, supported by governmental agencies, the scientific sector and environmental groups. Our approach of overlapping habitat suitability of multiple species, and types of rarity highlights the relevance of integrative studies of species complexes aiming to aid conservation of key ecosystems, as global biodiversity hotspots.
Acknowledgements Laura Salazar and Paula Mathiasen provided technical assistance on early stages of the MS. Cristian Echeverria shared his personal records for P. salignus. Emiliano Quiroga provided R scripts to test niche overlap. Marlon Cobos for technical support on kuenm analysis. Chrystian Camilo Sosa Arango for technical support on Gap-Analysis. Fernando Sánchez and Sandro Boi (our husbands, whom shared our respective quarantine offices) for providing observations on early stages of the MS.
Author contributions PQ and CS: have the equal participation in the construction and design of the study, in writing, in development of the analyses performed, and in the interpretation results.

Data Availability
All data is available on-line repositories as are expressed in main manuscript.
Code availability All software used in this study are referenced.

Conflicts of interest Not applicable.
Ethics approval Not applicable.

Consent to participate Not applicable.
Consent for publication Not applicable.