Materials and methods

Plant richness data. Data on plant richness, including that of ferns (non-seed) and seed plants, including gymnosperms and angiosperms, were collected and organized across different taxonomic levels (i.e., the total number of families, genera, and species) from published records regarding natural reserves in China. The list of published records was provided by the study by ^{Wang et al. (2017}, Figure 1). Based on the published records, we extracted data on the plant taxon richness found in core areas of protected areas. In China, the goal of core zones of protected areas is to protect relatively undisturbed natural vegetation, which has a long, uninterrupted history in this region, and the associated data thus represent an ideal dataset of plant richness (^{Tang et al. 2010}, ^{Zhang et al. 2017}).

Figure 1 Distribution of the protected areas in China used in our analysis. 

Based on previous studies (e.g., ^{Huang et al. 2016}, ^{Feng et al. 2017}, ^{Zhang et al. 2017}, ^{Liu et al. 2018}), we transposed the plant richness data from the core zones of the protected areas into grid data at a spatial resolution of 10 arc minutes (c. 16 × 16 km). Different groups of vascular plants (fern, gymnosperm, and angiosperm plants) were analysed in our study, as vascular plants include non-seed plants (ferns) and seed plants (gymnosperms and angiosperms). For the accurate nomenclature of scientific names, we followed the Plant List (http://www.theplantlist.org) and compared the lists of families, genera, and species based on The Plant List (http://www.theplantlist.org) and the Flora of China (http://frps.iplant.cn/) to identify the plant groups in our study. We deleted the wrong data on plant richness of protected areas. We found that the effect of the area of a nature preserve on plant richness could be excluded from further analyses because there was no significant relationship between reserve size and plant richness across the different taxonomic levels (i.e., the total number of families, genera, and species) based on linear regression modelling (P > 0.05). Finally, data from protected areas were obtained (detailed information in the Supplemental data and Figure 1).

Climate data. The MAT and MAP were used to assess the legacy effects of climate on the large-scale distributional pattern of plant richness (e.g., ^{Svenning et al. 2015}, ^{Feng et al. 2017}, ^{Blonder et al. 2018}). ^{Feng et al. (2017)} showed that historical and contemporary MAT and MAP could influence plant endemism in China. We downloaded the grid maps of historical and contemporary MAT and MAP at a spatial resolution of 10 arc minutes (ca. 16 × 16 km) from the WorldClim database (http://www.worldclim.org/). The extremely dry and cold climate during the Last Glacial Maximum (LGM; approximately 22,000 years ago) excluded tropical forests from China and caused other strong vegetational changes (^{Wang et al. 2012}, ^{Feng et al. 2017}). Paleoclimate (i.e., the LGM) has been shown to be the main driver of plant richness at large scales (e.g., ^{Kreft & Jetz 2007}, ^{Svenning & Skov 2007}, ^{Fang et al. 2012}, ^{Svenning et al. 2015}). The Holocene has not been long enough to have allowed speciation in most cases (^{Lister 2004}, ^{Svenning et al. 2015}); hence, we used average climate data from 1950-2000 AD to represent the contemporary scenario and paleoclimate data (i.e., the LGM) for the historical scenario. Paleoclimate data in regard to MAT and MAP were obtained from the CCSM4 general circulation model (http://www.cesm.ucar.edu/models/ccsm4.0/). The CCSM4 model consists of a coupled atmospheric, oceanic, and sea ice model with noninteractive vegetation and an atmospheric resolution of 10.0 arc minutes. The model is driven by variations in orbital configuration, greenhouse, ice-sheet topography, and coincident sea level changes and bathymetry for paleoclimates (^{Lawrence & Oleson 2012}). The paleoclimate data have the same coordinate system and resolution as the contemporary climate data. A paired-sample t-test coupled with a Bonferroni adjustment was used to evaluate the differences between the paleoclimate data and contemporary climate data across all the protected areas. The paired-sample t-test was conducted in JMP version 11.0 (SAS Institute Inc., Cary NC).

Data analysis. Spatial autocorrelation in ecological data can inflate Type I errors in statistical analyses (^{Diniz-Filho et al. 2003}). Hence, we used Moran’s I coefficient calculated on the basis of a distance matrix to assess the spatial autocorrelation in plant richness across the different taxonomic levels (i.e., the total number of families, genera, and species) (^{Diniz-Filho et al. 2003}). The default settings were used in SAM 4.0 (^{Rangel et al. 2010}), and the default number of distance classes was 17 with an equal number of pairs between different protected areas according to the available plant richness data (^{Rangel et al. 2010}). Pearson correlation coefficients were used to assess the correlations in plant richness between the different groups (i.e., vascular, fern, seed, gymnosperm, and angiosperm) and taxonomic levels (i.e., family, genus, and species) across the protected areas. The analysis of Pearson correlation coefficients was conducted in JMP version 11.0 (SAS Institute Inc., Cary NC).

Then, we used geographically weighted regression (GWR) coupled with the ordinary least squares (OLS) method to evaluate both the historical and contemporary climates and identify the large-scale distributional patterns of plant richness (^{Brunsdon et al. 1996}, ^{Mellin et al. 2014}, ^{Xu et al. 2016}). GWR is a local form of linear regression that is used to spatially model varying relationships based on the assessment of nonstationarity and the effects of spatial scale on ecological data (^{Brunsdon et al. 1996}, ^{Mellin et al. 2014}). Previous studies (e.g., ^{Foody 2004}, ^{Eiserhardt et al. 2011}, ^{Mellin et al. 2014}) have shown that GWR is useful in the investigation of spatially varying biodiversity-environment relationships because spatial autocorrelation and heterogeneity exist in ecological data. The variables pertaining to the historical and contemporary climates (including MAT and MAP) were regarded as explanatory variables, and plant richness was regarded as the dependent variable across the different groups of vascular, fern, seed, gymnosperm, and angiosperm plants. The specific settings for the GWRs were as follows: 1) the spatial function of the GWR was Gaussian; 2) the adaptive kernel was 10 % neighbours; and 3) optimization to minimize the AICc (corrected Akaike information criterion) was used for all bandwidths (^{Brunsdon et al. 1996}, ^{Eiserhardt et al. 2011}, ^{Xu et al. 2016}). The GWR analysis was conducted with SAM 4.0 (^{Rangel et al. 2010}).

