Climate change is creating a mismatch between protected areas and suitable habitats for frogs and birds in Puerto Rico

Climate change is altering the spatial distribution of many species around the world. In response, we need to identify and protect suitable areas for a large proportion of the fauna so that they persist through time. This exercise must also evaluate the ability of existing protected areas to provide safe havens for species in the context of climate change. Here, we combined passive acoustic monitoring, semi-automatic species identification models, and species distribution models of 21 bird and frog species based on past (1980–1989), present (2005–2014), and future (2040–2060) climate scenarios to determine how species distributions relate to the current distribution of protected areas in Puerto Rico. Species detection/non-detection data were acquired across ~ 700 sampling sites. We developed always-suitable maps that characterized suitable habitats in all three time periods for each species and overlaid these maps to identify regions with high species co-occurrence. These distributions were then compared with the distribution of existing protected areas. We show that Puerto Rico is projected to become dryer by 2040–2060, and precipitation in the warmest quarter was among the most important variables affecting bird and frog distributions. A large portion of always-suitable areas (ASA) is outside of protected areas (> 80%), and the percent of protected areas that overlaps with always-suitable areas is larger for bird (75%) than frog (39%) species. Our results indicate that present protected areas will not suffice to safeguard bird and frog species under climate change; however, the establishment of larger protected areas, buffer zones, and connectivity between protected areas may allow species to find suitable niches to withstand environmental changes.


