ORIGINAL RESEARCH article

Front. Ecol. Evol., 24 September 2020

Sec. Biogeography and Macroecology

Volume 8 - 2020 | https://doi.org/10.3389/fevo.2020.565708

Climate Vulnerability Assessment of the Espeletia Complex on Páramo Sky Islands in the Northern Andes

  • 1. Facultad de Ingeniería y Administración, Universidad Nacional de Colombia, Palmira, Colombia

  • 2. The Alliance of Bioversity International and the International Center for Tropical Agriculture (CIAT), Rome, Italy

  • 3. Escuela de Estadística, Facultad de Ingeniería, Universidad del Valle, Cali, Colombia

  • 4. Departamento de Ciencias Biológicas, Universidad de los Andes, Bogotá, Colombia

  • 5. Jardín Botánico de Cartagena “Guillermo Piñeres”, Turbaco, Colombia

  • 6. Corporación Colombiana de Investigación Agropecuaria AGROSAVIA, C.I. La Selva, Rionegro, Colombia

  • 7. Departamento de Ciencias Forestales, Facultad de Ciencias Agrarias, Universidad Nacional de Colombia – Sede Medellín, Medellín, Colombia

Abstract

Some of the largest impacts of climate change are expected in the environmentally heterogeneous and species rich high mountain ecosystems. Among those, the Neotropical alpine grassland above the tree line (c. 2,800 m), known as Páramo, is the fastest evolving biodiversity hotspot on earth, and one of the most threatened. Yet, predicting climate responses of typically slow-growing, long-lived plant linages in this unique high mountain ecosystem remains challenging. Here we coupled climate sensitivity modeling and adaptive potential inferences to efficiently assess climate vulnerability of Espeletia, Páramo’s most iconic, predominant and rapidly evolving plant complex. In order to estimate climate sensitivity, we first modeled the distribution of 28 Espeletia taxa under a niche conservatism scenario using altitude and five current (1970–2000) and future (2050 RCP 8.5) bioclimatic variables across 36 different Páramo complexes in the northern Andes (49% of the world’s Páramo area). As an alternative to range shifts via migration, we also computed the adaptive capacity of these Páramo complexes by considering three enhancing factors of the biodiversity’s adaptive potential as well as three environmental limiting factors of the populations’ plastic response. These predictors showed that diverse Páramos in the Eastern Cordillera were more vulnerable likely because the counteracting effects of the adaptive potential (r = −0.93 ± 0.01) were not sufficient to buffer higher distribution losses (r = 0.39 ± 0.01). Agriculture (r = −0.48 ± 0.01), mining (r = −0.36 ± 0.01), and rural population density (r = −0.23 ± 0.01) also weakened the adaptive capacity. These results speak for a limited persistence via migration in the short-term responses of Espeletia to climate change, even though the past population dynamics in concert with glacial cycling is indicative of a predominant role of range shifts. Furthermore, changing climate, together with a general inability to adapt, may eventually constrain the rapid diversification in the Espeletia complex. Our integrative modeling illustrates how future climate may impact plant populations in a mega diverse and highly threatened ecosystem such as the Páramo, and encourages carrying out similar estimates in diverse plant complexes across other high mountain and island-like ecosystems.

Introduction

How tropical alpine biodiversity hotspots are affected and will react to climate change is among the most pressing questions in current biological research. Most studies on responses to changing climate assume populations migration to higher altitudes so that organisms would track their ecological niche (; ), and in some cases even to lower altitudes due to competitive release (). Yet, when migration potential is restricted (; ), the only way left for organisms to persist is by adjusting to the new environmental conditions, a traditionally neglected alternative in niche modeling. Adjustment via phenotypic plasticity is presumably important in long-lived taxa because plastic responses can happen within the lifetime of an individual (). However, plasticity may be hampered when populations face novel anthropic conditions not found in their evolutionary history (; ), e.g., agriculture, mining and population density. Organisms can also adjust to climate change through adaptation (). Still, populations may not have enough adaptive potential (e.g. limited diversity) to keep up with the pace of climate change (). Therefore, assessing biodiversity hotspots’ vulnerability to climate change not only requires tracing ecological niche gaps under future scenarios – i.e., sensitivity (), but also needs to explicitly account for the adaptive potential () and plasticity to anthropic pressures. These improvements may make predictions more congruent with empirical data, even in potential climate refugia ().

One of the fastest evolving biodiversity hotspots is the Páramo (), a highly threated (; ) alpine ecosystem dominated by diverse grasslands above the treeline in the American tropics (at elevations of ca. 2,800–5,000 m). It is a main water supplier of wetland ecosystems and densely populated areas in the northern Andes (; ). Despite it has a relatively small surface area (35,000 km2), it contains over 3,000 plant species, several of which are endemic (; ) and emerged as unique adaptations to extreme environments that evolved during the Andean uplift in the last five million years (; ; ). Unparalleled diversification rates across plant groups at these tropical “sky islands” () has further been attributed to colonization by pre-adapted linages (), and Pleistocene glacial cycling (; ; ) in the last 2.4 Myr that led to repeated periods of connectivity and spatial isolation (i.e., “species pump hypothesis”). Population fluctuations in concert with past glacial cycling indicate a general inability to adapt and a major role of range shifts (; ; ; ). Yet, it remains to be explored whether the rapidly evolving Páramo can keep pace with the quick rate of climate change and human expansion.

The most iconic, abundant and diverse group in the Páramo is the Espeletia complex (; ; ). It originated in the Venezuelan Andes () from where it colonized southwards the Colombian Eastern Cordillera and the Ecuadorian Andes (), followed by a northwards migration into the Colombian Central and Western Cordilleras (). Espeletia taxa (ca. 120) are ecologically abundant across the Páramo landscape (), occurring in highly heterogeneous habitats. Espeletia populations are adapted to grow in the wet depressions of high valleys, in the dry exposed slopes or even within the forests at the tree line, therefore experiencing a wide range of climatic conditions within elevations and localities (). This environmental heterogeneity may have important implications for the reaction of Espeletia to changing climatic conditions. For example, habitat variation may provide suitable locations for migrants within only a few meters of their current locations (). However, populations adapted to a narrow range of conditions may respond poorly to future threats. In other words, environmental heterogeneity effects on the responses to climate change may be mixed, and sometime contradictory, especially for a long living tropical alpine plant with limited seed dispersal, such as Espeletia. Hence, a much-needed research is to balance climate sensitivity under the current ecological niche heterogeneity with the adaptive capacity of Espeletia to assess the climate vulnerability of this plant complex.