The correlation coefficients (r) and P-values from the GWR and OLS analyses were used to assess the associations between the historical and contemporary climates (including MAT and MAP) and plant richness. We used the adjusted R²_adj (%) from the GWR to determine the explanatory power of the historical and contemporary climate in regard to the large-scale distributional pattern of plant richness (^{Blonder et al. 2018}, ^{Liu et al. 2018}). Meanwhile, the R²_adj (%) of the OLS analysis was used to test the relationships between the historical and contemporary climates and plant richness. Then, we compared the correlation coefficient (r) and R²_adj (%) of the GWR with the OLS analysis to test whether the GWR performed better than the OLS method (^{Brunsdon et al. 1996}).

Finally, we used OLS to determine the best predictors of the large-scale distributional patterns of plant richness at different taxonomic levels (i.e., the total number of species, families, and genera) across the vascular, fern, seed, gymnosperm, and angiosperm plant groups in independent analyses (^{Nagelkerke 1991}, ^{Liu et al. 2018}). The adjusted R²_adj (%) from the OLS analysis was used to determine the explanatory power of climate in regard to the large-scale distributional patterns of plant richness. We conducted the OLS analysis in JMP 10.0 (SAS Institute Inc., Cary, NC).

Results

The ranges of vascular plant richness were 11-257, 30-1372, and 45-4543 at the family, genus, and species levels, respectively (Table 1). The average vascular plant richness was 132, 506, and 1117 from the family to species level (Table 1). Specifically, the average fern species richness was 20, 40, and 87 at the family, genus, and species levels, respectively, and the average seed plant species richness was 112, 467, and 1027 at the family, genus, and species levels, respectively (Table 1). The average family, genus, and species richness values were 4 (ranging from 1 to 10), 8 (ranging from 1 to 33) and 14 (ranging from 1 to 102) for gymnosperms and 108 (ranging from 22 to 203), 460 (ranging from 52 to 1,244) and 1,024 (ranging from 65 to 3,931) for angiosperm plants, respectively (Table 1).

