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.

Species distribution data

The distribution data of woody species were first obtained from the Atlas of Woody Plants in China: Distribution and Climate (Fang et al. 2011), which provides the best collection of county-level distributions of all woody species published before 2009 (see Appendix S1 for more details). We then substantially updated the distribution data using national and regional floras published between 2009-2018, peer-reviewed papers (https://www.cnki.net/, accessed in April 2018), and recently published specimen records (http://www.nsii.org.cn/, 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 Flora of China.

To further refine the species distributions, we compiled the elevational ranges and living habitats of each species from Flora of China 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². After removing the incomplete grid cells with less than 200 km² (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^o×4^o (see 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 WorldClim, 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 > 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 > 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 WorldClim database. Here, we used data projected by three GCMs: BCC-CSM1-1.1, MIROC-ESM and CCSM4. Previous studies suggest that 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 (see 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 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 > 0.5 (Allouche et al. 2006) were used 1) to estimate the thermal distribution curves (see 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. 2009). Both the actual and predicted species distributions were used to infer thermal niches of species, and the results were highly consistent (see 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 > 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 et al. 2012, Sentinella et al. 2020). Specifically, we used the quadratic to estimate the thermal distribution curves and AT_opt and AT_max of each species (see Deutsch et al. 2008):

where P(AT) represents the probability of occurrence of species at given ambient temperature (AT). For each species, the range of AT (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 (AT'_max) within species distribution ranges and compared them to the estimated thermal limits (AT_max).

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 <0.3 (or the rank of MAT importance value was >3) in more than half of its models (n > 5) based on each ENM algorithm (Thuiller et al. 2009). We also excluded species from the following analyses for calculating AT_opt and AT_max 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 AT_{hab_current} denote the current ambient MAT within each grid cell. TSM and WT were calculated as: TSM = AT_opt - AT_{hab_current}; WT = AT_max - AT_{hab_current}. Species with low TSM normally live in environments that are already close to their AT_opt, 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^o 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: WR = (AT_{hab_future} - AT_{hab_current})/WT, where AT_{hab_future} 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 >1.0 (termed as potential local extinction risk). The proportions of species with WR > 0.5 (hereafter CWR_0.5) and WR > 1.0 (CWR_1.0) 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_0.5 and CWR_1.0 were then mapped in ArcGIS 10.0 (ESRI, Redlands, CA), and the mean CWR_0.5 and CWR_1.0 within each 1^o 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_0.5 and CWR_1.0 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 > 0.5 and WR >1.0 (Appendix S1).

All statistical analyses were conducted in R 3.5.3 (R Core Team 2018).