In order to bridge this gap, our goals in this study were to (1) quantify climate sensitivity in the Espeletia complex by comparing current and future (2050 RCP 8.5) distributions under a niche conservatism scenario (i.e., range shift via migration), (2) compute the adaptive capacity by incorporating enhancing factors of the biodiversity’s adaptive potential (i.e., diversity, protected areas, and forest area in the buffer zone of the Páramo) as well as limiting factors of the plastic response (i.e., agricultural area, mining, and population density), and (3) estimate climate vulnerability of Espeletia by simultaneously coupling its climate sensitivity with its adaptive capacity. This work will help establishing how plant populations may respond to environmental change in a mega diverse tropical alpine region, where the largest impacts of climate pressures are expected (; ).

Materials and Methods

Study Area

The study area comprised the Colombian Andean Páramo biogeographical province (2,732–5,212 m, Figure 1), including the high-elevation northern Andes and the Sierra Nevada de Santa Marta, but excluding areas with similar altitudes in the Cordillera de Merida and Central American mountains (Figure 2). It was delimited in the north by the Sierra Nevada of Santa Marta (11°N), and in the south by the Colombia-Ecuador border (1 N). It spanned 36 Páramo complexes, as in and , denoting 49% of the world’s Páramos. The surveyed Páramo zones were covered in a 33% by protected areas1, and surrounded in a 30% by forest patches2. Meanwhile, agriculture and mining, respectively, disturbed 6 and 7%3 of these Páramos and their buffer belt (2 km around the tree line). The handling of these data sources is detailed below as part of the computation of the “Adaptive Capacity at the Study Area.”

FIGURE 1

for narrow-ranged species). Espeletia taxonomy is in accordance with .

FIGURE 2

) under a niche conservatism scenario (i.e., range shift via migration). Values close to zero in the sensitivity index (C) indicate few changes in the forecasted distribution of Espeletia due to climate change, while negative and positive values correspond to range losses and gains, respectively, based on the averaged presence probabilities. Five non-collinear WorldClim bioclimatic variables (Supplementary Table S3) were considered at a 30 s (1 km2) spatial resolution for the current and future scenarios, besides a Digital Elevation Model (DEM). Espeletia taxonomy follows . Páramo complexes are numbered in (A) according to and , specifically: (1) Perijá, (2) Jurisdicciones-Santurbán, (3) Tamá, (4) Almorzadero, (5) Yariguíes, (6) Cocuy, (7) Pisba, (8) Tota-Bijagual-Mamapacha, (9) Guantiva-Rusia, (10) Iguaque-Merchán, (11) Guerrero, (12) Rabanal y río Bogotá, (13) Chingaza, (14) Cruz Verde-Sumapaz, (15) Los Picachos, (16) Miraflores, (17) Belmira, (18) Nevados, (19) Chilí-Barragán, (20) Las Hermosas, (21) Nevado del Huila-Moras, (22) Guanacas-Puracé-Coconuco, (23) Sotará, (24) Doña Juana-Chimayoy, (25) La Cocha-Patascoy, (26) Chiles-Cumbal, (27) Paramillo, (28) Frontino-Urrao, (29) Citará, (30) Tatamá, (31) Duende, (32) Farallones de Cali, (33) Cerro Plateado, (34) Santa Marta, (35) Sonsón, and (36) Altiplano Cundiboyacense.

Climate Sensitivity in the Espeletia Complex

Climate sensitivity in the Espeletia complex was modeled by comparing current and future (2050 RCP 8.5) distributions under a niche conservatism scenario (i.e., range shift via migration). Occurrences of Espeletia samples for the 36 Colombian Páramo complexes were downloaded as georeferenced specimen data from providers served by GBIF4. From an initial set of 15,466 specimen records downloaded, a total of 1,432 (Supplementary Table S1), comprising 28 different Espeletia taxa according to Supplementary Table S2, were retained after filtering for georeferencing consistency (i.e., non-redundant records within the study area) and number of records. In order to avoid bias due to distributions with low number of input points (), specially for presumably rare taxa (), only putative species with at least 10 original records were retained (three taxa with more than 100 records, Supplementary Figure S1A). This threshold follows recommendation for narrow-ranged species, that is those with a moderate prevalence class (i.e., 27 taxa were present in less than seven different Páramo complexes, and only one taxon was present in 20 Páramos, Supplementary Figure S1B). This latter trend is expected in Espeletia species, for which distance (; ) and ecology () are the major speciation drivers.

Historical climate data was extracted from the WorldClim5 database, at a 30 s resolution, using the georeferencing of each record. All 19 bioclimatic variables were downloaded in the baseline period 1970–2000. The same variables projected to 2050, under the IPCC RCP 8.5 scenario, were gathered using CCAFS’s downscaled delta method (). This projection averages 33 different global climatic models at a 1 km2 resolution. A Digital Elevation Model (DEM) variable was further included. Prioritization of environmental variables () is crucial to enhance model complexity () for climate sensitivity assessments. By improving variable selection (), collinearity () and sampling bias () can be minimized. Collinearity in particular refers to a strong correlation between two or more predictor variables, and may cause instability in parameter estimation (). A measure to diagnose collinearity is the Variance Inflation Factor (VIF), the square of the multiple correlation coefficient resulting from regressing a predictor variable against all other predictor variables (). Hence, we performed variable selection by computing the VIF score using the usdm (Uncertainty Analysis for Species Distribution Models) package () in R v.3.3.1 (R Core Team), and drew Pairwise Pearson correlations (r) among all variables by means of the same software (Supplementary Figure S2). Five bioclimatic variables and DEM did not exhibit collinearity () at VIF < 10 (Supplementary Table S3), with absolute r-values below 0.7 from 0.05 to 0.65, and therefore were retained hereinafter.

