The Thermal Tolerances, Distributions, and Performances of Tropical Montane Tree Species

Introduction

With anthropogenic climate change driving rapid increases in global temperatures, many plant species will face greater risks of thermal stress (Saxe et al., 2001; Parmesan and Hanley, 2015). If severe enough, thermal stress can limit the performance of individuals and populations, eventually leading to the retraction of some species' ranges (or shifts in species' ranges if local extinctions are offset by simultaneous invasions into cooler areas) and local extinctions (Chen et al., 2011; Feeley, 2012; Lenoir and Svenning, 2015). These local extinctions can in turn lead to decreased diversity through biotic attrition (Colwell et al., 2008; Wiens, 2016), changes in community composition (Feeley et al., 2011, 2013; Duque et al., 2015; Fadrique et al., 2018) and potential changes in important ecosystem services such as carbon sequestration/storage (Clark et al., 2003; Brienen et al., 2015), regional climate regulation (Cox et al., 2000; Luo, 2007), and food production (Tito et al., 2018).

In plants, thermal stress can be caused by the increase in respiration and/or the decrease in photosynthesis due to decreased efficiency or damage to leaf metabolic processes. At high temperatures, photosystem II (PSII) can be permanently damaged leading to decreased electron transport rates, photosynthetic failure, and eventual leaf death (Baker, 2008). The temperature (generally 40-60°C) at which irrecoverable damage occurs to PSII, can be referred to as a plant's maximum physiological thermal tolerance. The heat tolerances of plants can vary markedly between species and previous studies have found that they are generally correlated with large-scale latitudinal gradients in temperature, such that plants growing in hot equatorial habitats (e.g., rainforests) have higher heat tolerances on average compared to plants from colder, higher latitude habitats (e.g., temperate or boreal forests) (O'sullivan et al., 2017; Zhu et al., 2018). However, the change in thermal tolerances across latitude is not as steep as the corresponding change in air temperature. This suggests that many tropical species may have reduced thermal safety margins (difference between temperatures and tolerances) compared to temperate species and thus will be at elevated risk of thermal stress due to global warming (Perez et al., 2016; O'sullivan et al., 2017).

A similar relationship between plant thermal tolerances and temperature should also exist along smaller-scale environmental gradients, for example across elevational gradients. In other words, plant species at lower, hotter elevations should have higher tolerances than species at higher, colder elevations (e.g., O'sullivan et al., 2013). However, this hypothesis remains almost entirely untested (but note that measurements of tree species in a single highland plot in Peru do show lower thermal tolerances than species in lowland forests at similar latitudes; O'sullivan et al., 2017). Consequently, we do not know if tropical plant species in general are at elevated risk of thermal stress or if highland tropical species have relatively greater thermal safety margins that can help protect them against the negative effects of rising temperatures. Likewise, it is unknown how leaf thermal tolerances relate to plant species' performances across large-scale spatial and temporal thermal gradients. For example, are the ranges of species with lower thermal tolerances spatially restricted to areas with lower regional temperatures? And, within specific sites, are species with higher thermal tolerances experiencing performance advantages, and hence increasing in their relative abundances compared to less-tolerant species, as temperatures rise under global climate change?

To answer these questions, we assess the thermal tolerances of >550 individuals of 164 tropical canopy tree species growing in 10 forest inventory plots, each 1-hecatare in area (= 189 species × site populations). Study plots are distributed across a steep elevation gradient ascending from near sea level to near treeline (~3,000 m asl) in the northern Andes Mountains of Colombia (Duque et al., 2015; Peña et al., 2018; Agudelo et al., 2019). We analyze changes of plant thermal tolerances between elevations at the species and community level (intraspecific patterns were not analyzed but will be investigated in future studies). We conduct two sets of analyses to investigate the relationship between the thermal tolerances of species and their ability to tolerate high temperatures. In the first, we test the relationship between species' thermal tolerances and their observed geographic distributions in relation to regional-scale mean annual and maximum temperatures. In the second set of analyses, we test the relationship between species' thermal tolerances and their changes in abundance through time as recorded in each of the study plots. These analyses are used to test the hypotheses that (1) species with higher thermal tolerances will have geographic distributions including hotter areas than species with lower thermal tolerances, and that (2) due to rising temperatures, species with higher thermal tolerances will experience performance advantages over co-occurring species with lower thermal tolerances and thus will increase in relative abundance over time. Together these analyses provide valuable insight into the importance of temperature and thermal tolerances in setting current and future species distributions.

Methods

Study Area

This study was conducted using data collected from 10 permanent forest inventory plots in the northwest Andes Mountains of Colombia, located between 5°50′ and 8°61′ N Latitude and 74°61′ and 77°33′ W Longitude. Plot elevations range from 192 to 2,880 m asl. The average total annual rainfall recorded in the plots varies from ~ 1,750 to 3,500 mm year⁻¹, and mean annual temperature (MAT, °C) varies from 10 to 27°C (Table 1; Figure 1) with an adiabatic lapse rate of ~5.7°C km⁻¹ elevation. Mean Annual Temperatures in this area have been increasing steadily at the rate of approximately 1°C per century since the early 1900s (Rohde et al., 2013).

Table 1. Location and climate of study plots.

Figure 1. Map of study area showing the locations of the 10 study plots (1-ha each) in the Antioquia department of Colombia. Plot elevations are Puerto Triunfo: 192 m; Támesis: 550 m; Maceo: 1,024 m; Porce: 1,033 m; Anorí: 1,794 m; Angelópolis: 2,118 m; Ventanas: 2,120 m; Jardín: 2,550 m; Caicedo: 2,635 m; Belmira: 2 800 m.

Plot Data

The study plots are all 1-ha in area and were initially established between 2006 and 2009. In the initial plot censuses, all trees, palms and tree ferns with stem diameter at breast height (DBH) ≥ 10 cm were tagged, mapped, measured, and identified. In all cases, the point of DBH measurement was painted on the stem to ensure that re-measurements in subsequent censuses were at the same stem location. The study plots were all recensused multiple times. During each plot recensuses, DBH growth, individual recruitment, and mortality were recorded. In cases where the recorded diameter growth was either < -0.1 cm year⁻¹ or larger than 7.5 cm year⁻¹, the DBH of the second census was adjusted to avoid growth values out of this range (Condit, 1998; Condit et al., 2004). The Angelópolis, Anorí and Belmira plots were each recensused three times (4 total censuses including initial) with the last census completed in 2017. The Caicedo, Jardín, Maceo, Porce, Puertro Triunfo, Támesis, and Ventanas plots were each recensused twice. In all censuses, all plants were identified to the lowest possible taxonomic level using voucher samples deposited at the University of Antioquia's Herbarium (HUA) in Medellin, Colombia.

Previous studies have indicated that these forests are undergoing a process of compositional change called “thermophilization,” possibly due to rising temperatures (Duque et al., 2015). In thermophilization, species with distributions centered on lower, hotter elevations (i.e., more-thermophilic species) increase in relative abundance and species with distributions centered on higher, cooler elevations (i.e., less-thermophilic or cryophilic species) decrease in relative abundance through time. This pattern mirrors that observed at other neotropical montane sites (Feeley et al., 2011, 2013; Fadrique et al., 2018).

Study Species

Thermal tolerances were estimated for a total of 15 to 22 species in each of the 10 study plots (Table S1). Species were chosen to include the 10 most-abundant canopy tree species in each plot as well as an additional ~10 species spanning a wide range of geographic-based thermal optima and observed changes in abundance through time (see below). To estimate thermal tolerances, we collected sun-exposed leaves from three mature individuals of each species. Leaf samples were obtained using specialized pruning equipment and/or professional tree climbers. Once collected, the leaf samples were placed immediately in plastic “ziplock” bags with moistened paper towels and transported in coolers to laboratories at the National University of Colombia in Medellin for testing. In total, we sampled 164 species. The majority of species were sampled from just one plot each; 21 species were sampled from two plots each and two species were each sampled from three plots (Table S1). Hereafter, all species by plot combinations are treated as “species.”

Estimating Photosynthetic Thermal Tolerances

We estimated the maximum thermal tolerance (T₅₀) of each of the study species to be the temperature at which photosystem II suffers ≥50% irrecoverable damage as determined through measurements of leaf fluorescence (F_V/F_M) in leaf samples exposed to different temperature treatments. Methods for estimating T₅₀ were based on standard protocols modified from Krause et al. (2010, 2013). Prior to tests, we first measured the initial status of each leaf based on their pre-treatment F_V/F_M. Leaves were dark adapted for ~20 min. After dark adaptation, we measured initial fluorescence emission (F₀) and maximum total fluorescence (F_M) using an OS30p⁺ handheld fluorometer (Opti-Science, Hudson, NH USA). F_V/F_M was then calculated as the ratio of maximum variable (F_V = F_M-F₀) to maximum total fluorescence. Subsequently, we cut 2 cm disks from each leaf (avoiding major veins). Leaf disks were placed in Miracloth fabric pouches (abaxial and adaxial leaf surfaces covered by three and one layer of fabric, respectively) to prevent anaerobiosis (Krause et al., 2010) and placed inside waterproof plastic bags. Air was removed from the bags and bags were completely submerged for 15 min in preheated circulating water baths maintained at set temperatures of 22.0, 40.0, 42.2, 44.5, 47.0, 49.3, 51.5, 53.2, and 56.5°C. During heat treatments, leaves were kept under dim light and thus were not likely exposed to adequate quantities of light needed to induce the production of violaxanthin and zeaxanthin, which can ameliorate damage at high temperatures and possibly lead to lower estimates of heat tolerances (Krause et al., 2015). After 15 min, the leaf disks were removed from the water baths and bags, and placed in petri dishes lined with moist paper towels where they were left under dim light (~1 μmol photons m⁻² s⁻¹) at room temperature (20–23°C) and allowed to recover for ~24 h. After this recovery period, the leaf disks were dark-adapted for 20 min and F_V/F_M was measured.

To estimate the T₅₀ for each species, we modeled the relationship of F_V/F_M vs. treatment temperature for each species using a logistic non-linear least squares model with the “nls” function in the R stats package [i.e., nls(F_V/F_M ~ θ₁/(1 + exp(–(θ₂ + θ₃ × Temperature)))] where θ₁ is the control value of F_V/F_M (≈ 0.8) and θ₂ and θ₃ are the intercept and slope coefficients of the logit(Fv/Fm) ~ Temperature relationship, respectively; R Core Team, 2018). T₅₀ was then calculated as the temperature modeled to cause a 50% reduction in F_V/F_M compared to the control (Figure S1). To generate bootstrapped mean and 95% confidence level estimates of T₅₀ for each species, we reiterated the model of F_V/F_M vs. temperature 100 times. During each iteration, F_V/F_M and temperature data were sampled randomly with replacement to recalculate T₅₀. For some species (n = 18), estimates of T₅₀ were deemed unreliable due to unreasonable values (i.e., T₅₀ >55 or <40^oC) or very large 95% confidence intervals (i.e., 95% CI > 5°C; Table S1, Figure S1). These species were excluded from any subsequent analyses that used T₅₀.

Species Abundances and Changes in Abundance

We calculated the relative abundance (BA_rel) of each study species in each plot census as the total basal area of the focal species (= summed cross sectional area of all conspecific stems in a plot at breast height [1.33 m above ground]) relative to the total basal area for all stems in the plot (including stems of species not assessed for T₅₀; Table S2). We then estimated the annualized rate of change in relative abundance through time for each species as the slope of the linear regression of BA_rel vs. census date (ΔBA_rel, % year⁻¹; Table S1). Using the ΔBA_rel estimates, we also classified each species as either increasing or decreasing in relative abundance through time.

Geographic-Based Thermal Niches

In order to assess how the study species' T₅₀ relate to the species' large-scale geographic distributions, we used natural history observation records to estimate the range of temperatures over which each study species is known to occur. Specifically, for each of the study species, we downloaded all available collection and observation records using the BIEN data portal (queried via R package BIEN in March 2019). The BIEN data portal provides access to a number of different public data sources including the Global Biodiversity Information Facility (GBIF), SpeciesLink and various inventory plot networks, and records have undergone a base level of cleaning to minimize taxonomic and georeferencing errors (Maitner et al., 2018). For each georeferenced observation, we extracted the mean annual temperature (MAT, BIOCLIM1) and mean maximum daily temperatures of the warmest month (MTWM, BIOCLIM5) from the CHELSA extrapolated climate database (Karger et al., 2017) at a spatial resolution of 30 arc seconds (~1-km resolution at the equator). Any records with obvious georeferencing errors (e.g., falling in large bodies of water or outside of the neotropics) were discarded. We also masked the records to include only those from within the “Tropical & Subtropical Moist Broadleaf Forests” and “Tropical & Subtropical Dry Broadleaf Forests” (Dinerstein et al., 2017) in the Andean countries of Colombia, Venezuela, Ecuador, Peru, and Bolivia, in order to minimize potential errors and to reduce the influence of non-forest or non-Andean congeners when working at the genus level (see below). For each of the study species with ≥10 “clean” georeferenced records from Andean forests, we then calculated the geographic-based thermal optimum (GT_opt) as the median of the extracted MAT values (analyses were repeated using the mean MAT values but results did not differ qualitatively). We calculated the geographic-based thermal maximum (GT_max) as the 95% quantile of the extracted BIOCLIM5 values (the 95% quantile was used rather than the absolute maximum in order to help minimize the influence of outliers or remaining georeferencing errors). For morphospecies (n = 37) and species with <10 records (n = 52), we estimated the GT_opt and GT_max at the genus-level using all available records for congeners (Table S1). For a small number of species (n = 8), there were not sufficient records available from Andean forests to estimate GT_opt and GT_max at either the species- or genus-level; these species were excluded from relevant analyses. We also conducted all analyses using just the species for which GT_opt and GT_max could be estimated at the species-level but the results were not markedly changed (Table S3).

Relationships between variables (elevation, temperature, T₅₀, GT_opt and GT_max, ΔBA_rel) were assessed using ordinary linear least squares regressions. In the analyses of ΔBA_rel, species were weighted by the strength of the coefficient of determination (R²) of their BA_rel vs. census date relationship. We also used logistic regression to test if the probability of increasing in relative abundance (ΔBA_rel > 0) was related to T₅₀, GT_opt, or GT_max. All analyses were conducted in R version 3.6.1 (R Core Team, 2018).

Results

In accord with a priori predictions, the photosynthetic thermal tolerances (T₅₀) of tree species were positively related to plot temperature and elevation, such that the species sampled from hot lowland plots tended to have higher T₅₀ than species sampled from cold highland plots (T₅₀ = 49.72°-0.0004 × Elevation, R² = 0.07, P = 0.00101; T₅₀ = 47.39° + 0.0827 × MAT, R² = 0.07, P = 0.00115; T₅₀ = 48.16° + 0.0372 × MTWM, R² = 0.02, P = 0.0701; Figure 2). At the plot level, there was a tendency for the average T₅₀ (with species weighted by their relative basal area) to vary with temperature and elevation, but the relationships were not significant (mean T₅₀ = 49.92°-0.0004 × Elevation, R² = 0.24, P = 0.149; mean T₅₀ = 47.47° + 0.0866 × MAT, R² = 0.23, P = 0.156; mean T₅₀ = 48.04° + 0.0501 × MTWM, R² = 0.12, P = 0.331; Figure 2).

Figure 2. The relationships between species' thermal tolerances (T₅₀) and (A) plot elevation, (B) mean annual temperature (MAT), and (C) mean maximum temperature of the warmest month (MTWM). Circle points and lines show species-level data; square points show plot-level averages with the species weighted by their relative abundances, and triangle points show unweighted plot averages. T₅₀ = 49.72°-0.0004 × Elevation, R² = 0.07, P = 0.00101; T₅₀ = 47.39° + 0.0827 × MAT, R² = 0.07, P = 0.00115; T₅₀ = 48.16° + 0.0372 × MTWM, R² = 0.02, P = 0.0701.

Similarly, tree species sampled from the lowland plots tended to have geographic distributions that were centered on hotter areas (i.e., higher GT_opt) and that included hotter areas (i.e., higher GT_max) than species from highland plots (GT_opt = 25.72°-0.0026 × Elevation, R² = 0.45, P < 2 e⁻¹⁶; and GT_max = 31.29°-0.0012 × Elevation, R² = 0.2; P = 1.07e⁻¹⁰; Figure 3). The plot-level average of GT_opt and GT_max both decreased significantly with elevation (mean GT_opt = 26.14°-0.0028 × Elevation, R² = 0.82, P = 0.000282; and mean GT_max = 31.64°-0.0016 × Elevation, R² = 0.64, P = 0.00567; Figure 3). Plot-level relationships did not change if the species were not weighted by their abundances.

Figure 3. The relationships between plot elevation and species' (A) geographic-based thermal optima (Gt_opt), and (B) geographic-based thermal maximum (GT_max). Circle points and lines show species-level data; square points show plot-level averages with the species weighted by their relative abundances, and triangle points show unweighted plot averages. GT_opt = 25.72°-0.0026 × Elevation, R² = 0.45, P < 2e⁻¹⁶; GT_max = 31.29°-0.0012 × Elevation, R² = 0.2, P = 1.07 e⁻¹⁰.

There was a significant positive relationship between T₅₀ and GT_opt (GT_opt = −2.93° + 0.4931 × T₅₀, R² = 0.05, P = 0.00858; Figure 4A). However, there was not a significant relationship between T₅₀ and GT_max (GT_max = 20.95° + 0.1680 × T₅₀, R² = 0.01, P = 0.165037; Figure 4B).

Figure 4. The relationships between species' physiological maximum thermal tolerance (T₅₀) and (A) GT_opt and (B) GT_max. GT_opt = −2.93° + 0.4931 × T₅₀, R² = 0.05, P = 0.00858; GT_max = 20.95° + 0.1680 × T₅₀, R² = 0.01, P = 0.165037.

Contrary to predictions, there was no overall pattern for tree species with higher relative T₅₀ (T_50,relative = T₅₀-Plot Average T₅₀) to increase in relative abundance through time (ΔBA_rel = −3.506 e⁻⁰⁶ + 1.046 e⁻⁰⁵ × T_50,relative, R² = 0.01, P = 0.218; Figures 5A, 6). A significant positive relationship between thermal tolerances and changes in relative abundance was observed in only three of the plots (Maceo: ΔBA_rel = −3.125 e⁻⁰³ + 6.368 e⁻⁰⁵ × T₅₀, R² = 0.36, P = 0.009; Porce: ΔBA_rel = −2.835 e⁻⁰³ + 5.843 e⁻⁰⁵ × T₅₀, R² = 0.60, P = 0.005; Ventanas: ΔBA_rel = −3.994 e⁻⁰³ + 8.320 e⁻⁰⁵ × T₅₀, R² = 0.37, P = 0.021; Figure 6D). Likewise, while there were positive relationships between species' changes in relative abundance through time and both their relative GT_opt (GT_opt,relative = GT_opt-Plot Average GT_opt) and their relative GT_max (GT_max,relative = GT_max-Plot Average GT_max), these relationships were not significant (ΔBA_rel = −5.300 e⁻⁰⁶ + 9.521 e⁻⁰⁷ × GT_opt,relative, R² = 0.02, P = 0.0677; and ΔBA_rel = −7.945 e⁻⁰⁶ + 9.9960 e⁻⁰⁷ × GT_max,relative, R² = 0.02, P = 0.120; Figures 5B,C). The probability of species increasing in abundance was not significantly associated with any measure of thermal tolerance [logit(ΔBA_rel >0) = 0.166 + 0.00268 × T_50,relative, P = 0.982; logit(ΔBA_rel >0) = 0.172 + 0.0100 × GT_opt,relative, P = 0.169; and logit(ΔBA_rel >0) = 0.137 + 0.0162 × GT_max,relative, P = 0.0821].

Figure 5. The relationships between annual changes in species' relative abundance (ΔBA_rel) and (A) relative T₅₀, (B) relative GT_opt, and (C) relative GT_max (see text). ΔBA_rel = −3.506e⁻⁰⁶ + 1.046e⁻⁰⁵ × T_50,relative, R² = 0.01, P = 0.218; ΔBA_rel = −5.300 e⁻⁰⁶ + 9.521 e⁻⁰⁷ × GT_opt,relative, R² = 0.02, P = 0.0677; and ΔBA_rel = −7.945 e⁻⁰⁶ + 9.9960 e⁻⁰⁷ × GT_max,relative, R² = 0.02, P = 0.120.

Figure 6. The within-plot relationships between annual changes in species' relative abundance (ΔBA_rel) and species' T₅₀. Horizontal lines indicate no change in abundance through time and vertical lines indicate plot-average T₅₀. Solid regression lines indicate significant (p ≤ 0.05) relationships (weighted regressions with species weighted by the coefficient of determination of their BA_rel vs. census date relationship). Plots are ordered by increasing elevation: (A) Puerto Triunfo: ΔBA_rel = 1.001 e⁻⁰³-2.161 e⁻⁰⁵ × T₅₀, R² = 0.06, P = 0.335; (B) Támesis: ΔBA_rel = 5.779 e⁻⁰³-1.169 e⁻⁰⁴ × T₅₀, R² = 0.47, P = 0.339; (C), Maceo: ΔBA_rel = −3.125 e⁻⁰³ + 6.368 e⁻⁰⁵ × T₅₀, R² = 0.36, P = 0.009; (D) Porce: ΔBA_rel = −2.835 e⁻⁰³ + 5.843 e⁻⁰⁵ × T₅₀, R² = 0.60, P = 0.005; (E) Anorí: ΔBA_rel = −8.976 e⁻⁰⁴ + 1.844 e⁻⁰⁵ × T₅₀, R² = 0.13, P = 0.150; (F) Angelópolis: ΔBA_rel = 2.400 e⁻⁰³-5.094 e⁻⁰⁵ × T₅₀, R² = 0.42, P = 0.004; (G) Ventanas: ΔBA_rel = −3.994 e⁻⁰³ + 8.320 e⁻⁰⁵ × T₅₀, R² = 0.37, P = 0.021; (H) Jardín: ΔBA_rel = 1.261 e⁻⁰³-2.678 e⁻⁰⁵ × T₅₀, R² = 0.06, P = 0.345; (I) Caicedo: ΔBA_rel = −7.950 e⁻⁰⁴ + 1.718 e⁻⁰⁵ × T₅₀, R² = 0.02, P = 0.567; and (J) Belmira: ΔBA_rel = 3.782 e⁻⁰⁴-8.006 e⁻⁰⁶ × T₅₀, R² = 0.01, P = 0.766.

Discussion

In this study, we tested several hypotheses about how the thermal tolerances of tropical montane tree species relate to their geographic distributions and performances under rising temperatures. First, we tested the hypothesis that species from cold, highland forest communities have lower thermal tolerances compared to species from lowland forests where ambient air temperatures are reliably hotter. Our results support this hypothesis since T₅₀, GT_opt, and GT_max all increased with plot temperature and decreased with plot elevation. Likewise, the plot-level average T₅₀ tended to decrease with elevation (and increased with plot MAT), and plot-level averages of both GT_opt and GT_max both decreased significantly with elevation.

We found a very high amount of variation in the T₅₀ among species co-occurring within each plot such that elevation or temperature only explained a small, albeit highly significant, proportion of interspecific variation. A potential explanation for these patterns is that thermal tolerances are evolved in response to leaf temperatures and that leaf temperatures can be decoupled from regional air temperatures due to the fact that (1) even within a single plot, different species can occur predominantly in different thermal or moisture microhabitats and (2) co-occurring species can have different thermoregulatory behaviors and mechanisms.

Even over small areas (e.g., within single hectare inventory plots), there may exist many different microhabitats that create different thermal regimes for the leaves of co-occurring species. For example, leaves of canopy species will be exposed to more direct sunlight and thus reach higher temperatures than leaves of understory species (Smith, 1978; Rey-Sánchez et al., 2016; Slot et al., 2019). Air and leaf temperatures can also be strongly influenced by small-scale changes in topography (Geiger et al., 2009; Dobrowski, 2011; Sears et al., 2011; Curtis et al., 2016; Graae et al., 2018; Lembrechts et al., 2019). Many studies have found that tropical species are non-randomly distributed within plots. These studies, which have focused primarily on the relationships between the distributions of species and patterns of soil nutrients and soil water availability, have found good evidence that more drought-tolerant species occur in dryer microhabitats and that drought-intolerant species occur in wetter microhabitats (Harms et al., 2001; Valencia et al., 2004; Engelbrecht et al., 2007; Condit et al., 2013; Cosme et al., 2017; Zuleta et al., 2018). It is possible that species with higher or lower T₅₀ are likewise selectively distributed in hotter or cooler microhabitats, respectively. Additional studies are needed to investigate the distribution of tropical tree species in relation to thermal microhabitats.

Another factor that can influence leaf temperatures is the thermal regulatory strategies of the plants themselves. Different leaf sizes, shapes, morphologies, and behaviors can all lead to markedly different leaf temperatures of co-occurring plants (Parkhurst and Loucks, 1972; Meinzer and Goldstein, 1985; Leigh et al., 2012; Fauset et al., 2018). For example, all else being equal, larger leaves will reach higher maximum temperatures than will smaller leaves (Leigh et al., 2017; Wright et al., 2017). Likewise, leaves that absorb more solar radiation (rather than reflecting or transmitting the light) will reach higher maximum temperatures (Smith, 1978; Rey-Sánchez et al., 2016). Another important leaf thermoregulatory mechanism is evapotranspiration, which cools leaves through latent heat loss (Nobel, 1999; Lambers et al., 2008). Evapotranspirative cooling requires access to water and thus is generally most effective in species that have developed mechanisms for the rapid uptake and transport of water (both to the leaf via xylem transport and out of the leaf via stomata) and/or species that occur in wetter microhabitats with abundant soil water availability (Gates, 1968; Oren et al., 1999; Meinzer et al., 2008). Because of habitat selection and thermoregulatory mechanisms, leaf temperatures can be decoupled from air temperature and may differ significantly between species under identical air temperatures (Michaletz et al., 2015). By extension, this also means that species under markedly different air temperatures can have similar or identical leaf temperatures (Helliker and Richter, 2008; Michaletz et al., 2016). If thermal tolerances are indeed adapted or acclimated to leaf temperatures rather than air temperatures, then this decoupling may explain the high amount of variation between species within plots as seen in our study, as well as in other studies looking at patterns of thermal tolerances across latitude (O'sullivan et al., 2017).

Another hypothesis that we tested in this study was that tropical tree species with higher thermal tolerances have geographic ranges that encompass hotter regional air temperatures (Zhu et al., 2018). In accord with this hypothesis, we did find a significant positive relationship between species' T₅₀ and their thermal optima as based on their large-scale geographic distributions. However, T₅₀ explained very little variation between species' GT_opt and there was not a significant relationship between T₅₀ and species' geographic-based thermal maximum (GT_max). As discussed above, microhabitat selection and thermoregulatory strategies may decouple leaf temperatures from air temperatures—especially when air temperatures are estimated at the scale of 1-km² pixels. If this is true, then it may explain how some species with relatively low T₅₀ are able to include areas with high mean annual or maximum temperatures in their distributions and likewise why some species with high T₅₀ are geographically restricted to cooler areas. In order to better assess how thermal tolerances relate to geographic distributions, finer-scale climate models will need to be coupled with leaf biophysical models to map the distributions of leaf temperatures across species ranges (Kearney and Porter, 2017). Another possibility is that thermal tolerances are locally adapted or acclimated to plot conditions, which would obscure any relationships between T₅₀ and large-scale species-wide distributions (Zhu et al., 2018).

The final hypothesis that we tested was that because of global warming, tree species that are more-thermophilic should increase in abundance through time within each plot relative to less-thermophilic species. Our results did not support this hypothesis. We did not find consistent or significant relationships between the direction or rates of changes in species' relative abundances (BA_rel) and any measure of the species' thermal tolerances or preferences (T₅₀, GT_opt, GT_max). As above, the lack of a strong relationship between thermal properties and their changes in abundance may be due at least in part to the decoupling of air temperature and leaf temperatures. In other words, some species may be better able to regulate their leaf temperatures and thus may experience less changes in leaf temperatures over time than other co-occurring species. It is also possible that other factors could be influencing changes in species' abundances besides rising temperatures, thus obfuscating patterns. More specifically, there is large stochastic component to tree mortality, especially in montane forests where landslides and other local disturbances are very common (Clark et al., 2016). These stochastic events lead to a large amount of demographic “noise.” To overcome this noise, large-scale spatial or temporal datasets are required (Wagner et al., 2010). It is possible that this study was simply too short of duration to detect a directional signal in abundance changes, especially since we were only looking at a small subset of focal species within individual 1-ha plots.

Conclusions

Global warming makes it imperative that we understand the ability of species to tolerate and perform at hotter temperatures. Here we assessed the maximum physiological thermal tolerance of photosystem II (T₅₀) for nearly 200 species of tropical montane trees growing at different elevations along the flanks of the Colombian Andes Mountains. We used our estimates of T₅₀ to test if thermal tolerances vary predictably across an elevation gradient of more than 2,500 m (corresponding to a ca 20°C gradient in mean annual and maximum temperatures). Our results showed that there is in fact a negative relationship between elevation and T₅₀ such that species occurring at lower elevations are generally capable of tolerating hotter temperatures than species occurring at higher elevations. However, our results also indicate that the relationship between T₅₀ and temperature is weak and extremely shallow (Figure 2). Indeed, the slope of the T₅₀ vs. MAT relationship was only 0.08, and the slope of the T₅₀ vs. MTWM relationship was only 0.05. Our findings mirror those of O'sullivan et al. (2017) who likewise found a shallower-than-expected relationship between plant thermal tolerances and temperature (specifically, T_max = 49.155 + 0.264 × MTWM; unfortunately, differences in the methods used to estimate thermal tolerances prohibit a direct comparison of results).

The shallow relationships between plant thermal tolerances and plot temperatures across both elevational and latitudinal gradients suggest that species from hot low-elevation and low-latitude forests will have smaller average thermal safety margins (difference between temperature and thermal tolerance) than either high-latitude or high-elevation species, and consequently that lowland tropical plants may be at especially high risk of damage from rising global temperatures (Perez et al., 2016). However, we also found a large amount of variation in the T₅₀ of species co-occurring in any given site (i.e., some species growing in low-elevation forests where MAT = 27°C can have lower T₅₀ than species from high elevation forests where MAT = 10°C, and vice versa) and also very little relation between species' thermal tolerances and their large-scale geographic distributions (i.e., some species with low T₅₀ have ranges that include hotter areas than species with higher T₅₀). Consequently, many lowland species may potentially have larger thermal safety margins than highland species, making it hard to make generalized predictions about risk.

One possible explanation for the high variation in thermal tolerances of co-occurring species is that T₅₀ (and other comparable measures) is evolved in response to extreme leaf temperatures that can be decoupled from air temperatures. A decoupling of air and leaf temperatures, and the potential ability of some plants to maintain relatively constant leaf temperatures despite global warming, could also explain why there was no significant relationship between species' T₅₀ and changes in species abundances over time. The thermal regulatory behaviors of tropical plants, and the ability of these behaviors to respond to climate change, remains poorly studied and deserves further attention.

Data Availability Statement

All datasets generated for this study are included in the article/Supplementary Material.

Author Contributions

KF and AD designed the study. JM-V, AS, and DT collected the field data. KF and TP analyzed the data. KF wrote the manuscript. KF, TP, and AD edited the manuscript.

Funding

This work was funded through grants from the Swiss Agency for Development and Cooperation, and the U.S. National Science Foundation (NSF, DEB-1350125).

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/ffgc.2020.00025/full#supplementary-material

Figure S1. Species-level changes in F_V/F_M recorded in leaves exposed to different temperature treatments and model fits used to estimate T₅₀.

Table S1. Study species, measures of thermal tolerance (T50, GTopt, GTmax), changes in abundance.

Table S2. Results of analyses using only species for which GT_opt and GT_max could be estimated at the species-level.

Table S3. Basal area (cm²) of focal species and total tree basal area recorded in each census of the 10 study plots.

References