Modelling the current and future biodiversity distribution in the Chilean Mediterranean hotspot. The role of protected areas network in a warmer future

Mediterranean Chile is part of the five recognized mediterranean‐type climates in the world and harbours a very rich floral diversity. Climate change has been reported as a significant threat to its biodiversity. We used the flora of Mediterranean Chile to analyse how biodiversity patterns, as measured by phylogenetic diversity, genus and species richness will respond to climate change scenarios and identify the areas that will harbour the greatest evolutionary potential and biodiversity richness. We also evaluated how these spatial patterns are depicted within the current network of protected areas.


| INTRODUC TI ON
Biodiversity loss has been reported as one of the most serious menaces to ecosystems, threatening human welfare and ecosystem resilience worldwide (Millennium Ecosystem Assessment, 2005;Rockström et al., 2009;Steffen et al., 2015). The Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES) includes the importance of biodiversity in a global context, defining Nature's Contributions to People (NCP) as all contributions, benefits or even negative contributions that people obtain from nature. These include ecosystem services, but also, for example the sense of identity and spiritual value of nature, and the evolutionary processes of biodiversity, that allow continued functioning of ecosystems (Pascual et al., 2017).
Despite the importance of biodiversity and its imminent threats, the evaluation of taxa and ecosystems for conservation purposes has been extremely complex, especially given limited resources (Wilson, McBride, Bode, & Possingham, 2006). One of the main challenges of any conservation strategy is the uncertainty caused by climate change phenomena, as any strategy needs to be able to predict the conditions that ecosystems will have to face in the future (Pressey, Cabeza, Watts, Cowling, & Wilson, 2007). To work with this uncertainty, some potential scenarios have been modelled to describe the changes that are expected to occur in the near future (Midgley, Hannah, Millar, Rutherford, & Powrie, 2002), and thus, biodiversity evaluations should promote the conservation of attributes that could sustain it given these future changes. An important attribute of biodiversity relates to the accumulation of diversity over evolutionary time (Jarzyna & Jetz, 2016), as the processes and interactions that model the distribution of biodiversity cannot be decoupled from their evolutionary history (Faith, 1992;Faith & Baker, 2006;Fisher & Owens, 2004;Forest et al., 2007;Pio et al., 2011;Sechrest et al., 2002;Vane-Wright, Humphries, & Williams, 1991). Thus, measures based on evolutionary history can be used as surrogates for feature diversity (sensu Faith, 1992), because taxa with shared ancestry usually share ecological and life-history traits (see, for example Cadotte, Cavender-Bares, Tilman, & Oakley, 2009;Wiens et al., 2010).
Indices based on evolutionary relationships have been instrumental in evaluating biodiversity across taxa and ecosystems (Winter, Devictor, & Schweiger, 2013). Phylogenetic diversity (PD), the most commonly used phylogeny-based index, is measured in a phylogeny as the sum of the branch lengths that connect a given set of taxa and provide an idea of how much evolutionary feature diversity would be lost if those taxa were not preserved (Faith, 1992).
When considering PD, one can identify areas that harbour a larger than expected evolutionary potential, and as proposed by Sgrò, Lowe, and Hoffmann (2011), this would allow identifying refugia that could maintain evolutionary features in the face of climate change. Previous studies have predicted changes in PD patterns, combining PD estimation with models of species distributions in several climate change scenarios Pio et al., 2011Pio et al., , 2014Thuiller et al., 2011). For example, Thuiller et al.
(2011) studied how general patterns of PD for three groups (plants, birds and mammals) were expected to change in Europe during this century under different climate change scenarios. The results were similar for all groups, providing possible scenarios that should be considered when planning conservation actions. Similarly, using a very diverse clade of eucalypts, González-Orozco et al. (2016) showed significant changes in future distribution of most species in the clade and a decrease in areas of palaeo-endemism.
Mediterranean climates are characterized by both high species richness and high levels of endemism (Ackerly, Stock, & Slingsby, 2014;Cowling, Rundel, Lamont, Kalin Arroyo, & Arianoutsou, 1996;Rundel et al., 2016) and, in the case of Chile, by having a larger increase in PD with species richness than for other mediterranean ecosystems (Morlon et al., 2011).
The Chilean Mediterranean region (MedCh) is considered a global hotspot of conservation priority due to its endemism and a high degree of habitat loss, mainly due to the conversion and replacement of the native vegetation (Armesto et al., 2010;Newton et al., 2012;Underwood, Klausmeyer, Morrison, Bode, & Shaw, 2009).
Because climate change has been reported as one of the significant threats to biodiversity , this study was set out to investigate how climate change scenarios may affect biodiversity patterns in a mediterranean biodiversity hotspot, using the Mediterranean area of Central Chile as a case study. For this area of the world, climate change predictions estimate a general decrease in precipitation and up to 4°C increase in temperature (Garreaud et al., 2017). The current protected area system for Chile protects 21% of the territory. However, this percentage is unbalanced throughout the country. Despite being inserted within one of the world's biodiversity hotspots, only 2% of the Mediterranean area of Chile is currently protected (see Marquet et al., 2004) and with national parks and reserves of reduced area and connectivity as compared to the larger areas concentrated in the Chilean Patagonia (Pliscoff & Fuentes-Castillo, 2011).
For this study, our main goal was to identify the areas, within the Mediterranean region, that will harbour the greatest evolutionary potential and richness of Chilean flora following climate change scenarios. We also aimed at evaluating whether the expected variation in future biodiversity patterns will be represented within the current PA system. We thus compared spatial patterns of richness at the species and genus levels, and PD at the genus level under current and future climate scenarios. We also evaluated gain, loss and turnover of species and genera richness towards future climate scenarios, and for Mediterranean endemic species, we related gain and loss to life-form. Finally, we estimated the capacity of the PA system to harbour future richness and PD.