The MAP and MAT in the contemporary climate were significantly higher than those in the paleoclimate (paired-sample t-test; P < 0.05). Specifically, the average historical MAT was 5.1 °C (ranging from -14.8 °C to 21.9 °C), and the contemporary MAT was 10.1 °C (ranging from -8.9 °C to 25.5 °C; Table 2). The average historical MAP was 830.7 mm (ranging from 17 mm to 2232 mm), and the contemporary MAP was 952.4 mm (ranging from 26 mm to 2,262 mm; Table 2).

Based on the Moran’s I coefficients, the spatial autocorrelation in the plant richness data was low across the different taxonomic levels (most Moran’s I coefficients were < 0.200 or > -0.200; Figure 2). We found that there were significant correlations in plant richness among the different taxonomic levels (i.e., families, genera, and species; P < 0.05; Table 3). Additionally, a significant relationship of plant richness among the different vascular, fern, seed, gymnosperm, and angiosperm plant groups could be detected (P < 0.05; Table 3). The correlation coefficients were the largest (0.9979; P < 0.05) between seed and gymnosperm plants at the family level, between vascular and seed plants at the genus level (0.9979; P < 0.05), and between seed and angiosperm plants (0.9979; P < 0.05) at the species level (Table 3).

Figure 2 Moran’s I coefficients across different taxonomic levels (i.e., families, genera, and species) and groups (i.e., vascular, fern, seed, gymnosperm, and angiosperm plants). 

All the correlation coefficients in the GWR between the historical and contemporary climate variables (including MAT and MAP) and plant richness were higher than 0.49 (P < 0.001), and those in the OLS analysis were higher than 0.20 across the different taxonomic levels (P < 0.001; Table 4). Values of R²_adj (%) represent the adjusted R²_adj (%) values obtained from the GWR and OLS analyses conducted to determine the explanatory power of historical and contemporary climate variables on the large-scale distributional pattern of plant richness. Based on R²_adj, the GWR showed better modelling performance than the OLS analysis, indicating that both historical and contemporary climate coupled with spatial autocorrelation can explain the distributional patterns of plant richness at large scales (Table 4).

The following results regarding the R²_adj value were obtained from the GWR. The combination of historical and contemporary climate could explain the richness of vascular plant families, genera, and species (R²_adj = 25.8 %, 26.1 %, and 26.4 %, respectively; P < 0.001; Table 4), and both historical and contemporary climate had the strongest explanatory power in regard to the richness of fern genera and species (R²_adj = 34.6 % and 37.6 %, respectively; P < 0.001; Table 4). The historical and contemporary MAT and MAP had the strongest explanatory power in regard to richness at the family level (P < 0.001; Table 4) but the smallest explanatory power in regard to species richness based on the R²_adj (%) from the OLS analysis (P < 0.001; Table 4).

Historical MAT explained plant richness at different taxonomic levels in a better way than contemporary MAT, and this explanatory power of the MAP was opposite to that of the MAT across the different taxonomic levels (Figure 3). Furthermore, the explanatory power of the contemporary MAP was the strongest for plant richness at the family, genus, and species levels (Figure 3). With the exception of historical and contemporary MAT, historical and contemporary climate variables better explained fern richness than seed plant richness across all the taxonomic levels (Figure 3). Regarding gymnosperm and angiosperm plants, the explanatory power of historical and contemporary climates was stronger for gymnosperm plant richness than angiosperm plant richness across all taxonomic levels (Figure 4). Furthermore, the historical and contemporary MAP had the strongest explanatory power in regarding to angiosperm plant richness at the family level (R²_adj > 60.0%; P < 0.001; Figure 4).

Figure 3 Changes in the large-scale richness of vascular, fern, and seed plants as a function of historical and contemporary climate across different taxonomic levels. The lines associated with the points in each panel indicate the relationships between historical and contemporary climate and plant richness across different taxonomic levels (i.e., family, genus, and species). Shaded areas over the dashed regression lines represent the 95 % confidence interval of the fitted values for each evaluated order. Values of R² represent the adjusted R²_adj (%) values from the OLS analysis, which indicate the explanatory power of historical and contemporary climate variables in regard to the large-scale distributional pattern of plant richness. We transferred the values of R²_adj to Figure 3 based on the unit “%". Plant richness represents the number of families, genera, and species. All the R²_adj (%) values obtained from the OLS analysis were significant (P < 0.05). 