The six non-collinear variables fed Maxent model v.3.4.1. (), a “machine learning” outline that uses maximum entropy and Bayesian inference to estimate occurrence probability distributions. Ten split-sample models were run, after checking for stabilization in the ROC curves (), with an average of 10 k pseudo-absences for each of the 28 putative taxa. Each replicated run type was set to “cross-validate” with 500 iterations (). Presence ≥ 95% was binary reclassified () and the average across repetitions was recorded in each case. This procedure was performed for the historical dataset (1970–2000). To predict distribution changes of Espeletia, we projected trained models on current presence localities and present background climate to 2050 under the RCP 8.5 scenario through Maxent’s “projection layer” tool. Current model validation considered AUC, an aggregate measure of performance across all possible classification thresholds. Minimum average testing AUC was above 0.98 (Supplementary Table S4). As AUC alone (Supplementary Table S5) may be deceiving (; ), the AIC and BIC criteria were also computed and checked for consistency (Supplementary Table S6).

All 56 distribution maps were converted to binary data, using a threshold of the 20th percentile, as suggested by the Maxent v.3.4.1 software (). These maps were then averaged to generate spatial information of current and future presence probabilities across all 28 taxa. Finally, the current distribution map was subtracted to the one obtained under a climate change scenario in order to estimate climate sensitivity in the Espeletia complex. Values close to zero in the sensitivity index meant few changes in the forecasted distribution of Espeletia due to climatic variation. On the other hand, negative and positive values corresponded to areas where the presence probabilities were, respectively, diminished or boosted by climate change.

Adaptive Capacity at the Study Area

We computed the adaptive capacity as the scaled additive effects of six main variables that expressed the capacity of the overall Espeletia biodiversity to respond without migration to the effects of climate change. Enhancing factors of the biodiversity’s adaptive potential were diversity, protected areas, and forest area around the Páramo. On the contrary, limiting factors of the plastic response were agriculture, mining, and population density.

Diversity may buffer the effects of climate change by enriching complex interactions, such as adaptive inter-specific introgression () and hybrid speciation (; ; ). Hybridization may allow for local adaptation in transitional environments and bi-directional transfer of adaptive variation through introgression (). Espeletia in particular exhibits weak species boundaries, rampant gene-flow and reticulate evolution (; ), leading to a high incidence of natural hybrids (; ). In order to account for the enhancing role of diversity on the adaptive potential, the number of Espeletia taxa was totalized in a 1 km grid using the occurrence data downloaded from GBIF (see footnote), as described in the previous section, and scaled from 0 to 1.

Another enhancing factor of the biodiversity’s adaptive potential is protected areas, where human intervention does not represent a direct threat, allowing biodiversity to be conserved while naturally coping with climate change. Protected areas selected for this analysis spanned 1,050 natural areas, according to the System of National Natural Parks (see footnote), and included 43 National Natural Parks, 663 Civil Society Reserves, 54 Regional Natural Parks, and 57 National Protected Forest Reserves. These categories were, respectively, rated as 1, 0.3, 0.1, and 0.1, following . Higher ratings were assigned to areas that presumably provide greater protection for biodiversity due to their low anthropogenic pressure and extent.

Forest patches around the Páramo complexes () may also contribute buffering climate change effects by facilitating () migration and gene-flow, which could increase the frequency of existing genetic variants adapted to particular conditions (). Additionally, forests provide conditions that may favor thermal and water regulation, decreasing the magnitude of climate change itself. Therefore, forest area in the buffer zone of the Páramo (2 km around the tree line) was gathered from the IDEAM (see footnote), and converted to binary points (1 = presence) in a 1 km grid.

On the other hand, major limiting factors of the plastic response are agriculture (), mining () and population density within the study area and in its buffer zone (). In order to account for these threats, agricultural and mining areas were downloaded from the Geographic Institute Agustín Codazzi (see footnote) for the period 2010–2012, and transformed to binary data (0 = presence) in a 1 km grid within the study area and in its buffer zone (2 km around the tree line). The population density indicator was computed as the inverse value of the normalized rural population density, scaled from 0 to 1, so that 1 meant absence of human settlements and 0 was the population density at the most populated rural zones in the study area. The population information was gathered from data projected to 2015 by the National Administrative Department of Statistics6.

Climate Vulnerability Index for Espeletia

In order to estimate climate vulnerability at the overall Espeletia complex, we simultaneously coupled its climate sensitivity with its adaptive capacity. Negative values of the sensitivity index meant areas where the presence probabilities of Espeletia were diminished due to climate change. Meanwhile, lower adaptive capacity scores spoke for a prominent role of limiting factors of the plastic response (i.e., agriculture, mining, and population density) as compared to enhancing drivers of the biodiversity’s adaptive potential (i.e., diversity, protected areas, and forest area around the Páramo). Therefore, the climate vulnerability score was computed as the complement value of the additive contributions of the sensitivity and the adaptive capacity indices. The climate vulnerability score was scaled from 0 to 1 for comparative purposes. Pairwise Pearson correlations (r) among all estimated scores were computed, and boxplots and ridgeline plots were drawn for the three indices in all 36 Páramo complexes. Computing and plotting were done in R v.3.3.1 (R Core Team).

Results

Climate Sensitivity of Espeletia Varied Widely Across Páramo Complexes

Climate sensitivity was estimated across all 36 Colombian Páramo complexes by comparing current (baseline period 1970–2000, Figure 2A) and future (2050 RCP 8.5, Figure 2B) average distribution probabilities of 28 different Espeletia taxa (nine had average sensitivity scores below zero). Values surrounding zero in the sensitivity index (−0.029 to 0.030, 42% of the study area, Figure 2C) stood for few changes in the forecasted distribution of Espeletia due to climate change, and were predominant in the north of the Eastern Cordillera, the Cruz Verde-Sumapaz Páramo complex south of the Bogotá montane savanna (the southwestern wetland complex of the larger Eastern Cordillera’s Andean highland plateau), and the small sky islands of the Western Cordillera.

