Climate change and tree cover loss affect the habitat suitability of Cedrela angustifolia: evaluating climate vulnerability and conservation in Andean montane forests

Study area

The study area covered the Andean montane forest region (Ecuador, Peru, Bolivia, and Argentina; Fig. 1), with elevations ranging from ≈1,800 to 3,300 m asl, as recorded by Pennington & Muellner (2010). The study was delimited by Tungurahua Province, Ecuador, in the north (1°S), and Catamarca Province, Argentina, in the south (28°S). The Servicio Nacional Forestal y de Fauna Silvestre-SERFOR approved the research outside Protected Natural Areas through the General Management Resolution RDG 007-2020-MINAGRI-SERFOR-DGGSPFFS.

Data sampling

We obtained occurrence data for Cedrela angustifolia from the Tropicos database (https://tropicos.org/home; Missouri Botanical Garden, 2022), the Global Biodiversity Information Facility database (https://gbif.org/; GBIF Secretariat, 2022), and scientific publications (i.e., Inza et al., 2012; Paredes-Villanueva, López & Navarro Cerrillo, 2016; Pennington & Muellner, 2010; Wong & Reynel, 2021) totaling 310 occurrences. In addition, we checked and excluded incorrectly georeferenced points, outliers and herbarium accessions with missing location data. We also used geographical distances to narrow our datasets. We clustered occurrences to a grid size of ∼one km² based on the environmental heterogeneity of the ecosystems in which C. angustifolia occurs, such as montane forests and inter-Andean valleys, to filter out closely duplicated records. Finally, 104 records were obtained from Ecuador, Peru, Bolivia, and Argentina (89 GBIF and Tropicos.org and 15 scientific studies).

Data preparation

We used 39 environmental variables: 22 bioclimatic variables for the present and future periods, elevation, 11 soil raster layers, one NDVI raster, and four forest change data rasters (Table S1). Under the assumption that it influences niche conservatism (Peterson et al., 2011), the climate sensitivity of C. angustifolia is related to its present and future distributions (2011–2040, 2041–2070, and 2071–2100 years). We obtained raster model bioclimatic data from the CHELSA database with a spatial resolution of 1.0 km² (https://chelsa-climate.org/; Karger et al., 2022) for the present and future periods. In the present and future models, we included 19 raster model bioclimatic variables and potential evapotranspiration (PET), climate moisture index (CMI), and near-surface relative humidity (HURS) (Table 1).

We used ‘Coupled Model Intercomparison Project Phase 6’ (CMIP6), which is part of the ‘Working Group on Coupled Modelling’ (WGCM), to assess projections of bioclimatic variables to three different times in the future (Eyring et al., 2016). We considered the model projections because the models are given a common set of future concentrations of greenhouse gases, aerosols, and other climate forcing to project what might happen in the future (Eyring et al., 2016).

We selected two future climate scenarios (shared socioeconomic pathways; SSP 3-7.0 and 5-8.5) for the five models to derive future climate projections (Table S1). We assume that SSP 3-7.0 is a less chaotic scenario because it represents a smaller reduction in global greenhouse gas concentrations by 2100 than SSP 5-8.5, which is the highest carbon emissions scenario and the most pessimistic view of the future. We considered a digital elevation model (DEM) variable (Karger et al., 2022) to assess the elevation effects.

In addition, we obtained 11 soil raster layers from SoilGrid (https://soilgrids.org/; Batjes, Ribeiro & Van Oostrum, 2020) with a spatial resolution of 250 m. Moreover, we obtained global forest change data (land cover, tree cover loss, loss and gain of forest raster layers from 2000 to 2022) from Global Forest Watch, with approximately 30 m per pixel (https://www.globalforestwatch.org/map/; Hansen et al., 2013). Furthermore, we used the Normalized Difference Vegetation Index (NDVI), which indicates the vegetation coverage and spatial distribution of vegetation, calculated as the annual mean NDVI for 2000–2020, using Earth Explorer data with one km² spatial resolution (https://earthexplorer.usgs.gov/).

For present and future models, we assumed that the forest loss, soil and NDVI layers were constant to determine the climate effect on vegetation loss and soil characteristics. Finally, we used the resample function so that all environmental variables had a resolution of ∼one km², so we used the bioclimatic variables as reference values (Y) and the soil, topography, and forest change layers as the layers to be resampled (X) considering a bilinear interpolation, which were included in the raster package (Hijmans, 2023) in R software v. 4.3.1 (R Core Team, 2022).

Data analysis

We determined the appropriate selection of environmental variables to perform statistical and ecological parameter analyses of suitable models, according to Ames-Martínez et al. (2022) (Fig. 2). We selected variables with multicollinearity based on measures, such as principal component analysis (PCA), to identify climatic variables potentially important in the geographic distribution of C. angustifolia. We utilized the entire environmental matrix and we selected the variables that were present in the first two principal components (Estrada-Peña et al., 2013) with eigenvalues > 1 and percentage of variance explained 53% and 31% of the variation respectively, for a total of 84%. In addition, we assessed the collinearity of the covariates with the variance inflation factor (VIF) in the pre-selected matrix of the PCA analysis, we sequentially excluded the one with the greatest VIF sequentially until all remaining predictors had VIFs < 10 (Cobos et al., 2019a). Furthermore, we selected variables with pairwise Pearson’s correlation < 0.7 from the VIF analysis matrix to be used for the species distribution model. Our approach involves the use of sdm (calibration function; Naimi & Araújo, 2016), fuzzySim (corSelect function; Barbosa, 2015), regclass (VIF function; Petrie, 2020), FactoMineR (PCA function; Husson et al., 2023), and virtual species packages (removeCollinearity function; Leroy et al., 2019).

Finally, after variable selection, we selected 16 environmental variables: isothermality (BIO3), temperature seasonality (BIO4), minimum temperature of the coldest month (BIO6), mean temperature of the warmest quarter (BIO10), precipitation seasonality (BIO15), precipitation of the warmest quarter (BIO18), precipitation of the coldest quarter (BIO19), climate moisture index (CMI), altitude (ALT), soil organic carbon stock (SOCS), organic carbon density (OCD), silt content (SILT), clay content (CC), deforestation (DEF), forest loss (LF), and NDVI. For the analysis of tree cover loss, we used the forest loss raster in some models and other models without this raster to compare forest loss in the present and future periods, assuming that the same current rate of forest loss will continue into the future.

Modeling and validation

Comprehensive evaluation

We derived the potential distribution from the best model and assessed the relative importance of each variable using the average performance evaluation indicators corresponding to the ‘area under curve’ (selected models with range of 0.97 < AUC < 1), and the partial receiver operating characteristic (partial ROC, selected models with a range of 0.90 < ROC < 1); to determine which of multiple models is most likely to be the best model we used Akaike Information Criterion corrected (AICc; as optimal complexity parameter), and Bayesian Information Criteria (BIC) (Hirzel et al., 2006; Peterson & Soberón, 2012; Jiménez & Soberón, 2020). A partial ROC was used on 50% occurrence data for bootstrap resampling, 100 iterations, and omission rate error (5%, maximum permissible omission error) (Cobos et al., 2019b). Success rate curve were further used to assess the performance of the MaxEnt and RF models in predicting species distributions for model validation (Rahmati, Pourghasemi & Melesse, 2016).

Finally, we used Schoener’s D index (Warren, Glor & Turelli, 2008) to compare the similarity of suitable distribution maps between the MaxEnt and RF models using the ENMTools package (Warren et al., 2021); this method calculates the random permutation by 100 times the occurrence of C. angustifolia in both models. If Schoener’s D index is lower (0, no overlap) or higher (1, complete overlap) than 95% of the simulated D values, both models were more dissimilar (or similar) from the suitable distribution than expected by chance (Sillero, Ribeiro-Silva & Arenas-Castro, 2022). Maps were averaged to generate spatial information on present and future presence probabilities across all C. angustifolia. We followed the ODMAP protocol for the modeling process (Overview, Data, Model, Assessment, and Prediction; Table S1; Zurell et al., 2020).

We determined the surface area variation in each climate model (km²) between present and future scenarios (Table 1). The final shapes were obtained using MaxEnt with maps processed using QGIS v. 3.18.3 (QGIS.org, 2021). We compared the variation in mean temperature, annual precipitation, and tree loss variation between the present and three future periods using raincloud plots from the packages ggdist (Kay & Wiernik, 2023), gghalves (Tiedemann, 2022), and ggplot2 (Wickham et al., 2021). Two-way ANOVA and post-hoc Tukey’s test were used to compare the means of these variables using the rstatix package (Kassambara, 2023).

Predicted refugges to climate change

We concatenated the present and future potential of Cedrela angustifolia distribution with predicted refuges to climate change, recognizing appropriate grids in scenarios SSP 3-7.0 and SSP 5-8.5. We used the consensus model between the present and future models, land cover (Hansen et al., 2013), and protected natural areas (PNA; IUCN, 2023) to recognize and estimate remnant patches outside the PNA. We performed a spatial distribution bias correction to avoid over-adjusting of future projections. We included 10,000 bias files and bioclimatic variables to assess potential refugees for climate change analysis. We implemented a Gaussian kernel analysis using the kernel function with QGIS software to avoid sampling bias and identify the highest potential refuges for climate change.

Model evaluation and contribution of predictor variables

We generated and compared all candidate MaxEnt and RF models and selected one model from each period that met the criteria of significance, predictive ability, fit and complexity (Table S2). Our results showed excellent performance for both the MaxEnt and RF models when assessed against the independent test dataset. Nonetheless, the MaxEnt model demonstrated a higher AUC ratio; however, RF exhibited better predictive performance than MaxEnt (Table S2).

For the present and future models, the relative contribution of each predictor variable to the SDMs was assessed by visualizing the percentage contribution and permutation importance (Table 2). In our analysis of the present and future models, we identified ten environmental variables as the most important factors in the model fit. These variables were precipitation seasonality (BIO15), soil organic carbon stock (SOCS), normalized difference vegetation index (NDVI), temperature seasonality (BIO4), organic carbon density (OCD), precipitation of the warmest quarter (BIO18), silt content (SILT), clay content (CC), loss forest (LF), and isothermality (BIO3). These variables showed a high percentage contribution to the model for both the present and future scenarios, as well as for the three periods (Table 2).

Present potential distribution

The potential distribution with the tree cover loss effect of C. angustifolia was approximately 13,080 km² (mean of MaxEnt and Random Forest values), and Schoener’s D index between the RF and MaxEnt models was 0.857. Nevertheless, excluding the effect of tree cover loss, the total distribution was 16,148.5 km², with a Schoener’s D index of 0.749. Notably, the areas in Peru and Bolivia were suitable for C. angustifolia according to the current records (Fig. 3).

With the loss forest effect, the distribution detected was 798 km² in Ecuador, 1,947 km²in Peru, 9,591 km² in Bolivia, and 744 km² in Argentina. Notwithstanding, without loss forest effect, we detected that in Ecuador exhibited 948.5 km², Peru with 2,579 km², Bolivia with 11,400.5 km², and Argentina 1,220.5 km² (Fig. 2).

We found that the presence of C. angustifolia was mainly influenced by these environmental variables: BIO3 (from 5.3 to 6 °C), BIO4 (20 °C to 31 °C), BIO18 (100 to 200 mm), BIO15 (750 to 830 mm), SOCS (0–5%), NDVI (0.95−1.00 units), OCD (0–5%), slit (0–10%), CCF (0–10%), LF (30–55%).

Future potential distribution

Under two scenarios, we detected a decrease in the distribution range of SSP 3-7.0 and 5-8.5 during the three periods (Fig. 3). For 2040, we estimated an extension equivalent to 8,934 km² (SSP 3-7.0) and 9,094 km² (SSP 5-8.5), with Schoener’s D index of 0.759 (Fig. 4A). Nevertheless, without the effect of tree cover loss, we detected 11,424 km² (SSP 3-7.0) and 11,078 km² (SSP 5-8.5), indicating 29.26% and 31.40% loss forest, respectively, with a Schoener’s D index of 0.786 (Fig. 4B).

For 2070, our models predicted a decrease in area from 29.91% (SSP 3-7.0) to 30.53% (SSP 5-8.5) without the loss forest effect, and from 33.93% (SSP 3-7.0) to 42.32% (SSP 5-8.5) with the influence of tree-cover loss, and Schoener’s D index was 0.846. Finally, for 2100, the area decreased by 30.43% (SSP 3-7.0) to 33.33% (SSP 5-8.5) without the loss forest effect; however, the loss forest effect decreased by 27.40% (SSP 3-7.0) and 38.47% (SSP 5-8.5) in the whole distribution, with Schoener’s D index of 0.872 (Figs. 4A, 4B).

The mean temperatures up to 2040 in the two scenarios showed no significant differences among the countries (Figs. 5A–5D). However, Argentina and Bolivia exhibit statistically significant differences by 2070. Similarly, both countries displayed similar annual mean temperature values (∼20.5 °C and ∼18.7 °C, respectively) for 2040 (Figs. 5A–5B). Ecuador and Peru showed similar annual mean temperature values (∼15.2 °C and 16.3 °C, respectively) between the present and 2040 (Figs. 5C–5D). In the four countries, 2100 presented statistically significant differences compared to the other periods (Figs. 5A–5D). In contrast, the annual precipitation showed no significant differences among the four periods or countries (Figs. 5E–5H).

Effect of forest cover loss in the present and future

The relationships between the predictor and response variables show how forest cover loss affect model predictions. Therefore, we further analyzed how the predicted probability of species occurrence changed with the forest cover loss effect using the marginal responses of the probability of presence suitability (Fig. 6A). We found that 30.59% of the total area was the most suitable habitat without the forest loss; nevertheless, Peru, Ecuador, and Argentina showed variations in the gain of forest loss area in 2040, 2070, and 2100 (Fig. 6B). For example, Ecuador and Argentina decreased by 51.25% and 14.08%, respectively, under the influence of the forest loss effect (Fig. 6B) or the total distribution, and Bolivia increased by more than 30% of the area distribution in all periods, without the forest loss effect. Nevertheless, Peru and Argentina are expected to decrease by 70% by 2070.

Potential refuges to climate change

Our analysis revealed variations in the habitat suitability distribution areas of C. angustifolia (9,449 km²) across the four countries. Nonetheless, Ecuador (724 km²), Peru (1,784 km²), and Argentina (683 km²) exhibited a high potential for refuge from climate change; whereas Bolivia (6,258 km²) displayed a decreasing habitat (Fig. 7). Only 24.28% of the current potential distribution is within protected areas, and this is expected to decrease to 25–30% by 2100.

Predictions of the core distribution regions from the MaxEnt and RF models showed a higher heterogeneity and stronger gradients. In addition, our analysis identified Azuay and Zamora Cinchipe as refugees in Ecuador; Cajamarca, Junín, Apurimac, and Cusco as refugees in Peru; Cochabamba, Chuquisaca, and Tarija as refugees in Bolivia; and Salta and Tucumán in Argentina as potential refugees to climate change.

Present potential climate sensitivity

Previous studies have systematically evaluated the predictive effectiveness of RF and MaxEnt models through automated parameter optimization (Cotrina et al., 2021; Mi et al., 2017; Zhao et al., 2022). Although the RF models typically provide robust and accurate predictions using default configurations (Freeman et al., 2015), MaxEnt models often require parameter refinement (Feng et al., 2019; Jiménez & Soberón, 2020). We selected the best feature classes and regularization parameters for MaxEnt and compared them with the RF model. Consistent with previous literature, both the RF and optimized MaxEnt models exhibited commendable predictions accuracies; however, the RF model showed a marginal superiority, which was evident in both cross-validation and external dataset evaluations. This is consistent with the findings of Mi et al. (2017) and Čengić et al. (2020), who described the improved performance of RF over standard MaxEnt configurations for species distribution prediction.

The distribution of Cedrela angustifolia is influenced by seasonal variations in precipitation and temperature, and the model performance is excellent (AUC > 0.98), demonstrating the efficiency and accuracy of the model (Warren & Seifert, 2011). These populations exhibit tolerance to colder temperatures and higher moisture conditions, maintaining their evolutionary climatic conditions, as detected by Muellner et al. (2010) and Koecke et al. (2013). Our analysis suggests that Bolivian montane forests provide a suitable ecological assemblage (81% present model), and ecosystem conservation for C. angustifolia (Pennington & Muellner, 2010), with more than half of the records of this species within the NPAs.

Effect of climate change in the future

The distribution patterns of different species are influenced by different ecological and evolutionary factors that allow them to survive in specific environments within diverse landscapes (Rahbek et al., 2019). Despite this, Tejedor Garavito et al. (2015) argued that land use change, particularly deforestation of AMFs, would be more detrimental to biodiversity than climate change, leading to the loss of relict-endemic tree species worldwide (Feeley & Silman, 2010).

Based on the evaluated scenarios, a temperature increase of 4.1 °C to 5 °C is expected compared the present temperature record. This would result in a significant reduction in suitable habitat in Bolivia (>20.26%) and Argentina (>28.99%). By 2100, C. angustifolia will not be able to find optimal sites, because the thermal limit will be exceeded (∼21.8 °C). However, Peru and Ecuador offer favorable climatic conditions that would allow C. angustifolia to migrate to cooler areas and persist through rapid changes in climate (Pearson, 2006). This confirms the climatic impact of C. angustifolia (Cotrina et al., 2021; Koecke et al., 2013; Rodríguez-Ramírez et al., 2022). Similar results have been reported for C. odorata (Sampayo-Maldonado et al., 2023) and other Cedrela species (Cotrina et al., 2021; Koecke et al., 2013).

We demonstrated that tree cover loss affects more than 30% of the range of C. angustifolia, which is a critical factor contributing to the reduction in C. angustifolia in the four countries. If the current rates of deforestation continue or increase, this will lead to reduction in distribution (Hansen et al., 2013). Similarly, the loss of tree cover increases the risk of soil erosion, leading to reduced soil fertility and increased sedimentation in the water bodies in these forests (Bax & Francesconi, 2018; Cuenca, Arriagada & Echeverría, 2016; Tapia-Armijos et al., 2015). Therefore, the climatic impact of deforestation on AMF weakens the relationship between atmospheric circulation and the hydrological cycle (Longobardi et al., 2016).

Potential habitat suitability and refugees

Suitable habitat for C. angustifolia will be maintained in all four countries; however, habitat suitability will decrease in all countries in contrast to the present and future models evaluated. Bolivia exhibited more illegal logging and forest fires than the other three countries, making it more vulnerable to climate change (Hansen et al., 2013). In contrast, Ecuador, Peru, and Argentina had areas with better habitat suitability than the present model, suggesting that unexplored forests with similar climatic conditions allow the species to adapt through natural or active restoration (i.e., C. angustifolia plantations in inter-Andean valleys; SERFOR, 2020).

Furthermore, their presence increases in NPAs, as described by Pennington & Muellner (2010), which shows that NPAs function as biodiversity reserves and buffers against the effects of changing climatic conditions, allowing the formation of refugia and providing ecological corridors for species to acclimate or migrate over the long term (Cuesta, Peralvo & Valarezo, 2009; Geldmann et al., 2013). In contrast, 75.72% of habitat suitability was detected outside the NPAs, indicating the need to develop forest management and monitoring strategies to protect these forests as they are more susceptible to selective logging and timber overexploitation, as we demonstrated by our results (Cotrina et al., 2021; SERFOR, 2020).