Figure 4 Changes in the large-scale richness of gymnosperm and angiosperm plants as a function of historical and contemporary climate across different taxonomic levels. The lines associated with the points in each panel indicate the relationships between historical and contemporary climate and plant richness across different taxonomic levels (i.e., family, genus, and species). Shaded areas over the dashed regression lines represent the 95 % confidence intervals of the fitted values for each evaluated order. Values of R² represent the adjusted R²_adj (%) values obtained from the OLS analysis, which indicate the explanatory power of historical and contemporary climates in regard to the large-scale distributional pattern of plant richness. We transferred the values of R²_adj to Figure 4 based on the unit “%". Plant richness represents the number of families, genera, and species. All the R²_adj (%) values obtained from the OLS analysis were significant (P < 0.05). 

Discussion

Our results showing that the average vascular plant richness was 132 (ranging from 11 to 257), 506 (ranging from 30 to 1,372), and 1,117 (ranging from 45 to 4,543) at the family, genus, and species levels, respectively, and that the plant richness was also high for fern, seed, gymnosperm, and angiosperm plants in protected areas of China (detailed information in Table 1) indicate that the protected area network in China has rich plant resources. In particular, Chinese protected areas have a high richness of angiosperm plants (Table 1). Hence, the ability of protected areas to conserve plant richness is strong in China.

We found that the contemporary climate variables (i.e., MAP and MAT) were significantly higher than those for the paleoclimate (paired-sample t-test; P < 0.05; Table 2), and previous studies (e.g., ^{Araújo et al. 2011}, ^{Keppel et al. 2015}, ^{Wan et al. 2018}) have shown that climate change has a high potential to threaten the effectiveness of protected areas in terms of conserving plant diversity at large scales. The exploration of the legacy effects of climate on the large-scale distributional patterns of plant richness is key for plant diversity conservation in protected areas in China. Based on the database of plant richness data from Chinese protected areas (with low spatial autocorrelation in plant richness data across different taxonomic levels; Figure 2), we examined the legacy effects of historical and contemporary climate on the large-scale distributional patterns of plant richness across different taxonomic levels.