Given a niche conservatism scenario, negative estimates of the climate sensitivity (≤−0.030, 32% of the study area, Figure 2C) predicted range losses based on the averaged presence probabilities. Páramo areas more likely to experience range losses were those in the Sierra Nevada de Santa Marta, and the Central Cordillera South of Las Hermosas (i.e., Nevado del Huila-Moras, and Guanacas-Puracé-Coconuco), as well as those in the northern corner of the Eastern Cordillera (i.e., Sierra Nevada del Cocuy, Pisba, and Tota-Bijagual-Mamapacha), where diverse summits surrounded the topographically intricate Chichamocha inter-Andean canyon. On the other hand, range shifts via migration enriched presence probabilities in other localities (positive estimates of the climate sensitivity ≥ 0.031, 27% of the study area, Figure 2C). Range gains were only predicted for Páramo areas in the Central Cordillera (Los Nevados, Chilí-Barragán, Las Hermosas, and La Cocha-Patascoy, Supplementary Figure S3) where the dominant taxon was E. hartwegiana, as well as in the summits of the northern tip of the Eastern Cordillera (Perijá).

Adaptive Capacity Was Jeopardized Across Most Páramo Complexes

The adaptive capacity across 36 Colombian Páramo complexes was computed using a scaled additive model of enhancing (Figures 3A–C) and limiting (Figures 3D–F) factors of the biodiversity’s adaptive potential and the plastic response. The adaptive capacity index (Figures 4B, 5A, Supplementary Figure S3A, and Supplementary Table S7) was higher (≥0.72, 34% of the study area) in species-rich Páramo areas (Pearson’s r = 0.19 ± 0.01, Figure 3A and Supplementary Figure S4), and matched extensive National Natural Parks (r = 0.59 ± 0.01, Figure 3B and Supplementary Figure S4), like the Eastern Cordillera’s Sierra Nevada del Cocuy, Chingaza and Sumapaz. The Páramo adaptive capacity index was also moderately superior (0.62–0.71, 29% of the study area, Figure 4B) in the most representative National Natural Parks in the Central Cordillera (Los Nevados and Nevado del Huila) as well as in the Sierra Nevada de Santa Marta. In turn, the western slopes of the Páramo complexes at the Eastern and Central Cordilleras exhibited lower adaptive capacity values (≤0.61, 37% of the study area), and were enclosed by less forest patches (30% on average, Figure 3C, r = 0.39 ± 0.01, Supplementary Figure S4).

FIGURE 3

). Forest patches (C) around the Páramo complexes (2 km around the tree line, 1 km grid) were recorded as presence points (1 = presence), while agriculture (D) and mining (E) for the period 2010–2012 were marked as null data points (0 = presence) within the study area and in its buffer zone (2 km around the tree line, 1 km grid). Páramo complexes numbers in (A) follow and : (1) Perijá, (2) Jurisdicciones-Santurbán, (3) Tamá, (4) Almorzadero, (5) Yariguíes, (6) Cocuy, (7) Pisba, (8) Tota-Bijagual-Mamapacha, (9) Guantiva-Rusia, (10) Iguaque-Merchán, (11) Guerrero, (12) Rabanal y río Bogotá, (13) Chingaza, (14) Cruz Verde-Sumapaz, (15) Los Picachos, (16) Miraflores, (17) Belmira, (18) Nevados, (19) Chilí-Barragán, (20) Las Hermosas, (21) Nevado del Huila-Moras, (22) Guanacas-Puracé-Coconuco, (23) Sotará, (24) Doña Juana-Chimayoy, (25) La Cocha-Patascoy, (26) Chiles-Cumbal, (27) Paramillo, (28) Frontino-Urrao, (29) Citará, (30) Tatamá, (31) Duende, (32) Farallones de Cali, (33) Cerro Plateado, (34) Santa Marta, (35) Sonsón, and (36) Altiplano Cundiboyacense.

FIGURE 4

and , that is: (1) Perijá, (2) Jurisdicciones-Santurbán, (3) Tamá, (4) Almorzadero, (5) Yariguíes, (6) Cocuy, (7) Pisba, (8) Tota-Bijagual-Mamapacha, (9) Guantiva-Rusia, (10) Iguaque-Merchán, (11) Guerrero, (12) Rabanal y río Bogotá, (13) Chingaza, (14) Cruz Verde-Sumapaz, (15) Los Picachos, (16) Miraflores, (17) Belmira, (18) Nevados, (19) Chilí-Barragán, (20) Las Hermosas, (21) Nevado del Huila-Moras, (22) Guanacas-Puracé-Coconuco, (23) Sotará, (24) Doña Juana-Chimayoy, (25) La Cocha-Patascoy, (26) Chiles-Cumbal, (27) Paramillo, (28) Frontino-Urrao, (29) Citará, (30) Tatamá, (31) Duende, (32) Farallones de Cali, (33) Cerro Plateado, (34) Santa Marta, (35) Sonsón, and (36) Altiplano Cundiboyacense.

FIGURE 5

However, positive effects due to protected areas and forest coverage on the adaptive potential were offset by limiting factors of the plastic response, primarily agricultural (6% of the study area and its buffer zone, Figure 3D, r = −0.47 ± 0.01, Supplementary Figure S4) and mining (7% of the study area and its buffer zone, Figure 3E, r = −0.41 ± 0.01, Supplementary Figure S4) pressures, mainly in the western slopes of the Páramo areas at the Eastern and Central Cordilleras. This let to the lowest adaptive potential estimates (≤0.52, 20% of the study area, Figure 4B) at the western slopes of the Páramo complexes in the Cauca province, in the Central Cordillera, and the Boyacá and Santander provinces (i.e., Almorzadero and Altiplano Cundiboyacense Páramo complexes), in the Eastern Cordillera. Rural population (Figure 3F) hampered the adaptive potential in all three Cordilleras and Páramo areas (r = −0.28 ± 0.01, Supplementary Figure S4).

