TMCF are defined as areas covered by tropical evergreen forests, without dry seasons, with annual precipitation up to 10,000 mm. The elevational range of the TMCF varies from 800 to 3500 m.a.s.l. depending on its location along the Andean Mountain ridge in Ecuador (Bruijnzeel et al., 2011). Thirty-two species from three taxonomic groups – ten birds, five amphibians and seventeen vascular plants (Table 1) – ecologically associated and geographically related with the eastern TMCF side of the Ecuadorian Andes were selected to be included in the ENM. These taxa were chosen because: (1) birds are a well-known group that contains the largest numbers of occurrence points per species, therefore containing most biogeographical information; (2) amphibians are ecological indicators of the status of the ecosystem or environmental stress due to their sensitivity to environmental changes (Beebee and Griffiths, 2005; Blaustein et al., 1994); and (3) vascular plants are often used to define vegetation types of an ecosystem (Blank and Blaustein, 2012; Prieto-Torres and Rojas-Soto, 2016; Prieto-Torres et al., 2016). Moreover, the selected species span a variety of life forms, and trophic guilds, are phylogenetically unrelated, locally endemic to TMCFs and have dissimilar geographical distribution patterns throughout the Neotropics. Only TMCF specialists that had most occurrences within an altitudinal range between 1300 and 2900 m.a.s.l. were selected. The number of species was defined based on previous studies indicating that a minimum of 20 species (Rojas-Soto et al., 2012), with more than ten occurrence records each (Pearson et al., 2007), to be able to get adequate ecosystem reconstructions.