The results of our GWR for both the historical and contemporary climate show significant explanatory power in terms of plant richness in China across various taxonomic levels (Table 4; Figures 3, 4), which hints that historical and contemporary climate affect the large-scale distributional patterns of plant richness across various taxonomic levels. Interestingly, the MAT was more important in the historical climate scenario than the present day scenario; however, the opposite was true of MAP according to the results of the explanatory power analysis (Figures 3, 4). The explanatory power of the historical and contemporary climate variables was higher at the family and genus levels and decreased at the species level based on the results regarding R² (Figures 3, 4). ^{O'Brien et al. (1998)} showed that the realized distributional limits of families and genera (unlike species) could be constrained by year-round or seasonally high ambient energy and by seasonally low ambient energy regardless of the water regime. Therefore, the historical and contemporary climate can better explain family and genus richness than species richness. Furthermore, historical and contemporary climate variables (including MAT and MAP) were found to better explain the richness of fern (non-seed) plants than that of seed plants (Figures 3, 4).

Some studies (e.g., ^{Wang et al. 2012}, ^{Svenning & Sandel 2013}, ^{Svenning et al. 2015}) have found that both the historical and contemporary climate can explain large-scale distributional patterns of plant richness. However, these effects may change across different taxonomic levels depending on the MAT and MAP (Figures 3, 4). There is ample evidence showing that Quaternary climatic change shaped the current patterns of plant richness and endemism across different regions of the world (e.g., North America, Europe, and Africa; ^{Svenning & Skov 2007}, ^{Normand et al. 2011}, ^{Svenning et al. 2015}, ^{Barnosky et al. 2016}, ^{Cotton et al. 2016}). The dynamics of plant richness may not follow the climatic equilibrium (^{Svenning & Sandel 2013}, ^{Svenning et al. 2015}). In addition, a shift to a new climatic equilibrium can cause time lags (^{Svenning & Sandel 2013}, ^{Svenning et al. 2015}). Plant species may experience a slower response to changes in temperature than to those of precipitation, and in many places, the migration of these species has shown a time lag in response to temperature changes (^{Svenning & Sandel 2013}, ^{Normand et al. 2011}, ^{Svenning et al. 2015}).

Plant community structure is strongly influenced by water under environmental change, and plant richness is strongly correlated with current water availability on a large scale (^{O'Brien 1998}, ^{O'Brien et al. 1998}, ^{Yang et al. 2011}). For example, the distributional pattern of plant richness is an important link to late Cenozoic precipitation trends, and a positive correlation between the mean annual rainfall and woody plant richness can be observed in southern Africa (^{O'Brien 1998}, ^{O'Brien et al. 1998}). Therefore, the response lags of plant richness to historical temperature and the effects of current precipitation on plant richness may drive the distributional pattern of plant richness at a large scale.

^{Blonder et al. (2018)} found that paleoclimate (i.e., MAT and MAP) is a better predictor of the spatial pattern of contemporary functional plant composition than contemporary climate predictors. The spatial pattern of contemporary functional plant composition is related to the distribution of plant richness at large scales (^{Petchey & Gaston 2002}, ^{Thompson et al. 2005}, ^{White et al. 2018}). Furthermore, plant species diversity may change more than functional-trait diversity because high levels of trait-based redundancy imply that the loss of a particular species should not affect ecosystem functions because of the maintenance of other species with similar traits (^{Dı́az & Cabido 2001}, ^{Petchey & Gaston 2006}). Relationships between functional-trait diversity and plant richness can still exist at various spatial and temporal scales (^{Petchey & Gaston 2002}, ²⁰⁰⁶, ^{Kraft et al. 2015}). Hence, historical MAT may have greater effects than contemporary MAT on the distributional pattern of plant richness, and contemporary MAP could also strongly affect the plant richness pattern.

The explanatory power of the historical and contemporary climate variables in regard to plant richness varied across the different taxonomic levels (Table 4; Figures 3, 4). Specifically, the historical and contemporary climate had the strongest explanatory power in terms of the family richness and the smallest explanatory power in regard to the species richness (Table 4; Figures 3, 4). The family taxonomic level can define the collective evolutionary distinctiveness of a set of species (^{Qian & Ricklefs 2007}, ^{Huang et al. 2016}). Furthermore, species’ abilities to thrive in an environment, resist and solve physiological problems, interact with other species and influence various ecosystem processes are determined by their functional traits (^{Dı́az & Cabido 2001}, ^{Fonseca & Ganade 2001}, ^{Rosenfeld 2002}, ^{Kuebbing et al. 2018}). In other words, while some species exhibit uncommon traits (functionally unique species), other species are functionally similar (i.e., represent redundant species) within one specific family (^{Naeem 1998}, ^{Rosenfeld 2002}). Hence, the effects of historical and contemporary climate on the large-scale distributional pattern of plant richness were found to be significant at the family level.

Furthermore, other recent environmental variables (e.g., human influences and usage of land) may influence the distributional patterns of plant species richness at large scales (^{Kier et al. 2005}, ^{Gerstner et al. 2014}). Such effects may escalate with increases in habitat areas and ranges (^{Lundholm 2009}). The habitat areas and distributional ranges of plants are generally wider at the family level than at the species level (^{O'Brien et al. 1998}, ^{Huang et al. 2016}). Therefore, the explanatory power of the historical and contemporary climate in regard to plant richness may rely on changes in taxonomic level (i.e., family, genus, and species) due to changes in habitat areas and distributional ranges.

Our results indicate that the influence of climate on plant richness at large scales differs between non-seed and seed plants (Table 4 and Figure 3). Fossil plant records (e.g., ^{Dubiel 1987}, ^{Collinson 2001}, ²⁰⁰², ^{Watkins & Cardelús 2012}, ^{Naugolnykh et al. 2016}) indicate that paleoclimates have affected the large-scale distributional pattern of ferns, while current bioclimatic variables, mainly those related to humidity (as water is an essential medium for fern reproduction), are closely associated with the variation in fern community composition. The physiological requirements and relative habitat restrictions of fern plants make them more sensitive to climate change than seed plants, and the effects of climate on plant richness may therefore differ between fern and seed plants (^{Schneider et al. 2004}, ^{Peppe et al. 2014}).

We found that the explanatory power of the historical and contemporary climate variables was higher for gymnosperm plant richness than for angiosperm richness (Table 4 and Figure 4). Such variation in explanatory power may be the result of evolutionary history and physiological adaptions to historical and contemporary climate (^{Wang et al. 2010}, ^{Yang et al. 2014}, ^{Lu et al. 2018}, ^{Xu et al. 2018}). For instance, ^{Lu et al. (2018)} found that herbaceous plants usually have higher molecular substitution rates than woody plants, partly due to their shorter generation times, apparently enabling herbaceous species in China to adapt quickly in response to climate change through increased genetic divergence and higher speciation rates.

Environmental heterogeneity and precipitation are the most important predictors of the diversity patterns of gymnosperms, followed by historical temperatures (^{Lü et al. 2018}). A number of gymnosperm plants are distributed mainly in western China (^{Lü et al. 2018}), and there is a large difference in historical temperature and contemporary precipitation between eastern and western China (^{Qin et al. 2015}, ^{Lu et al. 2018}, ^{Lü et al. 2018}). Furthermore, the temperature sensitivity of spring tree growth, water use, and successional strategies vary dramatically between the dominant angiosperm and gymnosperm plants (^{Bond 1989}, ^{Ma et al. 2016}, ^{Wan et al. 2017}, ^{Trugman et al. 2018}). The differences in the variables correlated with plant richness between gymnosperm and angiosperm plants may be related to their evolutionary histories and physiological adaptions to historical temperature and contemporary precipitation (^{Lu et al. 2018}, ^{Lü et al. 2018}). Hence, historical and contemporary climate influence the large-scale distributional patterns of gymnosperm and angiosperm plant richness to different extents.

Our study used protected area data to explore the effects of historical and contemporary climate on the large-scale distributional patterns of plant richness across different taxonomic levels, contributing to the conservation of plant diversity in China. China is the country/territory on Earth experiencing the greatest degree of land transformation (^{Liu et al. 2003}, ^{López-Pujol et al. 2006}, ^{Zhang et al. 2017}). Based on our results that the historical climate can shape the large-scale distributional pattern of plant richness, we predict that it will take a long time for the plant diversity in China to recover if the plant richness is damaged. Furthermore, it will be necessary to use large-scale data from protected areas to assess the effects of climate change on plant diversity in protected areas around the world (^{Araújo et al. 2011}, ^{Wan et al. 2014}, ^{Keppel et al. 2015}, ^{Wan et al. 2018}). Thus, plant diversity data from protected areas could be beneficial not only for scientists but also for decision makers and practitioners in other fields (^{Araújo et al. 2011}, ^{Wan et al. 2016}, ^{Zhang et al. 2017}). Hence, future studies should pay attention to the strategic value of different protected areas in the context of the plant diversity they protect, especially considering the current trends of land transformation around the world.

In conclusion, the distributional patterns of plant richness at large scales could be predicted across different taxonomic levels after the assessment of paleoclimate and contemporary climate data. Pleistocene temperature and current precipitation effects were studied to understand plant richness patterns, and such effects were found to be most important at the family level. In particular, the historical and contemporary climate data were better correlated with fern plant richness than with seed plant richness. Understanding the effects of historical and contemporary climate on the large-scale distributional patterns of plant richness across various taxonomic levels may help guide predictions of future plant diversity and facilitate the conservation of plant diversity under climate change.

Supplementary material

Data on plant richness, mean annual temperature (MAT; °C) and mean annual precipitation (MAP; mm). The data on plant richness refer to different taxonomic levels (i.e., family, genus, and species) based on the vascular, fern, seed, gymnosperm, and angiosperm plant groups.