Diverse Páramos in the Eastern Cordillera Were More Vulnerable

The scaled complement of the additive contributions of the sensitivity (Figures 2C, 5B and Supplementary Figure S3B) and the adaptive capacity (Figure 4B) indices were used to compute the climate vulnerability score (Figures 4C, 5C and Supplementary Figure S3C). The most vulnerable (≥0.57, 11% of the study area, Figure 4C) Páramo areas coincided with diverse Páramos in the Eastern Cordillera (Figure 3A) because of more distribution losses (r = 0.39 ± 0.01, Supplementary Figure S4). Agriculture (r = 0.48 ± 0.01, Supplementary Figure S4), mining (r = 0.36 ± 0.01, Supplementary Figure S4) and rural population density (r = 0.23 ± 0.01, Supplementary Figure S4) also boosted their vulnerability, mainly on the western slopes in the middle of the Central Cordillera (i.e., Guanacas-Puracé-Coconuco), and the northeast corner of the Eastern Cordillera (i.e., Almorzadero, Pisba, and Altiplano Cundiboyacense). On the other hand, the less vulnerable (≤0.38, 57% of the study area, Figure 4C) Páramo areas were those with the highest adaptive capacity (≥ 0.72, Figure 4B, r = −0.93 ± 0.01, Supplementary Figure S4) due to large National Natural Parks (Sumapaz, Chingaza, Sierra Nevada del Cocuy, Figure 3B), as well as those for which predicted range shifts were neutral (−0.029 to 0.030, Figure 4A, i.e., the Cruz Verde-Sumapaz Páramo complexes in the Eastern Cordillera) or positive (≥0.031, Figure 5B, Los Nevados, Chilí-Barragán, and Las Hermosas in Central Cordillera).

Discussion

The climate sensitivity score of Espeletia to climate change indicated range losses (negative values) widespread across most Páramo complexes. This suggests limited persistence of the Espeletia spp. populations via migration in the majority of Páramos, with some exceptions in the Central Cordillera possibly due to its higher connectivity. An alternative for populations to persist despite changing environments would be adaptation s.l. However, positive effects on the adaptive potential, like those due to protected areas and forests, are likely to be offset by limiting factors of the plastic response, such as rising agriculture and mining in the Páramo areas and their surroundings. Overall, this suggests that diverse Páramo areas of the Eastern Cordillera are the most vulnerable to changing climate. Climate change pushing Espeletia toward mountain summits where there is nowhere else to go (), together with a general inability of those populations to adapt, might constrain the rapid diversification that has made the Espeletia complex so unique.

Range Losses in the Espeletia Complex Will Widespread Throughout Most Páramos

Widespread range losses in the distribution of Espeletia are expected by 2050, a pessimistic prediction in line with previous modeling efforts (; ; ; ; ; ). This tendency is broadly defined as thermophilization, a progressive decline of cold mountain habitats and their biota (). Large reductions in modeled area and important upward shifts in the distribution of 28 Espeletiinae are already predicted by 2070 at the Venezuelan Páramo (). Despite our projection is carried out for a shorter time frame at another, wider, study area (Colombian Andes spanning the Eastern, Central and Western Cordilleras besides Sierra Nevada de Santa Marta), overall results point toward the same trend and are complementary. Similar predictions of climate sensitivity at the Páramo complexes of the Ecuadorian Cordilleras (1 sp.) were not considered as part of this study because of the difficulties in compiling and comparing data from different national repositories, specifically for some of the variables that compose the adaptive capacity score (i.e., agriculture and mining pressures). Yet, considering analogous approaches in the Ecuadorian Páramo (; ; ) is essential to gather a more comprehensive understating on the consequences of climate change at high Andean ecosystems north of the Huancabamba depression, a major biogeographical barrier for high-elevation plant taxa ().

There are two important exceptions for the widespread range losses in the distribution of Espeletia in the Colombian Andes. First, small and fragmented Páramo sky islands in the Western Cordillera exhibit neutral climate sensitivity, a counterintuitive projection that also pointed out for Espeletiinae in isolated Páramo areas at Cordillera de Mérida. This reiterative finding suggests that scarce topographical complexity (terrain ruggedness) in recently colonized small islands () may have favored habitat generalists, while intricate local-scale environmental heterogeneity at larger Páramo complexes could have triggered an explosive radiation (; ) of habitat specialists, each with a narrow ecological niche for migration to occur (). Modern colonization and limited topographical complexity, together with higher connectivity, may also account for the second exception, which is that range gains are only predicted for Páramo areas in the Central Cordillera (Los Nevados, Chilí-Barragán, Las Hermosas, and La Cocha-Patascoy), where the dominant taxon is E. hartwegiana. In line with this, Páramo areas more likely to exhibit range losses are those in the topographically complex Eastern Cordillera’s northeast corner, where diverse (; ; ) summits surround the deep (1,100 m) Chichamocha inter-Andean canyon.

Range shifts presuppose niche conservatism, but also depend on the intrinsic dispersal ability of the taxa (; ). In organisms with limited seed dispersal such as Espeletia, other ways for migration are key (). For instance, migration via sporadic pollen flow would not be as constrained as seed dispersal (). Alternately, hybrids with intermediate niche requirements () may serve as bridges for migration in a plant complex with rampant natural introgression (; ; ). Finally, sporadic episodes of habitat connectivity () in the lifespan of a long-lived organism such as Espeletia () could provide a window for migration across modern barriers.

Key traits may relax the niche conservatism assumption by conferring the ability to colonize new habitats. For instance, underground plant-soil interactions (; ; ) across different water levels () and cellular-anatomical physiological adaptations () may serve as pre-adaptations for a wider range expansion (). Plastic traits may also broad the spectrum of novel climates (). Coupling this modeling with eco-physiological studies in Espeletia (; ) is hence a key research line to better track future range shifts in the Páramo.

Climate Change May Constrain the Rapid Diversification of the Espeletia Complex

