DOI: 10.5061/dryad.rbnzs7hc3
title: Towards an understanding of the latitudinal patterns in thermal tolerance and vulnerability of woody plants under climate warming

description: Predicting spatial patterns in thermal tolerance and vulnerability of
 species under climate warming remains a challenge. Current knowledge is
 mainly from experiment-based thermal physiology of limited numbers of
 ectotherms, yet large-scale evaluations on plants remain elusive. Here,
 using distribution maps with spatial resolutions of 20×20 km for 5628
 woody species in China, we propose a novel approach, i.e. thermal
 distribution curves, to describe species' realized thermal niches,
 and then estimate their thermal tolerance and warming risks under
 projected climate warming in 2050s and 2070s. We find that species'
 vulnerability and potential local extinction risks within grid cells
 decrease with latitude and increase with aridity due to narrow thermal
 tolerance of species located at low latitudes and arid regions. Over 90%
 of species could still tolerate future warming in most areas, indicating
 relatively optimistic expectation of potential local extinctions. Our
 study presents a new framework to quantify climate warming impacts on a
 large number of species without sufficient physiological information, and
 provides fundamental references for conservation planning under climate
 change.
description: Species distribution
 data The distribution data of
 woody species were first obtained from the <i>Atlas of Woody Plants
 in China: Distribution and Climate </i>(Fang et al. 2011), which
 provides the best collection of county-level distributions of all woody
 species published before 2009 (<i>see </i>Appendix S1 for more
 details). We then substantially updated the distribution data using
 national and regional floras published between 2009-2018, peer-reviewed
 papers (<i>https://www.cnki.net/</i>, accessed in April 2018),
 and recently published specimen records
 (<i>http://www.nsii.org.cn/</i>, accessed in April 2018). A
 full list of the data sources is available in Appendix S4. To retain the
 accuracy of species distributions, we only compiled species distributions
 at county level or finer scales (e.g. villages, towns, etc.). Most
 counties in China are relatively small, while those located in the
 Taklimakan Desert and the northern Qinghai-Tibetan Plateau are large. To
 improve the accuracy of species distributions in these large counties,
 each of them was divided into 2-4 units during data compilation and
 species distributions in these counties were further identified to these
 small units using the information of species habitats and elevation
 ranges. The species names from different sources were standardized
 following <i>Flora of China</i>. To further refine the
 species distributions, we compiled the elevational ranges and living
 habitats of each species from <i>Flora of China</i> and other
 national and provincial floras. The habitat types of each species were
 matched with the 47 natural vegetation types recorded in the 1:1,000,000
 vegetation map of China (Editorial Committee of Vegetation Map of China
 2007). Then by integrating the county-level distribution of each species
 with the vegetation map and a digital elevation model (spatial resolution:
 1 km), we refined the county-level distributions of each species into an
 equal-area grid with a spatial resolution of 20×20
 km<sup>2</sup>. After removing the incomplete grid cells with
 less than 200 km<sup>2</sup> (i.e. 50% area of a grid cell)
 within national borders and coasts, a total of 23718 grid cells were
 remained for the following analyses. Since the performance of
 ENMs is largely affected by sample size of species occurrence, species
 with less than 20 occurrence records were not included following previous
 studies (Xu et al. 2019). The distributions of non-endemic species in
 China may not reflect their true thermal niches. To estimate species'
 niche more precisely, species whose major distribution ranges are not in
 China were not included for ENM calibrations. Specifically, we compiled
 the global distributions of Chinese woody species at a spatial resolution
 of ca. 4<sup>o</sup>×4<sup>o</sup>
 (<i>see</i> Appendix S5 for a list of data sources). Then we
 used the student-t test to examine the difference between species'
 temperature ranges (defined as the difference between the maximum and
 minimum temperature within a species distribution range) in China and the
 whole world. Species were excluded from the following analyses if their
 temperature ranges in China were smaller than 75% of those at the global
 scale. Finally a total of 5628 species were used for ENM
 calibrations. Current and future
 climate data Current climate at a
 spatial resolution of 2.5′ × 2.5′ was obtained from the
 <i>WorldClim</i>, including 19 bioclimatic variables. To
 reduce the multicollinearity among climate variables, we calculated
 Pearson's correlation coefficients (r) between each pair of variables
 and removed one from each pair with r &gt; 0.7. Finally, six
 bioclimatic variables were selected: mean annual temperature (MAT, bio1),
 mean diurnal range (bio2), temperature seasonality (bio4), precipitation
 seasonality (bio15), precipitation of warmest quarter (bio18) and
 precipitation of coldest quarter (bio19). Several temperature variables
 (e.g. MAT, Deutsch et al. 2008; summer temperature, Sentinella et al.
 2020) could be used to construct thermal response curves. MAT is highly
 correlated with growing season or summer temperature (Pearson's r
 &gt; 0.9). Previous studies suggest that MAT provides superior and
 robust estimations of species' thermal tolerance compared with other
 temperature measures (Deutsch et al. 2008, Diamond et al. 2012).
 Therefore, we used MAT to estimate species' thermal niches and
 warming impacts in our study. Future MAT in the 2050s
 (2041-2060) and 2070s (2061-2080) generated by General Circulation Models
 (GCM) were also obtained from the <i>WorldClim </i>database.
 Here, we used data projected by three GCMs: BCC-CSM1-1.1, MIROC-ESM and
 CCSM4. Previous studies suggest<b> </b>that<b>
 </b>BCC-CSM1-1.1 model could better reflect the magnitude of climate
 changes in China relative to other GCMs (Xin et al. 2013). Moreover, the
 results based on the three GCMs are highly consistent with each other
 (<i>see</i> Table S1 in Appendix S3). We hence reported the
 results based on BCC-CSM1-1.1 in the main text. For each GCM, we included
 three representative concentration pathways (RCP2.6, RCP6.0 and
 RCP8.5). ENM calibration and
 evaluation The
 relationships between species occurrences and the six selected climate
 variables were evaluated using two ENM algorithms: 1) generalized linear
 model (GLM) and 2) Maximum entropy model (MaxEnt) (Phillips<i>
 </i>et al. 2006). These two algorithms have been widely used and
 offer relatively good accuracy for modelling the relationship between
 species' occurrences and climate. For each species, we randomly chose
 80% of its distribution records for model calibration and the remaining
 20% for model evaluation, and this was repeated for 10 times for each ENM
 algorithm separately (i.e. 10 models). The true skill statistics (TSS) was
 used to evaluate the performance of these models, and those with TSS
 &gt; 0.5 (Allouche et al. 2006) were used 1) to estimate the thermal
 distribution curves (<i>see</i> the following section) and 2)
 to predict the occurrence probabilities of each species per grid cell
 under the current climate. Then for each species, the average occurrence
 probabilities per grid cell weighted by model TSS were estimated for each
 ENM algorithm separately, and were used to generate its predicted current
 binary distribution using the threshold that maximizes TSS values (Xu et
 al. 2019). The spatial
 extent for model calibrations and predictions of each species was the
 terrestrial range of China. ENMs tend to over predict species
 distributions due to the ignorance of the effects of dispersal limitation
 (Raxworthy et al. 2003). Therefore, a 200 km buffer around the observed
 range of each species was used to clip the predicted current distributions
 following previous studies (Kremen et al. 2008, VanDerWal et
 al<i>.</i> 2009). Both the actual and predicted species
 distributions were used to infer thermal niches of species, and the
 results were highly consistent (<i>see</i> below). We
 presented the results based on the actual distribution data in the main
 text. The model calibrations, evaluations and predictions were all
 conducted using BIOMOD2 package in R (Thuiller et al. 2009). Except the
 above-mentioned parameters, we adopted the default settings for both
 algorithms. Thermal distribution
 curves of species The response curves
 between the probability of occurrence and MAT after controlling for the
 effects of other predictors were estimated for each species using the ENMs
 with TSS &gt; 0.5 (Guisan and Zimmermann 2000), and were then averaged
 across the models of each ENM algorithm using model TSS as weights.
 Although previous studies on animal taxa suggest that thermal performance
 curves are often skewed (Deutch et al. 2008, Waldock et al. 2019), our
 preliminary analyses indicated that the thermal distribution curves of
 most species were symmetric (Fig. S1-S2 in Appendix S1), and hence could
 be estimated using quadratic functions (Diamond<i> </i>et al.
 2012, Sentinella et al<i>.</i> 2020). Specifically, we used
 the quadratic to estimate the thermal distribution curves and
 <i>AT<sub>opt</sub></i> and
 <i>AT<sub>max</sub></i> of each species
 (<i>see</i> Deutsch et al<i>.</i>
 2008): where
 <i>P(AT)</i> represents the probability of occurrence of
 species at given ambient temperature (<i>AT</i>). For each
 species, the range of <i>AT</i> (i.e. function domain) was
 determined by the MAT range within its actual and predicted current
 geographical distributions, separately. We also extracted the maximum
 temperatures (<i>AT'<sub>max</sub></i>)
 within species distribution ranges and compared them to the estimated
 thermal limits
 (<i>AT<sub>max</sub></i>). Importance values of different climatic variables in ENMs were estimated to evaluate their contributions in determining species distributions (Thuiller et al. 2009, Dyderski et al. 2018). Low importance values of MAT indicated that temperature was not a critical factor, and hence species distributions might not be sensitive to ambient temperature changes. To avoid overestimation of extinction risks induced by warming in the following analyses, we did not include species for the estimation of thermal distribution curves if the importance value of MAT was &lt;0.3 (or the rank of MAT importance value was &gt;3) in more than half of its models (n &gt; 5) based on each ENM algorithm (Thuiller et al<i>.</i> 2009). We also excluded species from the following analyses for calculating <i>AT<sub>opt</sub></i> and<i> AT<sub>max</sub></i> when their thermal distribution curves were (a) concave up (i.e. the probability of occurrence is the lowest at intermediate MAT; 134 species for GLM and 176 species for MaxEnt) or (b) incomplete (e.g. lacking a peak of occurrence probability; 385 species for GLM and 428 species for MaxEnt). Finally, thermal distribution curves were successfully estimated for 4239 and 3128 species based on GLM and MaxEnt, respectively. Although the mean importance value of MAT in GLMs was higher than that in MaxEnt (Fig. S1 in Appendix S2), the difference in the species numbers between the two algorithms did not bias our conclusions. Latitudinal patterns in warming impacts on species vulnerability and Local extinction risk Thermal safety margin (TSM) and warming tolerance (WT) have been previously used to evaluate the latitudinal patterns of warming impacts on the local extinction risks of terrestrial ectotherms (e.g. Deutch et al. 2008), endotherms (Khaliq et al. 2014) and marine fauna (e.g. Sunday et al. 2014, Stuart-Smith et al. 2015). Here, we applied these two metrics to the thermal distribution curves of woody plants. Specifically, let <i>AT<sub>hab_current</sub></i><sub> </sub>denote the current ambient MAT within each grid cell. TSM and WT were calculated as: TSM = <i>AT<sub>opt</sub></i> - <i>AT<sub>hab_current</sub></i>; WT = <i>AT<sub>max</sub></i> - <i>AT<sub>hab_current</sub></i>. Species with low TSM normally live in environments that are already close to their <i>AT<sub>opt</sub></i>, and hence small amount of warming could decrease their suitability. WT of a species provides an estimate of how much warming it could tolerate in a given habitat before its habitat suitability drops to zero. To demonstrate the geographic patterns in TSM and WT, we estimated the mean TSM and WT of each grid cell (i.e. community) as the average of all species within it. Then the values within each 1<sup>o</sup> latitudinal belt were averaged to explore the latitudinal gradient in TSM and WT using simple linear regressions. It is noteworthy that we calculated the mean TSM and WT per grid cell without considering species abundance due to lacking of data. Previous studies suggested that species abundance may not change the nature of the overall large-scale patterns in community mean temperature index (Stuart-Smith et al. 2015). To predict the latitudinal gradient in the vulnerability and local extinction risk of species within grid cells, we first estimated the warming risk (WR) of each species as the ratio of the magnitude of future climate warming to its WT in a given grid cell: <i>WR</i> = (<i>AT<sub>hab_future</sub></i> - <i>AT<sub>hab_current</sub></i>)/<i>WT</i>, where <i>AT<sub>hab_future</sub></i> represented the future ambient MAT of a grid cell. Species may become vulnerable to future warming in a grid cell if their WR was high (a threshold of 0.5 was arbitrarily used, representing half of WT), and may be exposed to potentially unsuitable temperature (i.e. population will experience sharp declines) and hence have high risk to go locally extinct if their WR was &gt;1.0 (termed as potential local extinction risk). The proportions of species with WR &gt; 0.5 (hereafter CWR<sub>0.5</sub>) and WR &gt; 1.0 (CWR<sub>1.0</sub>) in each grid cell were estimated to represent the magnitude of species vulnerability and potential local extinction risks within grid cells, respectively. Spatial patterns of CWR<sub>0.5</sub> and CWR<sub>1.0</sub> were then mapped in ArcGIS 10.0 (ESRI, Redlands, CA), and the mean CWR<sub>0.5</sub> and CWR<sub>1.0</sub> within each 1<sup>o</sup> latitudinal belt were estimated to evaluate their latitudinal gradients. We excluded grid cells containing less than five species to eliminate the influences of sample size. Studies suggest that dryland plants may be more vulnerable to climate warming than those living in humid regions (e.g. Suzuki et al. 2014). We therefore compared the mean CWR<sub>0.5</sub> and CWR<sub>1.0</sub> between regions with mean annual precipitation higher or lower than 200 mm. To examine the evolutionary dependence of warming risks across species, we also evaluated the phylogenetic signals in the proportions of species distribution ranges with WR &gt; 0.5 and WR &gt;1.0 (Appendix <i>S1</i>). All statistical analyses were conducted in R 3.5.3 (R Core Team 2018).