Species occurrences records were gathered from the Ecuadorian National Program for Biodiversity Monitoring (SINMBIO) database, a consortium of various institutions and universities of Ecuador (https://bndb.sisbioecuador.bio/bndb/collections/index.php) and from the Global Biodiversity Information Facility database (GBIF; www.gbif.org). We retrieved GBIF specimen records for plants and amphibians, and specimen and observation records for birds; additional species information on the species status was retrieved from the IUCN Red List (IUCN, 2021) and the Plant list (Alkin et al., 2013). Each species occurrence record was georeferenced to correct for imprecise geographic coordinates using the MapLink tool (www.maplink.com). Coordinates that were duplicated, inconsistent or from doubtful localities were eliminated. A total of 4373 unique occurrence data points (available at https://github.com/gamamo/CloudForests) were used to perform ENM using MaxEnt (version 3.3.3k) entropy algorithm (Elith et al., 2006, 2011; Phillips et al., 2006).

MaxEnt estimates the distribution probability by extracting a sample of background locations (i.e., climatic variables) that contrasts with the known species occurrences, performing iterations until a convergence limit is reached. This process results in a species suitability map ranging from 0 (unsuitable) to 1 (suitable) (Merow et al., 2013). Each species distribution was limited to its specific accessibility area. For each species, it was delimited by known historical and geographical barriers, as found in the literature and based on the locality of occurrence records (Barve et al., 2011; Peterson, 2011; Peterson and Soberón, 2012). The accessibility area boundaries were based on (1) the classification of terrestrial ecoregions (Olson et al., 2001); (2) the map of Biogeographical provinces of the Neotropics (Morrone, 2014) and (3) the altitudinal range limits, representing historical barriers to dispersal and tolerance limits (Barve et al., 2011; Soberón and Peterson, 2005). Nineteen climate layers from the CHELSA project (Karger et al., 2017) and the evapotranspiration layer from ESRI (Zomer et al., 2006) were used. These layers have a resolution of 30” (ca. 1 km2 at the Equator). Only low-correlated variables (−0.8 < r > 0.8) for each species’ model, according to the Jackknife test (all parameters at default) calculated by MaxEnt (Royle et al., 2012) were kept in the models to ensure reduction of over-fitting of the generated distribution models (Ortega-Andrade et al., 2015).

Species models were generated using a random sampling of 80% as training data and 20% as testing data. In addition, 10,000 iterations with no extrapolation to avoid false niche optima, with clamping (for the projection of future scenarios) were set. A total of five replicates per species were generated. A logistic response was chosen to get the values for habitat suitability (continuous probability from 0 to 1). The 10-percentile training presence threshold was chosen for the binary presence-absence maps. This threshold identifies suitable pixels that are predicted to have similar suitability as the pixels that contain the species occurrence record, as result of this rejecting 10% of the records to minimize over prediction caused by outliers (Jiménez-Valverde and Lobo, 2007; Prieto-Torres and Rojas-Soto, 2016); All other parameters were maintained at default settings. Model performance was evaluated by calculating commission and omission errors, the Area Under the Curve (AUC) of the Receiver Operating Characteristic (ROC) (Elith et al., 2011), as well by applying the partial ROC curves test. The later avoids an inappropriate weighing of the omission and commission components in the AUC analysis (Lobo et al., 2008; Peterson and Soberón, 2008). Partial AUCs were calculated using theNicheToolBox2 (http://shiny.conabio.gob.mx:3838/nichetoolb2/; Osorio-Olvera et al., 2020), with 50% of the unique occurrence testing data points, and represent the ratio of the AUC model to the null expectation (“AUC ratio”). We evaluated the statistical significance of AUCs compared with null expectations by bootstrap resampling 50% random points, replacing values 1000 times from the overall pool of data iterations. Significance (as compared with null expectation) was assessed by ranking the observed AUC ratio with the values from pseudo-replicates, following the proposal of Peterson et al. (2008).

Future Global Climate Models (GCMs) for the Ecuadorian TMCF were selected based on two Representative Concentration Pathway scenarios: RCP 4.5 is an optimistic response to stabilizing the anthropogenic components of radiative forcing; RCP 8.5 can be considered a pessimistic outcome with high levels of greenhouse gases (Thomson et al., 2011). Future scenarios were set to the year 2050 since TMCF faces conservation issues that threatens the persistence of the ecosystem in the short term. Two GCMs models were chosen: The ACCESS 1.0, developed by the Commonwealth Scientific and Industrial Research Organization and Bureau of Meteorology, and the HadGEM2-ES, developed by the Met Office Hadley Centre and Brazilian Institute for Space Research. The GCMs were downloaded from the CHELSA website (http://chelsa-climate.org/future/) as digital layers to be used in the ENM. The ACCESS 1.0 was chosen based on previous assessments of climate change on community-level species (Prieto-Torres et al., 2016), while the HadGEM2-ES was chosen as it provides better simulations of annual and seasonal precipitation in the Amazon compared to other GCMs (Sierra et al., 2015; Yin et al., 2013).

The reconstruction of the TMCF ecosystem was done by overlapping individual species predictive maps, which was subsequently constrained to the political limits of Ecuador. To analyze the geographical patterns of TMCF, the sum of species per 1 km pixels (from 1 to 32 species – independently of the species’ identity) were performed according to the methodology by Prieto-Torres and Rojas-Soto (2016). The numbers indicate the likelihood of TMCF to be present in this pixel. Higher number of species indicate a higher likelihood of being TMCF. Also, the mean range, elevation and standard deviation were calculated for each combination. This step was done independently for each of the five climatic models (from here on: Present, ACC45, ACC85, HGE45 and HGE85).

The current distribution of the TMCF in eastern Ecuador was evaluated by identifying areas that contain overlap of species distributions in the present model and thus higher certainty of TMCF presence. To assess the shifts of the TMCF area under each future climatic model, we calculated the proportional difference in the total number of pixels for each species combination between present and future models. The proportional decline of TMCF per species was calculated by dividing the TMCF area obtained from future models by the TMCF area obtained in the present model. The mean elevation and range were calculated from the total of each number of species present in the pixels (1–32). Elevation range differences (future models elevation range minus present model elevation range) and proportional elevation change (future models elevation range/present model elevation range) were also calculated. An ANOVA followed by a post-hoc Tukey test was used to detect significant differences between TMCF area and elevation between the future models.

We compared the overlapping differences between the predicted TMCF area and the distribution limits of the adjacent ecosystems, as defined by the Natural Regions of Ecuador map (Ron, 2020), per sum of species, to identify transitional areas. Furthermore, the proportion of predicted area that coincides with each ecosystem and the contribution of the overlap area to the total ecosystem size was calculated. To analyze the correlation of the modelled TMCF with the single-ecosystem TMCF map (ENE) defined by the Ministry of Environment of Ecuador (MAE, 2013), we calculated the overlap of pixels for each combination of species with the single-ecosystem map. The overlap of the TMCF distribution with the Ecuadorian protected areas (PA) was evaluated using the layer of protected areas downloaded from ProtectedPlanet.net (IUCN, 2012). We calculated the pixel overlap per sum of species and the proportion of pixels within the PAs. All analyses were performed using ArcMap 10.5 and RStudio (https://www.rstudio.com/).

Our future climate models showed that the current TMCF area would be pressured upon to higher elevations and that suitable areas for the individual species may largely decline. A possible reduction of 42–54% in TMCF area and an elevational shift ranging between 207 m and 429 m by 2050 were predicted. Our results suggest that mountain species will undergo drastic distributional changes, being pushed upwards to track their climatic niches as temperatures increase (Forero-Medina et al., 2011; Nogués-Bravo et al., 2007; Ortega-Andrade et al., 2015). Some species, however, may disappear due to the lack of favorable environmental conditions (Nogués-Bravo et al., 2007; Ortega-Andrade et al., 2015). Our predictions are consistent with the directional shift of TMCF species in Mexico (Peterson et al., 2002; Rojas-Soto et al., 2012), confirming the altitudinal movement of the ecological niches due to climate change impacts (Guevara, 2020; Forero-Medina et al., 2011; Nogués-Bravo et al., 2007) and that lowlands species may slowly become more frequent in higher elevations, while the environmental suitability areas of mountain species may decline or disappear (Bubb et al., 2004; Feeley et al., 2011; Forero-Medina et al., 2011; Nogués-Bravo et al., 2007).

The predicted movement of species towards higher elevations may lead to a redefinition of the spatial location of transition zones between montane ecosystems. If ecosystem boundaries can be defined based on species distributions, it is expected that these boundaries may shift over time and will need to be redefined whenever species distributions shift to track climate change (Guevara, 2020). For instance, in our future climate models, the TMCF that overlap with the paramos nearly doubled, while the lower limits showed a reduction of nearly almost half of its current area, which might indicate the increase of the TMCF in higher montane areas. This result coincide with the predictions of Helmer et al. (2019), that 70%–86% of the paramo in the Neotropics will be likely subjected to tree invasion from the adjacent ecosystems under climate change scenario. Lower parts of the TMCF area in turn, are predicted to be replaced by lower altitude ecosystems (Bubb et al., 2004).

The current PA areas may preserve roughly 33%–40% of the TMCF extent under future climatic scenarios, in contrast to the 8% that is currently under protection. The increase of TMCF in PA can be explained by the fact that several PA were created to protect the higher elevation ecosystems, such as the paramos. As the TMCF boundaries move towards higher elevations due to the altitudinal movement of suitability areas to species, the relative TMCF area can potentially replace current protected paramo areas. Thus, a corridor of PA along elevational gradients in montane areas could effectively minimize the negative effects of short-term climate change for montane ecosystems and transitional zones. Almost all the natural forests of the central Andean valley have been removed and only 4% remain intact on the west Andean slope (Aldrich et al., 1997). The synergistic effect of deforestation and agricultural expansion with climate change may result in short-term negative impacts in TMCF, especially along the southern part of Ecuador, where few PAs are located but TMCF area is extensive (Lessmann et al., 2014; Ortega-Andrade et al., 2015; Prieto-Torres et al., 2016). Because conservation efforts should also target also transitional borders with adjacent ecosystems and future climatic conditions, a reconsideration of the TMCF’s boundaries and the recognition of dynamic species movements are necessary to implement efficient conservation goals to this ecosystem.