The fact that diverse areas in the Eastern Cordillera are vulnerable to climate change is not only counterintuitive, but it also means that the long-term rapid diversification rate, striking in the Espeletia complex, could be jeopardized. Diversity is regarded as a buffer of the effects of climate change by enriching interactions (), such as facilitation (; ; ; ; ), adaptive inter-specific introgression (; ) and hybrid speciation (), of which there is abundant evidence in Espeletia (; ; ; ). Yet, these positive effects are offset by distribution losses in sensitive Páramo localities, besides rising mining, agriculture and population density. These are Páramo’s current major anthropogenic threats (; ) primarily in the western slopes of the Eastern and Central Cordilleras.

Permeable species barriers due to pervasive hybridization make the species concept highly debatable within the Espeletia complex. Therefore we refrained from emphasizing single-species responses that are likely to be overwhelmed by meta-population dynamics. suggested that the entire subtribe Espeletiinae requires an in deep taxonomic revision () because most of Cuatrecasas’ (2013) species are thought to be paraphyletic. This trend is not unexpected since the available taxonomy greatly relies on potentially plastic morphological and reproductive traits () that are susceptible to be driven by environmental effects, besides the fact that it disregards rampant hybridization (). Then, a more appropriate scenario for the young rapidly-evolving Espeletia complex with weak species boundaries would be that ecotypes in the same locality, despite depicting microhabitat divergence, are capable of sharing through gene flow and introgression adaptive genetic variants (). Gene swapping because of porosities in the species bounds would imply that (1) the effective size of “rare” taxa with few records is higher in terms of standing adaptive variation, and (2) concerted climate responses is feasible across lineages within the same Páramo sky island. Inter-mountain exchange is also plausible, as reported in other tropical sky islands (; ).

A possible way forward for populations to persist is the very same driver that is thought to contribute with the unbeatable diversification rate in the Páramo flora () – local-scale adaptive variation (; ). As part of this study we considered major enhancers and constrainers of the adaptive potential at a regional scale. However, populations may still harbor intrinsic genetic variation naturally selected to cope with environmental differences at the microhabitat level (). In highly heterogeneous mountainous environments, local-scale environmental effects are known for shaping genetic diversity and generating morphological diversity (), which may prove useful in the reaction to climate change within the same locality. There is already compiling evidence of ecological-driven divergence in the radiation of Andean Espeletia (), besides allopatric differentiation and isolation by distance (; , ). Ecological parapatric speciation is a likely consequence of local patterns of environmental variation () typically found at the Páramo ecosystem (). For instance, frost is more extreme in the wind-exposed high valley slopes than in the wind-sheltered depressions of those valleys or in the upper boundary of the cloud forest. Soil moisture content, higher in the well-irrigated high valley depressions and in the upper boundary of the cloud forest than in the drier slopes of the valleys, could also contribute to this divergence. It now remains to be explored whether this mosaic of environmental heterogeneity has enforced enough heritable polygenic () phenotypic variation as a pre-adaptation () to the new selective forces imposed by climate change.

One potential caveat of our results precisely concerns data availability at the micro-scale level. Bioclimatic data is reliably modeled throughout the Páramo range at a 1 km2 spatial resolution, which unfortunately overlooks micro-environmental drivers of the local-scale genetic adaptation (; ). This is of particular importance at mountain/alpine ecosystems (; ), where microhabitats may serve as refugia (), like across the treeline due to active landform processes (; ; ). Populations may rely on specific small-scale habitat attributes that are likely to be overlooked by SDMs (), such as topographic, geomorphic, and edaphic features, as well as the distribution of other taxa – e.g., competitors and facilitators (). Due to this, SDMs may exhibit substantial uncertainty () and lack of validation () at local scales. Therefore, niche models could be insensitive to micro-environmental processes, leading to an underestimation of the mechanisms governing populations’ responses. Even if bioclimatic models aim representing the realized niche, they may then disregard the fundamental niche (). Despite statistical imputation-like approaches such as environmental downscaling have been considered to model climate at a higher resolution for specific Páramo complexes (), consistently applying such techniques across cordilleras seems unfeasible because accuracy would be constrained by the accessibility of fine-scale bioclimatic data from regional weather stations. These restraints are likely to be overcome soon by incorporating microclimate into SDMs () through mechanistic algorithms () physiological models (; ), hydrological equations (; ), in situ data logging, and remote sensing (; ).

Another potential caveat of our inferences relate with data availability across the entire Páramo area in the northern Andes. Some components of the adaptive capacity score (i.e., protected areas, agriculture, and mining pressures) are downloadable only through national repositories, making standardizations and comparisons across all Páramo complexes in different countries impractical. Yet, it is worth to clarify that climate sensitivity, adaptation capacity and vulnerability indices are not meant to be absolute scores but rather relative values to compare among Páramo complexes. In this sense, the power of our inferences to rank Páramo areas is unlikely limited by the pessimistic predictions of the future climate scenario RCP 8.5, making our vulnerability index robust. Furthermore, our whole integrative analytical pipeline is in line with the most up-to-date approaches to carry out climate change vulnerability assessment (; ; ).

By assessing climate sensitivity and adaptive capacity across Páramo sky islands, we predicted climate change vulnerability in the rapidly evolving Espeletia complex, a syngameon () likely shaped by incomplete ecological speciation () and adaptive introgression (; ). However, in order to gather a more comprehensive view of climate change responses in the high mountain ecosystems of the northern Andes, other diverse and highly abundant plant genera in the Páramo () should be considered under a similar analytical framework as the described here. Immediate candidates to follow up this approach are Bartsia (), Chusquea (), Diplostephium (), Hypericum (), Loricaria (), Lupinus (; ; ), Oreobolus (; ), Puya (), and Senecio (; ). These combined efforts would ultimately reveal whether the fastest evolving biodiversity hotspot on earth, and in general tropical high mountain ecosystems (; ; ), have a chance to persist under current environmental and anthropogenic threats.

Perspectives

