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

Due to global warming, many species will face greater risks of thermal stress, which can lead to changes in performance, abundance, and/or geographic distributions. In plants, high temperatures above a species-specific critical thermal maximum will permanently damage photosystem II, leading to decreased electron transport rates, photosynthetic failure, and eventual leaf and plant death. Previous studies have shown that plant thermal tolerances vary with latitude, but little is known about how they change across smaller-scale thermal gradients (i.e., with elevation) or about how these thermal tolerances relate to species' local performances and geographic distributions. In this study, we assess the maximum photosynthetic thermal tolerances (T50) of nearly 200 tropical tree species growing in 10 forest plots distributed across a >2,500 m elevation gradient (corresponding to a 17°C temperature gradient) in the northern Andes Mountains of Colombia. Using these data, we test the relationships between species' thermal tolerances and (1) plot elevations and temperatures, (2) species' large-scale geographic distributions, and (3) changes in species' abundances through time within the plots. We found that species' T50 do in fact decrease with plot elevation but significantly slower than the corresponding adiabatic lapse rate (−0.4 vs. −5.7°C km−1) and that there remains a large amount of unexplained variation in the thermal tolerances of co-occurring tree species. There was only a very weak association between species' thermal tolerances and their large-scale geographic distributions and no significant relationships between species' thermal tolerances and their changes in relative abundance through time. A potential explanation for these results is that thermal tolerances are adaptations to extreme leaf temperatures that can be decoupled from regional air temperatures due to microclimatic variations and differences in the species' leaf thermoregulatory properties.


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(Feeley et al., , 2013Duque 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 largescale 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.

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 −1 , 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 −1 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).

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 −1 or larger than 7.5 cm year −1 , 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., morethermophilic 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(Feeley et al., , 2013Fadrique 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 50 ) 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 50 were based on standard protocols modified from Krause et al. (2010Krause et al. ( , 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 0 ) 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 0 ) 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 −2 s −1 ) 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 50 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 /(1 + exp(-(θ 2 + θ 3 × Temperature)))] where θ 1 is the control value of F V /F M (≈ 0.8) and θ 2 and θ 3 are the intercept and slope coefficients of the logit(Fv/Fm) ∼ Temperature relationship, respectively; R Core Team, 2018). T 50 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 50 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 50 . For some species (n = 18), estimates of T 50 were deemed unreliable due to unreasonable values (i.e., T 50 >55 or <40 o C) 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 50 .

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 50 ; 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 −1 ; 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 50 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 genuslevel 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 50 , 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 2 ) 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 50 , GT opt , or GT max. All analyses were conducted in R version 3.6.1 (R Core Team, 2018).

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 50 , GT opt , and GT max all increased with plot temperature and decreased with plot elevation. Likewise, the plot-level average T 50 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 50 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 50 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 50 and their thermal optima as based on their large-scale geographic distributions. However, T 50 explained very little variation between species' GT opt and there was not a significant relationship between T 50 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 2 pixels. If this is true, then  it may explain how some species with relatively low T 50 are able to include areas with high mean annual or maximum temperatures in their distributions and likewise why some species with high T 50 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 50 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 lessthermophilic 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 50 , 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 50 ) 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 50 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 50 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 50 and temperature is weak and extremely shallow (Figure 2). Indeed, the slope of the T 50 vs. MAT relationship was only 0.08, and the slope of the T 50 vs. MTWM relationship was only 0.05. Our findings mirror those of O'sullivan et al. (2017) who likewise found a shallowerthan-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 lowlatitude 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 50 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 50 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 50 have ranges that include hotter areas than species with higher T 50 ). 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 50 (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 50 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).

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 50 .
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 2 ) of focal species and total tree basal area recorded in each census of the 10 study plots.