Introduction
Changes in species distributions are one of the most evident fingerprints of climate change. Species range shifts towards higher elevations and latitudes have been documented for a variety of functional and taxonomic groups in different regions of the world (Parmesan and Yohe 2003;Chen et al. 2011Chen et al. , 2012. Changes in species distribution due to climate change have been linked to population declines (Rosenberg et al. 2019; Sánchez-Bayo and Wyckhuys 2019), local extirpation (Wiens 2016; Campos-Cerqueira and Aide 2017) and extinctions (Monzón et al. 2011). In addition, projections based on meta-analyses indicate that this trend will continue in the future, threatening the conservation of fauna worldwide (Nunez et al. 2019). Reflecting the global nature of climate change, population declines have been documented even in pristine and protected areas (Blake and Loiselle 2015). As climate changes, species that presently occur in protected areas may also shift their ranges outside the boundaries of the protected areas, increasing their vulnerability to other human impacts (Monzón et al. 2011). There is thus an urgent need to understand how climate change is altering species distributions relative to distributions of protected areas.
While protected areas are still our primary strategy to conserve species, communities, and ecosystems (Rodrigues et al. 2004;Saout et al. 2013), the recognition that our current network of protected areas may not efficiently protect some species under climate change scenarios dates back to the late 1980s (Hannah 2008). Despite an increase in studies testing and quantifying the role and effectiveness of protected areas in conserving species under future scenarios (e.g., Araújo et al. 2004Araújo et al. , 2011Welch 2005;Hannah et al. 2007;Mackey et al. 2008;Hole et al. 2011;Alagador et al. 2014;Gaüzère et al. 2016), much debate still exists regarding the ability of protected areas to mitigate the effect of climate change on species distributions. Some studies have shown that existing protected areas will be able to mitigate the effects of climate change on bird species distributions in France (Gaüzère et al. 2016), to offer native species refuge from invasive species under climate change worldwide (Gallardo et al. 2017), to facilitate adjustments by bird species to rapid climate conditions in Finland (Lehikoinen et al. 2019), to highly overlap (> 45%) with high diversity of endangered anuran species in Costa Rica (García-Rodríguez et al. 2012), and to slow bird species decline in Finland (Virkkala et al. 2019). Nevertheless, species occurrence inside protected areas are expected to change because protected areas worldwide will fail in retaining more than 90% of their available climate (based on mean annual temperature and annual precipitation) at current levels (Elsen et al. 2020). In Europe, an in-depth study encompassing a network of protected areas and 75% of terrestrial vertebrates (n = 585) found that current unprotected areas will be more effective in retaining climate suitability for species than existing protected areas (Araújo et al. 2011). An assessment of the distribution of 1200 plant species in Europe found that 6-11% of species would be absent in the existing reserves in the future (Araújo et al. 2004). In the tropics, the distribution and suitable habitat of bird and frog species in the Cerrado and Atlantic Forest (Lemes et al. 2014;Vale et al. 2018;Borges et al. 2019;Velazco et al. 2019) are projected to decline inside protected areas. Given that the population of many species are expected to decline inside protected areas in the future, studies assessing long-term suitable areas for species occurrence outside the boundaries of protected areas are urgently needed. By identifying these suitable areas that are stable through time for different species and taxonomic groups, we can adapt our conservation strategy by creating, expanding, moving, and connecting protected areas to accommodate population range shifts in the landscape.
To assess the role and the effectiveness of protected areas in the future, studies have applied species distributions models (SDM) to evaluate the effect of climate change on species distributions and conservation strategies. Use of SDM is the most common approach used to evaluate the effectiveness of protected areas in conserving species under climate change scenarios (Seo et al. 2009;Araújo et al. 2011;Hannah et al. 2020), and the implementations of these models, as well as their uncertainties, have been extensively reviewed (Araújo and Guisan 2006;Heikkinen et al. 2006;Elith and Leathwick 2009). One common approach to assess the impact of climate change on species distribution is identifying suitable areas in both present and future scenarios (Araújo et al. 2004;Wang et al. 2017;Báez et al. 2019;Beaumont et al. 2019). Nevertheless, the current distribution of a species is only a snapshot of a dynamic process that holds a historical legacy (Maguire et al. 2015). By identifying areas that remained suitable from past climate change, and that will stay suitable in the future, we can prioritize areas for conservation that are more likely to stay stable against climatic fluctuations. Unfortunately, few of these studies have included historical projections to assess habitats that maintain suitable conditions for the persistence of species in the long-term. Moreover, most studies rely on presence-only data (Hao et al. 2019) and often focus on a single species hindering our understanding of community response to climate change. Presence-only data do not provide information about the environmental variables associated with areas where the species are absent, cannot estimate prevalence (proportion of occupied sites) and suffer from sampling bias (Guillera-Arroita et al. 2015). Thus, the dynamic response of species to climate change still poses a challenge for the establishment of successful long-term protected areas networks.
This study focuses on areas that maintain suitable areas for several species of birds and frogs over past, current, and future climate regime periods. These always-suitable areas (ASA) may facilitate species' persistence and should be prioritized for conservation. Given that most protected areas were not established considering climate change or non-static species ranges, we asked the following question: Do the distributions of always-suitable areas for species of birds and frogs overlap with the distribution of existing protected areas? To answer this question, we determined the species distributions of eight frog and 13 bird species based on historical, current, and future climate models and detection/nondetection data from ~ 700 sites and > 900,000 1-min audio recordings in Puerto Rico.

Study area
Puerto Rico is an archipelago comprising the islands of Puerto Rico, Culebra, Vieques, and several small islands. Puerto Rico is the smallest of the Greater Antilles (~ 8900 km 2 ) and has a rugged topography with an elevational range from sea level to 1300 m. The spatial distribution of rainfall in Puerto Rico is variable, being highest in the Sierra de Luquillo rainforest in the eastern part of Puerto Rico (4305 mm/year), and lowest near Guánica in southwestern Puerto Rico (768 mm/year). Forests covered less than 10% of the main island in the late 1940s when the human population and agriculture (sugarcane, coffee, plantain, cotton, corn, and rice) were widespread in Puerto Rico (Grau et al. 2003). Since then, most agricultural areas were abandoned, the human population is declining (Castro-Prieto et al. 2017), and forests have recovered (Grau et al. 2003). Currently, forests cover more than 54% of the main island (Marcano-Vega 2017).

3
Most of the protected areas in Puerto Rico (~ 71%) are in the interior mountains coinciding with the location of the original remnants of the 10% of forests that have persisted after the arrival of the Spanish in 1493. In addition, most forest recovery occurred at higher elevations (Grau et al. 2003). Protected areas in Puerto Rico are mostly small (1-115 km 2 ), cover 16.1% of the total land surface (1436 km 2 ), and include public and privately-owned land (Castro-Prieto et al. 2019). We focused our comparison on 150 terrestrial protected areas that are co-managed by the Department of Natural Resources of Puerto Rico, US Fish and Wildlife Services, USDA Forest Service and the Conservation Trust of Puerto Rico.

Data collection
We used acoustic data stored in the RFCx-ARBIMON platform (https:// arbim on. rfcx. org/ home) to acquire detection/non-detection data on 13 bird and eight frog species. We selected audio recordings from 11 ecological projects conducted in and around several protected areas and including ~ 700 sampling sites from 2015 to 2019 (Table 1, Fig. 1). All 11 projects aimed to understand bird and frog distribution patterns in different regions of the island. MCC, TMA, JAC, and BAM were directly responsible for the deployment of audio recorders and the execution of 10 of the 11 projects. The selection of sampling sites within each project followed an approximate systematic sampling design to capture the full range of conditions found in each region. The only exception was the El Yunque National Forest project, where sites were deployed in three transects along an elevational gradient (see Campos and Aide 2017 for more details on sampling design). Approximately 50% of the sampling sites were inside protected areas and in high elevation areas (> 400 m).
All audio recordings were collected using portable acoustic recorders (L.G. Android and Audiomoth), which were programmed to record at a sampling rate of 48 kHz with 24 kHz bandwidth, using a medium gain (30.6 dB) and following a schedule of 1-min 1 3 audio recording every 10 min (i.e., in total 6 min of audio recordings per hour), 24 h per day for 1-2 weeks per sampling site. All recorders were placed at least 200 m from the nearest recorders on trees at the height of 1.5 m, and ~ 50 m from any road with traffic. Previous field tests conducted by our team indicate that the recorders can detect vocalizations of all focus species up to ~ 100 m. Because the historical and future climate data used in the analyses have a spatial resolution of 2 × 2 km, sampling sites were collapsed at the scale of the grid cell. In this way, we transferred values associated with the sampling points (detection/non-detection) to the spatially overlapping raster cells so that each point is assigned to a grid cell. The value of a grid cell was determined by the maximum values associated with the sampling points. In that way, a grid cell with points with detections and non-detections was assigned as occupied.

Species of greatest conservation need (SGCN)
We selected species (Table 2) that were listed as SGCN by the most recent Puerto Rico State Wildlife Action Plan (PRSWAP 2019). Most SGCN have low and/or declining populations and are in need of conservation action. Teams composed of species experts used a standardized decision process to evaluate rarity, threats, and trends for each species according to NatureServe Conservation Status Assessment. The main objective of the PRSWAP is to identify and address the greatest conservation needs of Puerto Rico`s fish and wildlife populations and to prioritize efforts for species with the greatest conservation needs. Due to modeling limitations, we removed species that occurred in less than six sampling sites.
Automated classification algorithms were used in the RFCx -ARBIMON II platform to determine the presence and absence of each species in the audio recordings. Classifications of all recordings were based on a template matching procedure. This procedure searches through the audio data (~ 900,000 1-min recordings) for acoustic signals and detects regions that have a high correlation with a template that has been selected by the user. All regions of interest (ROIs) with values above the selected correlation threshold (0.1) were presented as potential detections for posterior validations (LeBien et al. 2020). The selection of a low threshold resulted in a high number of false positives though the number of false negatives was negligible. MCC and two other experts in Puerto Rico fauna (Noelia Nieves and Nashally F. Mercado) manually reviewed the template matching results. In this step, we annotated the results as either positive or negative, indicating the corresponding

Climatic data and projections
Monthly precipitation estimates for the period 1980-2018 (39 years) were derived by blending data from radar imagery from the National Oceanic and Atmospheric Administration (NOAA; https:// water. weath er. gov/ precip/) and rain gauge measurements from reporting stations obtained from the Global Historical Climatology Network (GHCN; Menne et al. 2012). This method attempted to reconcile two shortcomings in these datasets: (1) the radar data have a short period of record (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) and can suffer from biases compared to physical measurements of rainfall and (2) the rain gauge data are sparsely located and have significant observational gaps in the period of record. Given these two constraints, we developed an approach that leveraged aspects of both datasets. The radar data were used to estimate a seamless precipitation climatology, providing remotely-sensed observations in areas where rain gauge measurements are absent. The rain gauge data were used to estimate monthly precipitation anomalies (relative to the short climatological period) over the 39-year period of interest. We used a 4-step procedure to assimilate data from these two independent sources, resulting in a single harmonized estimate of gridded (2 km × 2 km resolution) monthly precipitation. First, monthly climatologies were calculated for both the radar data and the rain gauge network for the overlapping 12-year period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) and were then used to calculate differences between the rain gauges and the co-located radar data grid cells (i.e., the bias). Second, a geostatistical model (ordinary kriging with an elevation co-variate predictor) of the monthly climatological bias was built to predict the level of bias in the radar data at each 2 km × 2 km grid cell. Note that in our bias calculations we are assuming that the rain-gauge measurements of precipitation represent 'truth'. We then subtracted the predicted monthly climatological bias from the gridded radar climatology to obtain the 'bias-corrected' precipitation climatology. Third, using the R package 'fields' (Furrer et al. 2017), we fit a Thin Plate Spline (TPS) model to interpolate the variancestandardized monthly rain gauge anomalies (see Wilks 2011), and included elevation as a covariate to explain spatial variation. Finally, the interpolated rain gauge anomalies were then added to the bias-corrected radar climatology to produce estimated monthly precipitation amounts for the 2 km × 2 km resolution grid. A similar but simpler approach was used to estimate monthly maximum and minimum temperature over the same time period . A 'smart-interpolation' method was used where an assumed lapse rate of − 6.8 °C/1000 m was first applied to the observed temperature dataset (~ 25 high-quality weather stations from the GHCN dataset) to estimate a monthly climatology of 'sea level' temperatures. We then used a TPS model to interpolate these climatological values to a 240 m resolution grid. We chose to first interpolate to a higher resolution (240 m) before aggregating to the common 2 km resolution grid to better resolve the fine-scale terrain effects on temperature. Prior work (Daly et al. 2003) has shown that with the available station network density on the island, accurate interpolation could be accomplished at this resolution. This is expected to reduce bias in the final estimated 2 km resolution used in the full historical and future climate analysis.
The TPS model included a coastal-distance index to account for a moderating influence of the ocean on minimum temperatures during cooler months (Daly et al. 2003). Next, an inverse-distance weighting (IDW) model was used to interpolate monthly temperature anomalies from the raw station data. These anomalies were then added to the TPS-modeled climatologies after the lapse rate adjustment was re-applied to obtain estimated temperatures at the DEM-derived elevations. Finally, the monthly gridded temperatures were aggregated to the common 2 km × 2 km resolution grid.
Projected precipitation and temperature were obtained from the outputs of two dynamically downscaled GCMs for the period 2040-2060 and compared to 1985-2005, using the Weather Research and Forecasting (or WRF) model (Bowden et al. accepted). These climate model simulations were run under the RCP8.5 scenario, a higher greenhouse gas emissions scenario at the upper range of what is considered the plausible bounds of future emissions trajectories. Dynamically downscaled climate models were used because the high spatial resolution (2 km × 2 km) makes it possible to resolve small-scale climatic processes and the influence of local topography and land surface energy exchanges on the island's climate. In addition, dynamical downscaling was used because the explicit resolution of atmospheric processes improves representation of physical feedback processes, such as cloud formation and moisture fluxes. These processes are difficult to account for in empirical models (Khalyani et al. 2016) that assume stationary relationships between models and observations through time. We note that with only a single climate scenario and two GCMs, there is likely to be a large portion of the plausible future climate space that will not be sampled by the model space. Nevertheless, for the time period under consideration, the climate model results are consistent with the majority of climate model projections for the region, which suggest widespread warming and drying in the eastern Caribbean (Chadwick 2017;Bhardwaj et al. 2018). From the model outputs, we obtained the absolute temperature change and the relative precipitation change and applied them to our interpolated dataset of historical monthly temperature and precipitation to obtain a bias-corrected estimate of the projected change in these variables. Finally, 19 bioclimatic variables were derived from the temperature and precipitation estimates from 1980-2018 and 2040-2060 using the R package "dismo" ). Nevertheless, we further reduced the number of bioclimatic variables after excluding highly collinear variables for the analyses. We have also subdivided the 1980-2018 time period into past (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989) and present (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) periods. The historical and current periods were selected based on periods that depicted an average representation of past and current climate without large natural disturbances (e.g., droughts and hurricanes). Because elevation is a good proxy for habitat types, intensity of land use, and animal and plant communities we have also included elevation along with the bioclimatic variables as predictors in the species distribution models.

Variable selection
To avoid multicollinearity among the 19 candidate bioclimatic and physical covariates in the species distribution models, we used the VIF function (Variance Inflation Factor) from the "usdm" package (Naimi 2015) in R version 3.6 (R Core Team 2017). This function progressively excludes collinear covariates through a stepwise procedure. VIFs were calculated using two methods: VIFcor (threshold = 0.7) and VIFstep (Threshold = 10). VIFcor method first find a pair of variables which has the maximum linear correlation (greater than threshold) and exclude the variable with greater VIF. This procedure is repeated until no variable with a high correlation coefficient (greater than threshold) with other variables remains. VIFstep calculates VIF for all variables, then excludes the variable with the highest VIF (greater than threshold). The procedure is repeated until no variables with VIF greater than threshold remains. Only variables that were selected by both the VIFcor and VIFsteps approaches were used in the analyses. The same variables, at the same spatial resolution, were used to project species distributions into past and future climate scenarios (Table 3).

Species distribution models
The detection/non-detection data of the 21 species were used to create species distribution models for the current time period using nine methods commonly used in species  (Hao et al. 2019): generalized linear models (GLM), random forest (RF), support vector machine (SVM), boosted regression trees (BRT), classification and regression trees (CARS), maximum likelihood (MAXLIKE), multivariate adaptative regression spline (MARS), multi-layer perceptron (MLP) and radial basis function network (RBF). Models were created with default tunings (i.e., default settings for each fitted model). Models were built using ten runs using the subsampling data splitting method (n = 10), so for each run, 80% of the data were used for training, and 20% percent for testing. Model performance was evaluated by two metrics (Allouche et al. 2006): (a) the area under the curve (AUC) of the receiver operating characteristic (ROC), a threshold-independent statistic that ranges from 0 (model equivalent to a random prediction) to 1 (perfect model discrimination between presence and absence records); and (b) the true skill statistic (TSS), a threshold-dependent statistic that ranges from − 1 (model equivalent to a random prediction) to 1 (perfect model fit). Species distribution models with AUC < 0.6 or TSS < 0.3 were excluded from the ensemble procedure because their predictive powers were similar to random predictions. Predictions were made through ensemble models using the weighted averaging method based on TSS. Despite some recent critics about the use of ensemble models (Hao et al. 2020), this approach is often recommended to address the uncertainty associated with using multiple modeling methods (Araújo and New 2007;Thuiller et al. 2019). The use of multiple modeling methods to assess species distributions is reasonable, considering that no single model will capture the complexity of environmental factors influencing species occurrence in time and space. By taking advantage of different modeling methods, we assume that each modeling method can contribute to assessing the species' distribution. The models created for the current time period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) were then used to project species distributions for the historical and future periods using ensemble models and the climatic variables from 1980-1989 and 2040-2060, respectively. The historical and current time periods were selected based on time periods that depicted an average representation of past and current climate without large natural disturbances (e.g., droughts and hurricanes). The relative importance of each variable was evaluated through a randomization procedure that measures the correlation between the covariates' predicted values and predictions in which the covariate under investigation is randomly permuted. If a variable's contribution to the model is high, then it is expected that the prediction is more affected by a permutation, and therefore the correlation is lower.
We generated binary maps for each species and each time period based on the mean threshold value among the different algorithm predictions of the sensitivity-specificity equality approach (Liu et al. 2005). We also overlaid the binary maps for all three climate scenarios (historical, current, and future) to obtain distribution maps for always-suitable areas (i.e., areas with suitable habitats across the three climate scenarios). Contrary to the more common use of the current and future distribution to assess climate refugia, our ASA approach must agree in all three time periods (past, present, and future), and therefore we provide a more conservative approach in identifying potential areas for protection. ASA can be better understood as areas with high habitat stability that can allow the longterm persistence of the species, and therefore, ASA can be interpreted as a habitat refugia (Ashcroft 2010). By adding the historical component, we are taking into account possible historical range contractions due to biotic (e.g., diseases) and anthropogenic (e.g., urbanization) constrains that may underestimate species niche and influence future projections (Nogués-Bravo 2009;Maguire et al. 2015;Faurby and Araújo 2018;Santini et al. 2020). In that way, ASA focus on areas that consistently provide safe havens for species under changing climate and may provide the best chances for the survival of species under future climate change. The identification of these suitable habitats can, therefore, help in prioritizing conservation areas, thus providing relevant information to policymakers. In addition, we created maps displaying the overlap of always-suitable areas for (1) birds, (2) frogs, and (3) birds and frogs combined and we calculated the percent and extent (km 2 ) of total and accumulative ASA that overlap with protected areas and the percent and extent of total and accumulative ASA that do not overlap with protected areas.
All spatial analyses were carried out using the packages "raster," "usdm," and "sdm" implemented in R 3.6 program language (R Core Team 2017) and QGIS software (Version 3.6).
One limitation of our study is that our modeling approach does not account for species detectability. Therefore, our models may be unrealistic and may underestimate species distributions (Lahoz-Monfort et al. 2014;Guillera-Arroita et al. 2015). Nevertheless, we have compared outputs from both occupancy models that take into account species detectability and ensemble models using the SDM package, and we have found similar patterns in species distributions (Campos-Cerqueira et al. unpubl). Therefore, we have chosen to use the SDM package because it provides a user-friendly platform that allows quick analyses of multiple species simultaneously using several modeling approaches.

Results
The distributions of 21 bird and frog species were based on detection/non-detection data from > 900,000 1-min recordings from ~ 700 sites across Puerto Rico. Model performance assessed by average AUC varied from 0.63 (Icterus portoricensis) to 0.94 (Eleutherodactylus unicolor) (Fig. 2, Table S1). Model performance accessed by averaged TSS varied from 0.31 (Melanerpes portoricensis) to 0.89 (E. unicolor).
Precipitation in the warmest quarter was the most important variable affecting frog distributions whereas temperature seasonality and precipitation in the warmest quarter were the most important variable affecting bird distribution (Table S2). Precipitation in the warmest quarter increased from the historical (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989) to current time period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) across the whole island (Fig. S1). The increase was less than 25 mm for most of the island, but there were areas in the northern karst region and the Central Cordillera where the increase was greater. In contrast, there was an overall drying trend from current conditions to the future . The models suggest that precipitation in the warmest quarter will decline for much of the island by 150 mm or more, and in parts of western Puerto Rico precipitation could decline by at least 325 mm. An exception to this trend is a portion of the Luquillo Mountains, which is projected to have a slight increase in precipitation (Fig. S1). The distributions of birds varied greatly (Fig. S2). Vireo altiloquus was the most widespread species in all periods, occurring in ~ 80% of Puerto Rico, followed by Loxigilla portoricensis and Melanerpes portoricensis. All other species had limited distributions, occurring in less than 30% of Puerto Rico in any period. Three species (Geotrygon montana, Nesospingus speculiferus, and Setophaga angelae) had very limited distributions occurring in less than 10% of Puerto Rico. ASA were very small (less than 10% of Puerto Rico) for seven bird species (Buteo platypterus, Antrostomus noctitherus, Geotrygon montana, Icterus portoricensis, Nesospingus speculiferus, Setophaga angelae, and Vireo latimeri). The stable areas of three of these seven species (Antrostomus noctitherus, Vireo latimeri, Icterus portoricensis) occurred mostly in the west portion of Puerto Rico. In contrast, the ASA for the other four species were mainly concentrated in montane rainforest habitats.
Most anuran species had limited distributions occurring in less than 10% of Puerto Rico (Fig. S2). In addition, ASA were extremally small (< 2% of Puerto Rico) for the majority of frog species. The grass coqui (Eleutherodactylus brittoni) was the most widespread species occurring in over 39% of Puerto Rico. Eleutherodactylus unicolor was detected only in the Luquillo Mountains in both historical and current time, and its distribution in all scenarios is restricted to this region. Eleutherodactylus gryllus, and E. portoricensis are both currently restricted to the Luquillo and Cayey Mountains. Although historically these species occurred in other regions (e.g., Central Cordillera), the Luquillo Mountains are the only ASA for these species. Eleutherodactylus locustus was detected in the Luquillo and Cayey Mountains in both historical and current time periods, and the always-suitable areas encompass both regions and a small portion of the Toro Negro State Forest. Eleutherodactylus hedricki was detected in the Luquillo, Cayey, and Toro Negro Mountains, both in historical and current time periods, but always-suitable areas encompass only the Luquillo Mountains. Although historically E. richmondi occurred in several locations in Puerto Rico, currently, this species is mainly restricted to the Luquillo and Cayey Mountains and Maricao State Forest. The ASA for this species encompasses a small portion of the Cayey Mountains and several small localities along the Central Cordillera.
A large portion of ASA is outside of protected areas in Puerto Rico (≥ 80%) for birds, frogs, and birds and frogs combined (Table 4). Despite the large extension of ASA outside protected areas, approximately 85% of protected areas overlap within ASA for at least one species (Table 5). These results show that the existing protected area network 1 3 is located in appropriate areas for species protection, but the extent of these areas captures < 20% of the ASA. Furthermore, the protected area network captures a greater percent of the ASA of birds (78%) than for frogs (39%). In general, the relationship between the accumulative number of overlapping distribution of species ASA and protected areas was the same for birds and frogs (Fig. 3). As the number of overlapping ASA species distribution increased there was a dramatic reduction in the extent of the area occupied, meaning that areas with a high number of overlapping species ASA distributions are relatively small compared with areas with just a few  Fig. 3 The relationship between the accumulative number of overlapping ASA species distributions and A-C extent (km 2 ) of the overlapping area and D-F percent of the overlapping distribution that occurs within protected areas for A and D birds, B and E frogs, and C and F birds and frogs combined 1 3 overlapping species ASA distributions. In contrast, the area with the greatest number of overlapping species had the highest level of protection. The area of overlap of ASA for more than six species was ~ 96 km 2 for frogs and 388 km 2 for birds. Areas with a high number of overlapping ASA for frog species (> 6 species) were restricted to few sites in the east (Luquillo Mountains), southeast (Cayey Mountains), and the Central Cordillera (Toro Negro) (Fig. 4). In contrast, areas with a high number of overlapping ASA for bird species were restricted to the west (Maricao State Forest), northwest (Karst Mountains), and to a few sites in the east (Luquillo Mountains) (Fig. 4).
The Maricao State Forest appears as an always-suitable area for the majority of bird species (12 species), followed by the northwest Karst mountains (9 species) and the Luquillo Mountains (7 species) (Fig. 4). In contrast, the Luquillo Mountains appear as an alwayssuitable area for the majority of frog species (8 species), followed by a portion of Toro Negro (6 species) and Cayey Mountains (6 species) (Fig. 4).

Discussion
Our findings indicate that existing protected areas in Puerto Rico will not suffice to safeguard Species of Greatest Conservation Need in the context of projected climatic changes that could occur by mid-century under a very high greenhouse gas emissions scenario. A larger portion of the always-suitable areas for several species are outside protected areas, which may compromise the long-term survival of these species. In addition, protected areas are relatively small and sparse, and this can hinder species movement and increase local extinction risk. Therefore, the establishment of new protected areas, buffer zones, and natural corridors can substantially improve the likelihood of species conservation and persistence even in an 'extreme' mid-century climate scenario. By combining species distribution projections from different time periods (i.e., historical, current, and future) and identifying the always-suitable area (i.e., ASA) for each species, we can assess habitat stability for species and identify and prioritize areas that are more likely to promote long-term survival of multiple species to mitigate the impacts of climate change.
Despite more than 80% of the always-suitable areas for several species are outside protected areas, three protected areas (El Yunque National Forest, Maricao State Forest and Carite State Forest) had high ASA overlap. Although these areas may play an important role in preserving bird and frog populations, they are currently isolated from other protected areas and land use intensity in the surrounding matrix may not support the wellbeing of remaining bird and frog population inside the protected areas. In addition, the extension of these areas captures < 20% of the ASA.
Our study and others focused on the effects of climate change on species distributions (Lemes et al. 2014;Noce et al. 2017;Vasconcelos et al. 2018;Báez et al. 2019;Velazco et al. 2019) are limited by the assumption that the current species distribution is at equilibrium with the environment. In contrast, the current distribution may reflect past biotic (e.g., competition, predation, diseases, dispersal) or anthropogenic (e.g., urbanization, habitat degradation) effects, (Nogués-Bravo 2009;Boulangeat et al. 2012;Maguire et al. 2015). In addition, species niches are not static, and current distributions may reflect the adaptation of species to past and current climate conditions (Nogués-Bravo 2009). There is also much discussion about using binarization to determine suitable and unsuitable areas. Although some studies have advocated against this approach (Guillera-Arroita et al. 2015;Santini et al. 2020), there are few alternative solutions. To address these limitations, we provide a conservative approach that requires suitable areas to agree in the three time periods and include several species from two taxonomic groups.
Another limitation of our study is a sampling gap in part of the Central Mountains (Fig. 1). Nevertheless, our work provides the most comprehensive effort in acquiring detection/non-detection data using a standardized methodology, and current species distribution projections are in agreement with our experience in the field. Furthermore, there are few species records for this region in other databases (e.g., Global Biodiversity Information Facility-GBIF).
Overall, the distributions of individual frog species were more limited than those of the bird species, and consequently, patterns of occurrence in the protected areas and ASA are different for frogs and birds. Frogs had lower representativeness in the protected areas than bird species, a pattern also observed at the global scale (Rodrigues et al. 2004), and unfortunately most of the always-suitable areas for both frogs and bird species occurred outside of protected areas. Always-suitable areas for frogs are mostly concentrated in the eastern and central mountains of the island, which are more humid and receive more precipitation than the western part of the island (Daly et al. 2003). In contrast, always-suitable areas for birds are mainly concentrated in the western part of the island. The limited distribution of the frog species associated with a small overlap with protected areas and poor dispersion abilities may increase frog species vulnerability to extinction due to natural and human disturbances. Given that almost all frog species in Puerto Rico are already under a threat category (IUCN 2020), expansion of protected areas to incorporate the present and future distributions of these species may be necessary to ensure continued persistence.
Although birds have higher dispersal abilities than frogs, increasing connectivity between protected areas may also benefit bird populations. For instance, two endemic and 1 3 forest specialist bird species (Nesospingus speculiferus and Setophaga angelae) currently occur only in a few isolated populations in montane forests suggesting either poor dispersal abilities or specific habitat requirements. Expanding protected areas and connectivity can also mitigate the effects of natural disturbances such as hurricanes and droughts on fauna (Campos-Cerqueira and Aide 2021; Wunderle 1995) by providing a better representation of diverse landscapes. Alternatively, and complementary, the development of buffer zones around protected areas may offer a cheaper solution to increase connectivity. Several studies have shown that an increase in forested coverage around protected areas can increase their capacity to conserve viable populations (DeFries et al. 2005;Hull et al. 2011;Xu et al. 2016).
Most of the always-suitable areas are in higher elevations (> 400 m), with a clear spatial west-east pattern along the central cordillera. It is not a surprise that most of the overlap of always-suitable areas occur in regions at high elevations because these areas coincide with forest refugia during the intense agricultural phase that occurred in Puerto Rico between 1800 and 1940. In addition, many species are moving toward high elevation areas in Puerto Rico (Campos-Cerqueira and Aide 2017; Campos-Cerqueira et al. 2017) and other regions of the world (Chen et al. 2011) due to climate change. Therefore, given the essential role that high elevation areas play in species distribution any strategic plan to improve the conservation of SCGN will need to focus on the protection of areas around the Central Cordillera.
Secondary forest with different ages now covers more than 50% of the island. However, these new forests are still spatially fragmented and mainly concentrated in mountainous regions, whereas high-density urban areas and the few remaining agropastoral activities remain concentrated in the lowlands (Gould et al. 2012). Despite having suitable climatic conditions for the occurrence of species in the long-term, the ability of these secondary forest to provide essential resources (e.g., nesting cavities) to long-term species occurrence remain to be tested. Although forest cover has remained stable since 2004, the human population is still declining, and the likelihood that forest cover will increase is high. In this optimistic scenario, always-suitable areas outside protected areas will remain in land cover states less modified by humans and will be able to provide stable habitats for SGCN, provided that these secondary forests allow long-term species occurrence. Nevertheless, a recent study shows a 90% increase in residential development around protected areas in Puerto Rico between 2000 and 2010 (Castro-Prieto et al. 2017), despite the human population decline, which may compromise the future of always-suitable areas outside protected areas. In this scenario, the establishment of new protected areas in locations with a high number of overlapping ASA from species of great conservation need may be the best conservation strategy. Alternatively, and complementarily, assisted migration, captive breeding, and reintroduction may also be part of future conservation planning (Loss et al. 2011), especially with species with low dispersal abilities like frogs.
The approach presented in this study to assess the effectiveness of protected areas to promote species persistence can be easily applied in other parts of the world to identify always-suitable areas for several taxonomic groups and evaluate the role of protected areas in the context of climate change: (1) historical, current and future climate data from around the world are currently available online (e.g., Worldclim; Fick and Hijmans 2017), (2) audio recorders are more accessible and affordable (e.g. Audiomoth; Hill et al. 2019), and (3) accessible analytical platforms to store, manage and analyses audio recordings and extract data for species distribution models already exist (i.e., RFCx-ARBIMON; LeBien et al. 2020). Furthermore, several user-friendly R packages are available to create species distribution models. Given the speed of projected climatic changes, the variability in species-level responses, and the fact that only 12% of the Earth's surface (and 16% in Puerto Rico;Castro-Prieto et al. 2019) is currently under some form of legal protection (Jenkins and Joppa 2009;Joppa and Pfaff 2011), effective conservation actions will also likely require efforts to systematically improve the extent, frequency, and duration of longterm monitoring of fauna.

Conclusions
Some protected areas in Puerto Rico provide always-suitable habitats for many bird and frog species, helping to conserve these species in a future, which may include extreme climate change. Nevertheless, long-term protection of these species may not be successful given that (1) regions with a high number of overlapping ASA species are relatively small, (2) protected areas are small and sparse (Castro et al. 2019), and (3) large portion (~ 80%) of always-suitable areas is currently outside protected areas. Focusing special attention to the large proportion of always suitable areas currently outside protected areas may be effective in planning investments for regional conservation. To our knowledge, no new protected areas have been established in Puerto Rico to address the impact of climate on species distributions; nevertheless, we already have all the tools to identify and prioritize areas that are more likely to safeguard species in a global change scenario.