Besides range shifts (; ; ; ; ) and adaptive responses (; ; ; ), introgression, hybridization (), and polyploidization (; ; ) which have not been explicitly explored in this study, may also provide innovative () adaptive sources in plants (; ). The Espeletia complex is known for a high incidence of natural hybrids as told by phenotypic and genetic markers (; ). Hence, assessing the role of reticulate evolution and adaptive introgression for the Espeletia complex to face a changing environment will require sampling of hybrids and their putative parental populations across the Páramo’s current micro-ecological mosaic.

Even though hybrids could occupy intermediate niches, persistence of populations in alpine environments that face climate change is regarded as mostly mediated by local-scale (i.e., microhabitat) variation (; ). For example, environmental heterogeneity may provide new suitable locations for migrants within only a few meters of their current locations (), driving in this way plant responses to warming (). More localized environmental data must then be gathered (), since current climate repositories lack the resolution needed to describe the microhabitat heterogeneity. Topographic and soil properties, often overlooked by climatic models, are presumably key in defining Páramo’s wet and drier microhabitats in the wind-sheltered high valley depressions, and in the more wind-exposed slopes. This local-scale environmental trend has diversity implications for other endemic plant species in the northern Andes (; ). Besides high-resolution environmental data, eco-physiological (; ) and ecological adaptive traits (e.g., ; ) must be noted more systematically at local scales in natural surveys and field tests, like reciprocal transplant assays of ecotypes among microhabitats (as in ), and space-by-time substitution trials across altitudes () and habitats (, ). Ultimately, micro-scales keep genetic variance (, ; ; ; , ; ; ), morphological diversity, and adaptive trait variation (; ).

Conservation efforts may also buffer the impact of climate change on Páramo ecosystems by enhancing their adaptive capacity. Vast areas of Páramo land now safe from conflict are not being rapidly protected, making them susceptible to deforestation for cattle and agriculture (), and illegal logging () and mining (). Thus, improving the current network of protected areas would help minimizing ongoing threats (), especially in the highly diverse Páramo complexes of the Eastern Cordillera where putatively new species are still being discovered (; ). Alongside adoption of mitigation policies, adaptation policies () may help safeguarding Páramos. Finding more sustainable land uses by empowering local communities and developing ecotourism are then as essential to relieve human impact on tropical-alpine plant diversity in the northern Andes.

Last but not least, extending the modeling approach implemented in this work to other key pressures (e.g., fire, ; ; ) and plant groups () in the Páramo, as well as to other sky islands around the mountains of the world (; ; ; ), and more broadly to other island-like systems (; ; ; ), will help understanding climate change effects on unrelated taxa experiencing similar evolutionary processes (; ; ). Such systems offer natural experiments to assess the role of colonization and adaptation (; ; ) in the past and ongoing responses to climate change, which undeniably will complement ecological predictive modeling.

Statements

Data availability statement

All datasets generated for this study are included in the article/Supplementary Material and in the Dryad repository under https://doi.org/10.5061/dryad.gf1vhhmmp. Further inquiries can be directed to the corresponding author.

Author contributions

AC, JV, and SM conceived the study. JV compiled the data. JV and JM analyzed the data with contributions from AC and JL. AC, JV, and SM interpreted the data. AC wrote the manuscript with contributions from JV, and further edits made by the other co-authors. All authors contributed to the article and approved the submitted version.

Funding

Grants 4.1-2016-00418 from Vetenskapsrådet (VR) and BS20170036 from Kungliga Vetenskapsakademien (KVA) are thanked for funding AC. The University of California (CA, United States) Conservation Grant, the Wildlife Conservation Society (WCS, NY, United States) and Facultad de Ingeniería y Administración from Universidad Nacional de Colombia (Palmira, Colombia) are also acknowledged for supporting JV’s presentation at ESRI User Conference during 2018 at San Diego (CA, United States). The editorial fund from the Colombian Corporation for Agricultural Research (AGROSAVIA) is thanked for financing this publication.

Acknowledgments

We thank early discussions on Espeletia ecological divergence with A. Antonelli, C.D. Bacon, I. Calixto-Botía, N. Contreras-Ortiz, I. Jiménez, J. Mavarez, and M. Pineda. J. Ramirez-Villegas, and C.C. Sosa assistance in the prioritization of bioclimatic variables is recognized. V. Bermúdez, L. Gaona, and A. Salamanca are acknowledged for administrative support. M.J. Torres-Urrego aid while drafting this work is deeply appreciated, too. Some of the analyses and ideas presented in this manuscript were polished thanks to comments from participants of the Facultad de Ingeniería y Administración (UNAL) FiaInnova Workshop, and the International Center for Tropical Agriculture (CIAT)’s Decision & Policy Analysis (DAPA) Annual Meeting, respectively, held in November 2018 and April 2019 at Palmira (Colombia). Topic editors and reviewers are also thanked for envisioning a well-timed Research Topic on the “Consequences of Climate Change for Plant Biodiversity in High Mountain Ecosystems,” and for improving an early version of this manuscript through their thoughtful recommendations.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2020.565708/full#supplementary-material

Supplementary Figure 1

Frequency distributions of the number of records per taxon (A), and the number of Páramo complexes (B) spanned by each of the 28 Espeletia taxa across Páramo sky islands in the Colombian Andes. Number of records, and number of Páramo complexes spanned for each taxon are drawn from Supplementary Table S2.

Supplementary Figure 2

Pairwise Pearson correlation coefficients (r) among 19 WorldClim bioclimatic variables (in the baseline period 1970–2000), and Digital Elevation Model (DEM) at a 30 s (1 km2) spatial resolution. Five bioclimatic variables and DEM did not exhibit collinearity problems (Supplementary Table S3) and were retained.

Supplementary Figure 3

Ridgeline plots for the adaptive capacity (A), climate sensitivity (B), and vulnerability (C) scores across 36 Páramo complexes in the Colombian Andes. For comparative purposes, Páramo complexes are sorted in (A) from those with a low adaptive potential to anthropic pressures according to the enhancing and limiting factors depicted in Figure 3, in (B) from those likely to experience a range loss (distributions below zero) to those predicted to exhibit a range gain (distributions above zero) by 2050 (under RCP 8.5), while in (C) Páramo complexes at the top are the most vulnerable. Ridgeline plots summarize grid estimates shown in Figure 4.