| Species and genera database
We compiled a database for the Chilean Mediterranean flora based on herbaria locality information (Chilean National Museum of Natural History-SGO, Department of Botany of the University of Concepción-CONC and EIF-Herbarium at University of Chile).
Additionally, the Mediterranean portion of the Chilean flora database collected by Scherson et al. (2017) was considered. Occurrences were corrected eliminating outlier's records following known ranges from the literature when possible. Uncertain records were also excluded. Our database was composed of 94.048 unique records in Chile, including 2.661 species (72% of total Chilean flora) and 660 genera (82% of total Chilean genera). This database was then filtered to get the species list in the Mediterranean macrobioclimatic zone in Chile (Luebert & Pliscoff, 2017). Species life-form and origin were sourced from Zuloaga, Morrone, Belgrano, Marticorena, and Marchesi (2008).
We applied a distance filter in our database by removing species/ genera occurrences that were closer than 4 km from each other. The disaggregation function included in the Ecospat r package (Di Cola et al., 2017) was used to remove these occurrences. After the distance filter was applied, we obtained a final database with 1.727 species and 571 genera belonging to the Mediterranean region, considering a minimum occurrence number of 10 records. From the species selected, 79% had more than 15 unique records (see Appendix S2

| Current and future climate dataset
Current bioclimatic surfaces for Chile were sourced from Pliscoff, Luebert, Hilger, and Guisan (2014) considering a spatial resolution of 1 × 1 km and representing a period of 50 years . Global future climate surfaces are usually obtained from the WorldClim database. However, given the availability of the climate baseline developed by Pliscoff et al. (2014), with more spatial accuracy for this area, we incorporated this baseline to infer future climate predictions, and RCP8.5 (Riahi et al., 2011). We selected the CSIRO-MK3 Global Circulation Model (GCM) because its forecasts in annual mean temperature and annual precipitation for the Chilean Mediterranean region across all RCPs (2.6, 4.5, 6 and 8.5) are the most similar to the average predictions from the CMIP5 ensemble of models (Taylor, Stouffer, & Meehl, 2011). Differences between GCMs to identify the most suitable model for the region were assessed by using the shiny package 'GCM compareR' (see details in Appendix S1- Table   S1.1 and Figure S1.1; Fajardo, Corcoran, Roehrdanz, Hannah, & Marquet, 2018). The study was defined following the limits of the Mediterranean macrobioclimatic zone in Chile (Luebert & Pliscoff, 2017), according to two bioclimatic variables-annual mean temperature and annual precipitation. This definition considers this area to be located between 23-39° Lat. S and 74-70°Long.W.

| Species distribution models (SDMs)
The Maxent SDM algorithm (Phillips, Anderson, & Schapire, 2006) based on presence-only data was used for estimating Species Distribution Models at the species and genus level.
We selected the less correlated predictor variables to fit the models: annual mean temperature (bio1), temperature seasonality (bio4), maximum temperature of warmest month (bio5), annual precipitation (bio12) and precipitation seasonality (bio15). Variable selection was done by using the Pearson correlation between pairs of variables. Those variables that presented correlation values over 0.7 were discarded from the analysis. Pairwise correlation was done using the ENMTools package in r (Warren, Glor, & Turelli, 2010) (see heatmap correlation results in Appendix S2 as Figure   S2.1).
We used climate values from 10,000 randomly sampled points records as background data to account for sampling intensity, as well as to reduce spatial sampling bias (see Thornhill et al., 2017). We partitioned the occurrence locations data randomly into 80%-20%, for training and testing, respectively.
We used the 'distance hybrid' method and R script described in Thornhill et al. (2017) to constraint SDMs projections. This approach applies a Gaussian distance decay function to the species occurrence distribution to fit them in climatically suitable areas. In this study, a sigma of 500 km was applied for each species in current and future SDMs.
Binary range SDMs were estimated applying the specificity-sensibility threshold, converting continuous predictions to binary. The modelling extent was restricted to the complete area of Chile, which was used for all species and genera models. After that, results were restricted to the Mediterranean macrobioclimate zone.

| Future changes in biodiversity patterns
Biodiverse v.2.0 (Laffan, Lubarsky, & Rosauer, 2010) was used to estimate genus richness (GR) and phylogenetic diversity (PD) at the genus level. The phylogenetic tree of Chilean flora was sourced from Scherson et al. (2017), which was based on 756 genera, representing 87% of accepted Chilean native vascular plant genera, and pruned to include only the genera in the database. We considered a 10 × 10 km grid cell as it represents a fine scale to recognize spatial patterns. To be able to compare PD patterns with richness patterns obtained from the SDMs, we resampled SDM's spatial results to a 10 km × 10 km grid cell. Species richness (SR) was estimated by the sum of all binary range SDMs.
Differences between PD and GR were calculated according to Pio et al. (2011Pio et al. ( , 2014, using discrepancy values. These values were obtained per grid cell, normalizing PD and GR separately into z-scores, by subtracting each observed value from the mean, and dividing them by the standard deviation for all grids. Normalized PD was then directly subtracted from normalized GR for each grid cell (or Normalized Discrepancy according to Pio et al., 2014). Positive values indicate cells where the species present account for a higher discrepancy between PD and richness, meaning a larger amount of PD (see Pio et al., 2014), and the inverse applies for cells with negative values.
We calculated gain, loss and turnover of species and genera based on binary ranges to evaluate future changes on taxa distribution. Taxon gain, loss and turnover (Berteaux et al., 2018) were estimated as: Also, we conducted three linear regressions to describe the relation of genus measurements (GM i,j ) being explained by species measurements (SM i,j ) in both scenarios (i), and three measurements (gain, loss and turnover, each represented by a j). For all of them, we tested three possible relationships: Linear regressions and R 2 were estimated using r software (R Development Core Team, 2018). For each scenario and measurement, we fitted all three models and selected the relationship through AICc.

The percentage of gain and loss was calculated individually for
Mediterranean endemic species and related to life-form (herb, liana, shrub, ferns, succulent, tree and vine).

| Protected area network analysis
All metrics of future patterns of change were evaluated within protected areas (PAs) and outside of them, by considering their average variation towards future scenarios. We analysed the official Chilean state protected areas categories: Natural Parks, Natural Reserves, Natural Monuments and Nature Sanctuaries (Government of Chile, 1977, 1984, 1994, compiled from the national database of protected areas available at http://areas prote gidas.mma.gob.cl/. Furthermore, we evaluated the linear relation between protected area size and metrics of future patterns of change.

| Future patterns of change
Differences between current and future values for the biodiversity patterns measured (PD, GR and SR) are shown for two climate Taxon gain = 100 × taxon future gain∕taxon current richness Taxon loss = 100 × taxon future loss∕taxon current richness Taxon turnover = (taxon future loss + taxon future gain) ∕ (taxon current richness + taxon future gain) However, a higher general decrease in SR than in GR is projected in the Mediterranean zone of Chile for both future scenarios. Figure 2 shows the pattern of normalized discrepancy between PD and GR. For current climate, both the northern as well as the southern areas of MedCh have higher discrepancy implying that the assemblages found there exhibit a disproportionately larger fraction of PD than accounted by richness (i.e. areas of high evolutionary diversity and uniqueness). The opposite is observed in the central area of MedCh, where normalized PD is lower than normalized GR in the coastal area and central depression but remains higher than GR in the Andes mountain range. For both future scenarios, the difference between north-south and the central area intensifies, and the Andes mountain shows in these cases lower PD than expected by GR. This pattern appears more intense for the RCP8.5 scenario.
Even though gain, loss and turnover for species and genera under climate change exhibited different magnitudes for the RCP8.5 and RCP2.6 scenarios, the general spatial patters observed were similar.
Results are shown in Figure 3 for the RCP8.5 scenario and Figure   S2.5 in Appendix S2 for the RCP2.6 scenario.
Gain in both species and genera was higher towards the centraleast part of the region, in the Andes mountain range; meanwhile, their loss concentrated in the southern area of MedCh. Species turnover and genera turnover were higher in the Andes mountain ranges in the northern area of MedCh (see Figure 3).

| Protected area network analysis
Protected areas (PAs) in the Mediterranean region of Chile are present mainly in the Andes mountain ranges (Figure 5a) and coastal F I G U R E 2 Normalized discrepancy between phylogenetic diversity (PD) and genus richness (GR) for current (a), RCP2.6 (b) and RCP8.5 (c). Positive values (red scale) indicate cells where normalized PD is larger than normalized GR, implying that the species present in a cell account for a disproportionately large amount of PD. Negative values indicate cells where normalized PD is lower than normalized GR (blue scale) areas, with very few of them located in the central depression.
The total area occupied by PAs does not exceed 2% of the entire region. Figure 5 shows the spatial distribution of species and genus gain, loss, turnover and future PD at PAs in the most pessimistic scenario, where it can be observed that PAs located at the Andean mountain ranges will harbour species and genus gain. Taxon gain and turnover will increase at larger PAs; instead, species loss, genus loss and PD will decrease at larger PAs. Conservative scenario RCP2.6 maintained these tendencies (see Figure S2.8 in Appendix S2).
The average differences between taxon gain, loss and turnover inside and outside PAs are shown in Figure 6. It can be observed that average species gain and loss in both scenarios were higher than the ones observed for genera gain and loss, both within as well as outside the PAs. However, this difference was only significant between species and genus loss outside the PAs (Kruskal-Wallis test p < .001 for both RCP8.5 and RCP2.6 scenarios), and species and genus gain outside the PAs (Kruskal-Wallis test p < .001 RCP2.6 scenario). In general, gain was higher inside PAs, and loss was higher outside of them. Taxon turnover was similar for both species and genera, both within and outside the PAs.

Taxon richness (GR and SR) and PD inside and outside the
PAs-under current and future scenarios-are shown in Figure S2.9 in Appendix S2. Results show a slightly higher GR and SR average inside PAs for both future scenarios. Outside PAs, taxon richness and PD did not exhibit significant differences between current and future scenarios.
Finally, we estimated these trends for each PA within MedCh (see Appendix S3 as Tables S3.2 and S3.3).
F I G U R E 3 Gain, loss and turnover for species and genus richness.

| Current and future patterns of biodiversity
Our study allowed to compare the current and future patterns of PD, GR and SR for mediterranean ecosystems in Chile and the identification of important areas in terms of evolutionary potential and vulnerability to climate change. We found that current PD and taxon richness patterns are similar as has been suggested in previous studies (Forest et al., 2007;Scherson et al., 2017). Nevertheless, under different climate change scenarios the differences between current and future climatic conditions followed similar patterns only for PD and GR. This response can be explained by the tight correlation (0.99 correlation coefficient for current PD ~ GR) between PD and GR (Forest et al., 2007;Scherson et al., 2014). However, by computing the discrepancy of PD-GR we identified the areas where PD is higher or lower than predicted by GR. Areas in the north-south coastal range of MedCh show higher PD than predicted by richness in future scenarios. Moreover, differences in PD, GR and SR consistently identify the east of the region, upward in the Andes mountain ranges, as an area of large changes. This pattern of change, coincides with those reported for Mediterranean-type climate zones worldwide (Ackerly et al., 2014;, where the shifts in the range of climates has been reported to occur poleward and upward and are in agreement with the current trend in climate change within MedCh, where harsher conditions (annual and summer precipitation downward trends) have been identified from 1975 to 2015 (Deitch, Sapundjieff, & Feirer, 2017).
Our findings for future biodiversity patterns showed a considerable difference between climate change scenarios, being RCP8.5 where the current-future differences were larger and more evident.
For central Chile, it has been recently shown that the RCP2.6 and RCP8.5 scenarios can present large differences by the end of the century for projected warming and drying, and a larger range of precipitation uncertainty has been documented for the RCP2.6 scenario (Bozkurt, Rojas, Boisier, & Valdivieso, 2018). In general, we found a decrease in SR for the entire Mediterranean Region according to both scenarios, although this is stronger for RCP8.5. F I G U R E 5 Spatial distribution of species and genus gain, loss and turnover in protected areas of MedCh for the RCP8.5 scenario. Protected areas are shown in panel a. Taxon gain is shown in blue graduate circles (b & f), and taxon loss is shown in red graduated circles (c & g). Additionally, taxon turnover is shown in green graduated circles (d & h). Finally, genus phylogenetic diversity is shown in orange graduated circles (e). Linear regressions between protected area size and taxon gain/loss/turnover and PD are shown at the left of each map Ideally, analyses of future changes in PD spatial patterns should be done at the finest taxonomic and spatial scales possible. However, the data problem is a very complicated setback for most areas in the world (Cardoso, Erwin, Borges, & New, 2011;Scherson, Fuentes-Castillo, Urbina-Casanova, & Pliscoff, 2018).
SDMs have accounted for paucity of georeferenced data (Fois, Cuena-Lombraña, Fenu, & Bacchetta, 2018), and as this study shows the possibility of using them to extrapolate occurrences in areas where fewer locality information is found, taking into account limitations regarding number of independent occurrences and uneven distribution of species information in the datasets.
However, not enough DNA data are yet available-at least for plants-for doing phylogenies of floras of very large areas at the species level. As a consequence of that, most studies of spatial patterns of PD done so far have been either done in small areas at the species level (for example Zhang & Zhang, 2017), or in large areas at larger taxonomic scales, genera (Forest et al., 2007;Scherson et al., 2017;Thornhill et al., 2016) or OTUs .
For this study, we used a previously generated phylogeny at the genus level (Scherson et al., 2017). We do not yet have the phylogeny of this flora at the species level. Thus, the data presented in this study for changes in PD patterns can only be interpreted at the genus level. Even though we did observe positive correlations among changes in gain, loss and turnover between genera and species, these correlations are not high enough to assume that we can predict patterns of species PD in the present or future, given the results obtained with genera.

| Protected area network analysis
Land use change, climate change and the occurrence of extreme climatic events are expected to become more severe at the end of the century (Sala, 2000). MedCh is located within the Chilean hotspot (Arroyo et al., 2006;Mittermeier, Turner, Larsen, Brooks, & Gascon, ;Myers, Mittermeier, Mittermeier, DaFonesca, & Kent, 2000), and as such, having a resilient system of PAs should be a national priority. Mediterranean ecosystems in Chile have been shown to be affected by several conservation issues. Intense land F I G U R E 6 Protection level variation in species and genus gain, loss and turnover, for RCP8.5 and RCP2.6 scenarios for 2080, inside and outside the protected areas of MedCh. For the Box plots, red points indicate average values and p-value refers to Kruskal-Wallis analysis use change is concentrated in the Central Depression of MedCh, where native vegetation replacement by exotic tree plantations and agriculture (Armesto et al., 2010;Miranda et al., 2017) at accelerated rates in both space-time have been documented (Echeverria et al., 2006;Schulz, Cayuela, Echeverria, Salas, & Rey Benayas, 2010). The increase in fire events (Urrutia-Jalabert, González, González-Reyes, Lara, & Garreaud, 2018) as a consequence of anthropic influence and synergistically affected by a megadrought (González, Gómez-González, Lara, Garreaud, & Díaz-Hormazábal, 2018) has been another cause of the loss of native vegetation (Bowman et al., 2018). Forest fires pose large pressures over PAs and the matrix among them. Along with invasive species and land restoration, controlling forest fires has been identified as a priority action to improve conservation . Despite these problems, the area has very few PAs; instead, smaller national parks and reserves are concentrated in the Andes mountain ranges mainly (Pliscoff & Fuentes-Castillo, 2011).
Our study analysed current and future biodiversity patterns of PAs in MedCh and quantified expected taxon (both genera and species) gains and losses under future climate change scenarios. Our main findings show that on average, taxon gain is to be expected within PAs, and higher average loss will occur outside of them. However, patterns of gain and loss among PAs show a complex geographical pattern. PAs located along the Andean mountain ranges showed higher taxon gain, whereas PAs located in the central depression will concentrate higher taxon loss. Regarding PAs size, taxon gain and turnover will increase in larger PAs. Two PAs showed the high- As SDMs assume niche conservatism, the fact that the PA system will gain species suggests that most PAs in the study area are going to maintain a mediterranean-type climate at least until 2080 in both RCPs scenarios. Klausmeyer and Shaw (2009), analysing future climate patterns for all Mediterranean areas in the world, reported a relatively stable projected status for MedCh. Our results mostly agree within the Andean mountain range, where most PAs are located. However, our results show a decrease in PD and taxon richness in the Central Depression of MedCh, thus suggesting future climatic instability for this area. The same authors suggested that PAs in mediterranean-type climates have the capacity to provide ecosystem stability to mitigate climate change effects (see ).

| Conservation implications
In this study, we identified the areas within the Mediterranean region that will harbour the greatest evolutionary potential-through PD-and richness of Chilean flora, following climate change scenarios. We also evaluated whether this expected variation will be protected by the current PA system. Future biodiversity patterns as well as taxon gain, loss and turnover have been reported here for the first time in the Chilean Mediterranean Region. In this sense, our findings enhance the importance of the current PAs to harbour this future variation, despite their reduced number and size along the region. However, one of the expected effects of climate change is the reshuffling of ecological species (Alexander, Diez, Hart, & Levine, 2016), which is likely to be stronger in zones where a higher species gain is predicted to occur. This can generate new conservation issues that will need to be addressed.

ACK N OWLED G EM ENTS
We thank CONICYT grants 3190433, PII20150091, 1171586 and 1181677, AFB17008 and a grant from the Global Environmental Facility (GEF-5810). We also thank the curators of herbaria in Chile (CONC, SGO and EIF) who provided locality information.