Supplementary Figure 4

Pairwise Pearson correlation coefficients (r) among climate sensitivity (S), adaptive capacity (AC), climate vulnerability (V) and the components of the adaptive capacity in the Espeletia complex across Páramo sky islands in the Colombian Andes. The adaptive capacity enhancing components were diversity (D), protected areas (PA) and forest area around the Páramo (F), while its limiting factors were agriculture (A), mining (M) and rural population density (RPD). The area of each of the 36 Páramo complexes is also considered in hectares (Ha). Correlation estimates and 95% confidence intervals are presented above the diagonal, while below the diagonal circles are colored and sized accordingly. Estimates are gathered from Supplementary Table S7 and are negatively transformed for agriculture (A), mining (M) and rural population density (RPD) in order to aid interpretation.

Supplementary Table 1

Occurrence data of 28 Espeletia taxa across Páramo sky islands in the Colombian Andes. Occurrences were gathered as georeferenced specimen data from GBIF. Given a starting set of 15,466 specimen records, a total of 1,432 were kept after cleaning for georeferencing consistency (i.e., non-redundant records located within the study area) and number of records per species (to avoid distributions with less than 10 points, as suggested by for narrow-ranged species). Espeletia taxonomy follows .

Supplementary Table 2

Number of records for each of the 28 Espeletia taxa across Páramo sky islands in the Colombian Andes, and altitudinal range for each taxon. Number of records, number of Páramo complexes spanned, and minimum and maximum altitudes for each taxon are gathered from Supplementary Table S1.

Supplementary Table 3

Variance Inflation Factor (VIF) for five non-collinear (VIF < 10) WorldClim bioclimatic variables considered for modeling at a 30 s (1 km2) spatial resolution for the current and future scenarios, besides a Digital Elevation Model (DEM). Absolute r-values ranged from 0.05 (BIO18∼DEM) to 0.65 (BIO15∼BIO18).

Supplementary Table 4

Current model validation based on average Area Under the ROC Curve (AUC) and porcentual contribution of the six non-collinear environmental variables for each of the 28 Espeletia taxa. AUC in the training and testing datasets are an aggregate measure of performance across all possible classification thresholds. Pseudo-absences (number of background points) are also averaged for each of the 28 putative taxa. Un-averaged statistics are kept in Supplementary Table S5. The six non-collinear variables included five WorldClim bioclimatic variables at a 30 s (1 m2) spatial resolution. Environmental variables (Supplementary Table S3) are abbreviated as follows: (BIO2) Mean Diurnal Range, (BIO4) Temperature Seasonality, (BIO15) Precipitation Seasonality, (BIO18) Precipitation of Warmest Quarter, (BIO19) Precipitation of Coldest Quarter, and (DEM) Digital Elevation Model.

Supplementary Table 5

Current model validation based on Area Under the ROC Curve (AUC) and relative contribution of the six non-collinear environmental variables for each of the 10 runs in the 28 Espeletia taxa. The following overall summary statistics are kept: optimum training subset, regularized training gains, un-regularized training gains, effective iterations, training AUC, optimum testing subset, test gains, test AUC, AUC standard deviation, pseudo-absences (number of background points), and entropy. The following relative contribution scores are kept for the six environmental variables: contribution, permutation importance, training gain, test gain and AUC without and with only each variable. The six non-collinear environmental variables (Supplementary Table S3) included five WorldClim BIO variables at a 30 s (1 km2) spatial resolution, and are coded as: (BIO2) Mean Diurnal Range, (BIO4) Temperature Seasonality, (BIO15) Precipitation Seasonality, (BIO18) Precipitation of Warmest Quarter, (BIO19) Precipitation of Coldest Quarter, and (DEM) Digital Elevation Model.

Supplementary Table 6

Current model validation based on Akaike (AIC) and Bayesian (BIC) Information Criteria per repetition of each of the 28 Espeletia taxa. Optimum training and testing subsets are summed under NTR+TS.

Supplementary Table 7

Climate sensitivity (S), adaptive capacity (AC), and climate vulnerability (V) in the Espeletia complex across Páramo sky islands in the Colombian Andes. The adaptive capacity enhancing components were diversity (D), protected areas (PA) and forest area around the Páramo (F), while its limiting factors were agriculture (A), mining (M) and rural population density (RPD). Diversity (D), which is the number of Espeletia taxa totalized in a 1 km grid based on the occurrence data downloaded from GBIF, and the inverse value of the normalized rural population density (RPD) were scaled from 0 to 1. A total of 1,050 protected areas (PA) were rated as 1, 0.3, 0.1, and 0.1 for National Natural Parks, Civil Society Reserves, Regional Natural Parks and National Protected Forest Reserves, respectively, following (). Forest patches (F) around the Páramo complexes (2 km around the tree line, 1 km grid) were recorded as presence points (1 = presence), while agriculture (A) and mining (M) for the period 2010–2012 were marked as null data points (0 = presence) within the study area and in its buffer zone (2 km around the tree line, 1 km grid). The area of each of the 36 Páramo complexes is also presented in hectares (Ha), as well as its name and ID number following Figure 2A, and .

References

Summary

Keywords

Neotropical alpine region, biodiversity hotspots, species distribution modeling, niche conservatism, range shifts, adaptation, Espeletiinae Cuatrec. (Asteraceae)

Citation

Valencia JB, Mesa J, León JG, Madriñán S and Cortés AJ (2020) Climate Vulnerability Assessment of the Espeletia Complex on Páramo Sky Islands in the Northern Andes. Front. Ecol. Evol. 8:565708. doi: 10.3389/fevo.2020.565708

Received

26 May 2020

Accepted

26 August 2020

Published

24 September 2020

Volume

8 - 2020

Edited by

Domenico Gargano, University of Calabria, Italy

Reviewed by

Rodolfo Gentili, University of Milano-Bicocca, Italy; Mauro Fois, University of Cagliari, Italy

Updates

Copyright

*Correspondence: Andrés J. Cortés,

These authors share senior authorship

This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics