Climatic Change Can Influence Species Diversity Patterns and Potential Habitats of Salicaceae Plants in China

Salicaceae is a family of temperate woody plants in the Northern Hemisphere that are highly valued, both ecologically and economically. China contains the highest species diversity of these plants. Despite their widespread human use, how the species diversity patterns of Salicaceae plants formed remains mostly unknown, and these may be significantly affected by global climate warming. Using past, present, and future environmental data and 2673 georeferenced specimen records, we first simulated the dynamic changes in suitable habitats and population structures of Salicaceae. Based on this, we next identified those areas at high risk of habitat loss and population declines under different climate change scenarios/years. We also mapped the patterns of species diversity by constructing niche models for 215 Salicaceae species, and assessed the driving factors affecting their current diversity patterns. The niche models showed Salicaceae family underwent extensive population expansion during the Last Inter Glacial period but retreated to lower latitudes during and since the period of the Last Glacial Maximum. Looking ahead, as climate warming intensifies, suitable habitats will shift to higher latitudes and those at lower latitudes will become less abundant. Finally, the western regions of China harbor the greatest endemism and species diversity of Salicaceae, which are significantly influenced by annual precipitation and mean temperature, ultraviolet-B (UV-B) radiation, and the anomaly of precipitation seasonality. From these results, we infer water–energy dynamic equilibrium and historical climate change are both the main factors likely regulating contemporary species diversity and distribution patterns. Nevertheless, this work also suggests that other, possibly interacting, factors (ambient energy, disturbance history, soil condition) influence the large-scale pattern of Salicaceae species diversity in China, making a simple explanation for it unlikely. Because Southwest China likely served as a refuge for Salicaceae species during the Last Glacial Maximum, it is a current hotspot for endemisms. Under predicted climate change, Salicaceae plants may well face higher risks to their persistence in southwest China, so efforts to support their in-situ conservation there are urgently needed.


Introduction
Understanding geographical variation in biological diversity and its underlying mechanisms has re-emerged as a research hotspot in ecology and biogeography [1][2][3].With accelerating specimen digitization and more high-resolution environmental data, biodiversity can now be gleaned from online databases (e.g., Global Biodiversity Information Facility, GBIF).Over the last 20 years, such data Forests 2019, 10, 220 3 of 22 so the simulation results were able to fully capture and convey this family's nationwide geographical distribution patterns [35,36].In summary, the above three considerations validate our use of ENM to simulate the potential distribution of Salicaceae in China at different times (past, present, and future).
Here, we took advantage of a comprehensive georeferenced occurrence dataset of Salicaceae in China coupled to high-resolution environmental data, to which we applied Maxent algorithms to build species distribution models.We hypothesized that dry and cold climatic conditions during the glacial period would have decreased the total suitable distribution area for Salicaceae plants and that with more climatic warming in the future, suitable habitats in the present would also face the risk of being lost.Based on these presumptions, we had three goals: (1) to project and quantify the Salicaceae species' changes in the extent of their potential suitable distribution under a variety of past, present, and future climate scenarios/years; (2) to map the spatial patterns of species diversity and endemism for Salicaceae members; (3) to determine which factors play key roles in shaping the patterns of Salicaceae plant species richness, weighted endemism, and corrected-weighted endemism.

Spatial Data
Species distribution data relating to the occurrence of Salicaceae plants were extracted from China National Specimen Information Infrastructure (NSII) (http://www.nsii.org.cn/), the Chinese Virtual Herbarium (CVH) (http://www.cvh.org.cn/), and the Global Biodiversity Information Facility (GBIF) (https://www.gbif.org/).All the species names were verified and rectified by using Taxonomic Name Resolution Service (TNRS) and by taxonomic experts.We obtained 10,322 Salicaceae specimens in total, of which 2107 were Populus and 8215 were Salix; the latter included specimens of 5 subgenera: 2848 of Chamaetia, 42 of Chosenia, 2120 of Salix, 61 of Triandra, and 3144 of Vetrix.However, because most specimens we collected lacked information on their latitude and longitude, we georeferenced these records based on their described location information.When collecting plant specimens, more intensive sampling is typically carried out in one area or region over others [37], and this sampling bias will skew the data representativeness for the species distributions vis-à-vis climate variables.To reduce and correct this spatial bias, we applied a systematic sampling method [38] using a resolution of 1 × 1 km, due to the spatial error in georeferenced specimen records and to match the resolution of environmental data (WorldClim).Our study region was divided into a grid of equally-sized square areas of 1 km 2 , which ensures that each grid cell has only one distributed record.After filtering, 2673 uniquely located records remained (Table S1), each with an accurate location, were used to simulate potential distributions of the Salicaceae family: 1025 Populus, 2306 Salix, 760 Chamaetia, 29 Chosenia, 1136 subgenus Salix, 43 Triandra, and 1114 Vetrix records.To obtain the species richness pattern for Salicaceae, we simulated the potential distribution of each species.Species currently occupying less than 5 grid cells were removed from our analyses [39,40].A total of 6678 georeferenced specimen records belonging to the 5 Salix subgenera representing 215 species (Table S2) were used for these niche models.To avoid too much data stifling the model runs, the simulation results were resampled with a spatial resolution of 20 × 20 km for analyzing the relationships between species diversity patterns and environmental variables.

Environmental Parameters
The growth and distribution of Salicaceae species are susceptible to various interacting factors, namely climate (e.g., temperature and precipitation), UV-B radiation, and soil [41][42][43].For example, UV-B radiation and temperature caused gender imbalance in Salicaceae plants, driving divergent evolutionary trends in Populus and Salix genera and changes in their population numbers [44,45].Hageer et al. [46] and Chen et al. [47] studies showed that in addition to climate variables, the effects of soil variables should also be taken into account in the simulation of the potential distribution of species.Hence, for our modeling, we chose four types of environmental variables: climate, UV-B radiation, soil, and topography.
Initially, 39 environmental variables were selected to build a species distribution model.These consisted of elevation, 19 bioclimatic variables, 15 Food and Agriculture Organization (FAO) soil variables, and 4 UV-B variables.Before any modeling, to avoid multicollinearity of variables-this leads to over-fitted species distribution models, with the highly correlated variables affecting model accuracy-we calculated Spearman's rank correlation coefficients among them, selecting those variables with r s -values < 0.75 for inclusion (Tables S3 and S4).This key filtering step resulted in 15 environmental variables in total used to predict the patterns of Salicaceae (5 bioclimatic, 1 topographic, 5 soil, and 4 UV-B factors; as listed in Table 1).Data on the 5 bioclimatic factors and elevation were obtained from the WorldClim website (https://www.worldclim.org/),while those for the 5 soil factors were obtained from the China Soil Map Based Harmonized World Soil Database (v1.1) (http://westdc.westgis.ac.cn/).Data on the 4 UV-B factors were taken from the gIUV database (http://www.ufz.de/gluv/).All data were resampled to a spatial resolution of 1 km 2 .

Historical and Future Climate Scenarios
The Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC) shows that, over the past 100 years (1880-2012), the average annual temperature at the Earth's surface has increased by 0.85 • C, most of which (0.072 • C rise) has occurred in the last 60 years (1951-2012) [24].The IPCC's AR5 adopted a climate-coupled model of the international coupling model for phase 5 (CCMIP5) and discussed four representative greenhouse gas emission scenarios: RCP2.6,RCP4.5, RCP6.0, and RCP8.5.RCP stands for "representative concentration pathway" (i.e., the typical concentration trajectory); its assigned number conveys the amount of radiative forcing (RF, in W/m −2 ) in 2100 relative to its value in 1750.We selected the two most typical concentration trajectories, RCP4.5 and RCP8.5, since they respectively represent an intermediate stable trajectory and a highly concentrated one.
To better capture the changes in historical climate, we calculated the difference between temperature and precipitation during the present and LGM (i.e., bio1-anomaly = Bio1 present − Bio1 LGM ; Table 3).

Building the Species Distribution Model (SDM)
The species distribution model (SDM) is widely used to simulate the suitable distribution pattern of species under climate change scenarios [48][49][50].The maximum entropy (Maxent) model is a type of machine-learning algorithm model.Specifically, its algorithm adopts the Maxent principle, of not constraining any unknown distribution information yet also preserving the information constraint of the distributed environment variable data, to predict scientifically the potential distribution of one or more species.A logistically attractive feature of Maxent is that it can predict species' potential distributions by their presence-only data and environmental variables.Compared with other SDMs, the predictions of Maxent are better when the number of distribution points is uncertain and the correlation between different climatic factors is unclear [51].In recent years, SDMs using the maximum entropy model have played a major contributing role to species diversity conservation in the future [52].
Using the environmental variables at present associated with the distribution points of Salicaceae species, we constructed the model based on the maximum entropy theory to simulate this family's realized niche.Next, these simulated real ecological niches were projected to different times (i.e., past or future) for calculation to build the model [53].After modeling the current suitable habitat area for Salicaceae with current climate data, we conducted modeling projections for past (i.e., LIG and LGM) and future climatic scenarios (i.e., RCP4.5-2050,RCP4.5-2070,RCP8.5-2050, and RCP8.5-2070) to predict their respective availability of suitable habitat.These SDM projections were made in Maxent software.Based on the 15 environmental variables and our species distribution data, we built an SDM model for the Salicaceae family for past, present, and future climate scenarios.For 215 species, we also simulated the potential distribution of each by using the maximum entropy model in MaxEnt v3.3 software (http://biodiversityinformatics.amnh.org/open_source/maxent/).Only those species with at least 5 distribution records were used to establish the species distribution model, as anything less is considered insufficient for model predictive ability [40].MaxEnt was run with these rules: use only linear features for SDM with distribution data of <5 records; add quadratic features when occurrence records consisted of ≥5 and <15 records; add hinge features for ≥15 but <80 records; include product and threshold features when species records exceeded 80 [54].In China, Salicaceae species have rather specific environmental niches and geographical distribution ranges; hence, we applied the clamping function of Maxent.We calculated the contribution of each environmental variable via Jackknife testing and the set cross-validation function, with logistic outputs.This operation was repeated 10 times, with 75% of species occurrence records used for model training and the remainder (25%) for model testing.To reduce model over-fitting to low levels, we set the regularized multiplier to a value of 2.00 and the number of background points to 10,000 [55].Next, an ASCII grid layer was produced with the largest obtained AUC (area under the receiver operating characteristic (ROC) curve).The ROC curve was employed to test the model's predictive accuracy by judging the AUC value (range of 0-1): a prediction was perfect if its AUC value = 1.This rarely happens, yet when AUC > 0.7 the simulation results are considered valid [56,57].
Next, the arithmetic results from MaxEnt were loaded into ArcGIS v10.2 to carry out a fitness classification and visualization expression, which enabled us to generate the potential distributions of Salicaceae plants.The source of the map was the National Catalogue Service for Geographic Information (http://www.webmap.cn/main.do?method=index).It was critical to choose an appropriate probability threshold of species existence when converting the continuous suitability index maps to binary habitats.Extensive research has demonstrated that a "maximum training sensitivity-plus-specificity" approach leads to highly accurate predictions and is thus superior to other threshold division methods [58,59].Therefore, to demarcate the habitat and non-habitat of Salicaceae species, we used a "maximum training sensitivity-plus-specificity" threshold in our study.

Biodiversity Pattern Indices
We used three common biodiversity metrics-species richness, weighted endemism, and corrected-weighted endemism-to quantify the species diversity patterns.Once the probability threshold of species existence had been determined, the matrix of all presence/absence layers was generated.Its rows represented the 26,137 grid cells covering China's landmass while its columns indicated the presence of the 215 Salicaceae species we modeled.
Species richness (SR) was defined as the sum of unique species occurring per cell.It is one of the most commonly used indice in large-scale biodiversity conservation research, characterized by its simplicity and intuition.Endemism is the limitation of a biological group unit (species, genus, or family) to a particular geographical area.In this respect, it is of great value to understand the nature, occurrence, and evolution of a regional flora.Weighted endemism (WE), the sum of the reciprocal total number of cells each species in a grid cell is found in, emphasizes areas having a high proportion of species with restricted ranges.Corrected-weighted endemism (CWE) is simply the weighted endemism divided by the total number of species in a cell [60].Hence, the CWE emphasizes areas that have a high proportion of species with restricted ranges, yet these are not necessarily species-rich areas.Both the WE and CWE can be used to reveal the distribution center of endemism phenomena.The spatial distribution pattern of the WE value will depends largely on the species richness pattern.The CWE fully embodies the spatial distribution characteristics of endemic species, but it lacks objectivity because it uses thean artificial boundary.Therefore, these three indices complement each other.SR, WE, and CWE were mapped in ArcGIS v10.2, and calculated as follows: where C is the number of grid cells in which each endemic species occurs, and K is the total number of species in a grid cell.

Changes in Core Distribution Centers
To examine changing trends from the past to the future, we tested and compared the changes in potential habitat area and distribution centers of the past, present, and future suitable areas by using the SDM toolbox [61].This tool calculates the distributional changes between two binary SDMs (i.e., past vs.present SDMs), providing a table output depicting the predicted contraction, expansion, and areas of no change in a given species distribution [62].The tool was also used to reduce Salicaceae' distributions to a centroid and to form a vector file to describe the changes in a centroid's range and direction.

Statistical Analysis
To explain the species diversity patterns of Salicaceae, we conducted variation partitioning analysis.It is useful to explore the relative explanatory power of independent variables, especially between two groups of environmental factors-i.e., contemporary environments vs. historical climate change-in forming species diversity patterns.This explanatory power from environmental variables amounted to 60.9% (Figure S2), so their selection should help to better explain biodiversity patterns.Generalized linear models (GLM) were used to explore the relationships between the species diversity indices and environmental variables.We built three models: (i) species richness and environmental variables model; (ii) weighted endemism and environmental variable model; and (iii) corrected-weighted endemism model.However, spatial autocorrelation in species diversity patterns might inflate Type I error rates [63].To eliminate this influence on our significance testing, we used modified the t test applied to all the models' coefficients [64,65].The geographical coordinate system used for the central point of each grid cell was WGS1984.All statistical analyses were done in R v3.5.1 (http://www.r-project.org/).

Species Distribution Model and Its Accuracy
China was divided into 26,137 grid cells, with 2673 of them (10.22%)having exact geographic coordinates.Models for the Salicaceae family with a cross-validation AUC close to 0.75 were considered useful in our study.The training and testing of the AUC value for the 215 species models was >0.8, which indicated these models were non-random and could be reliably used for modeling the relationship between diversity patterns and environmental factors.Potential suitable habitat areas are currently widely distributed in China's subtropical and warm temperate zones, but occurred mainly in its central and southwestern regions.Highly suitable areas were found in the Hengduan Mountains, Daba Mountains, and Qinling Mountains (Figure 1).corrected-weighted endemism model.However, spatial autocorrelation in species diversity patterns might inflate Type I error rates [63].To eliminate this influence on our significance testing, we used modified the t test applied to all the models' coefficients [64,65].The geographical coordinate system used for the central point of each grid cell was WGS1984.All statistical analyses were done in R v3.5.1 (http://www.r-project.org/).

Species Distribution Model and its Accuracy
China was divided into 26,137 grid cells, with 2673 of them (10.22%)having exact geographic coordinates.Models for the Salicaceae family with a cross-validation AUC close to 0.75 were considered useful in our study.The training and testing of the AUC value for the 215 species models was >0.8, which indicated these models were non-random and could be reliably used for modeling the relationship between diversity patterns and environmental factors.Potential suitable habitat areas are currently widely distributed in China's subtropical and warm temperate zones, but occurred mainly in its central and southwestern regions.Highly suitable areas were found in the Hengduan Mountains, Daba Mountains, and Qinling Mountains (Figure 1).

Changes in the Potential Range of Salicaceae
Compared with the LIG, the niche model predicted losses in current suitable habitat areas mainly distributed in the middle temperate zone of northern China, namely Inner Mongolia, Heilongjiang, Jilin, Xinjiang, Qinghai, and Gansu (Figure 2A).The projected area of reduction amounted to 139.72 × 10 4 km 2 , which accounts for 38.67% of the currently suitable area (Table 2).This range contraction was mainly concentrated in high latitude areas.Conversely, the gains in suitable range area totaled 62.79 × 10 4 km 2 (17.38% of the currently suitable area), mostly distributed in provinces of Ningxia, Shaanxi, Gansu, Henan, and Hebei.Expansion occurred mainly in the midlatitudes.Overall, from the last glacial period to the present, we saw a wide range of contractions, with the total suitable habitats reduced by ca.76.92 × 10 4 km 2 (Figure 2A, Table 2).In contrast to the

Changes in the Potential Range of Salicaceae
Compared with the LIG, the niche model predicted losses in current suitable habitat areas mainly distributed in the middle temperate zone of northern China, namely Inner Mongolia, Heilongjiang, Jilin, Xinjiang, Qinghai, and Gansu (Figure 2A).The projected area of reduction amounted to 139.72 × 10 4 km 2 , which accounts for 38.67% of the currently suitable area (Table 2).This range contraction was mainly concentrated in high latitude areas.Conversely, the gains in suitable range area totaled 62.79 × 10 4 km 2 (17.38% of the currently suitable area), mostly distributed in provinces of Ningxia, Shaanxi, Gansu, Henan, and Hebei.Expansion occurred mainly in the mid-latitudes.Overall, from the last glacial period to the present, we saw a wide range of contractions, with the total suitable habitats reduced by ca.76.92 × 10 4 km 2 (Figure 2A, Table 2).In contrast to the LIG, the suitable habitats increased significantly-the total suitable habitat area increased by ca.69.63 × 10 4 km 2 -from  2).These gains amounted to 106.88 × 10 4 km 2 , or 29.58% of the currently suitable area, and were concentrated in parts of Tibet, Sichuan, Qinghai, Gansu, Ningxia, Inner Mongolia, Shannxi, Hebei, Liaoning, Beijing, Tianjin, Jilin, and Heilongjiang.Expansion in habitats increased with increasing latitude.The predicted losses occurred mainly in the subtropics of southern China, namely Guangxi, Guangdong, and Fujian, covering an area of ca.37.25 × 10 4 km 2 (10.31% of the currently suitable area; Figure 2B, Table 2).LIG, the suitable habitats increased significantly-the total suitable habitat area increased by ca.69.63 × 10 4 km 2 -from the LGM to the present (Figure 2B, Table 2).These gains amounted to 106.88 × 10 4 km 2 , or 29.58% of the currently suitable area, and were concentrated in parts of Tibet, Sichuan, Qinghai, Gansu, Ningxia, Inner Mongolia, Shannxi, Hebe, Liaoning, Beijing, Tianjin, Jilin, and Heilongjiang.Expansion in habitats increased with increasing latitude.The predicted losses occurred mainly in the subtropics of southern China, namely Guangxi, Guangdong, and Fujian, covering an area of ca.37.25 × 10 4 km 2 (10.31% of the currently suitable area; Figure 2B, Table 2).Under all four future climate scenario/year combinations of RCP4.5-2050,RCP4.5-2070,RCP8.5-2050, and RCP8.5-2070, the predicted area of suitable habitats shrunk (Table 2).The parts of Salicaceae's range with high risks of habitat loss were distributed in subtropical zones, including southern and eastern regions of China, Jiangxi, Fujian, eastern Hunan, southwestern Yunnan, northern Guangxi and Guangdong, southeast Hubei, the border of Chongqing and Sichuan, the junction of Anhui, Henan, Jiangsu, and Shandong, as well in limited regions of Hebei, Tianjin, and Beijing (Figure 2C-F).The Maxent model predicted total losses of 20.21 × 10 4 km 2 ~82.06 × 10 4 km 2 (5.59%~22.71% of the currently suitable area; Table 2).These losses were mostly concentrated at low elevation and the contraction generally increased at lower latitudes.Most of the predicted gains in suitable habitat area were restricted to Salicaceae's range in the north of China's warm temperate zone (Tibet, Sichuan, Qinghai, Inner Mongolia, Heilongjiang, and Jilin), which amounted to ca. 20.18 × 10 4 km 2 ~78.06 × 10 4 km 2 (5.58%~21.60% of the currently suitable area; Table 2), all located in high latitude regions.The total potential suitable habitat area decreased under the four future climatic scenario/years, with contractions in Jiangxi, Fujian, and eastern Hunan most apparent (Figure 2C-F).

Core Distributional Shifts
The center of predicted suitable areas in the present was located in the Hubei province, with geographical coordinates of 110.62 E and 32.79 N (Figure 3).The centroid of potential suitable habitats in the LIG shifted to northern regions (111.11 E, 35.12 N) and high latitudes.By contrast, the centroid in the LGM period shifted to southern regions and low latitude (110.28E, 29.52 N).Under the future climate scenario/years combinations, the center of potential suitable habitats was located at more northerly areas with higher latitudes, except under RCP8.5-2070.Evidently, the distribution center of Salicaceae plants is likely to shift to high latitudes from the LGM into the future (Figure 3).geographical coordinates of 110.62 E and 32.79 N (Figure 3).The centroid of potential suitable habitats in the LIG shifted to northern regions (111.11 E, 35.12 N) and high latitudes.By contrast, the centroid in the LGM period shifted to southern regions and low latitude (110.28E, 29.52 N).Under the future climate scenario/years combinations, the center of potential suitable habitats was located at more northerly areas with higher latitudes, except under RCP8.5-2070.Evidently, the distribution center of Salicaceae plants is likely to shift to high latitudes from the LGM into the future (Figure 3).

Relationships between Species Richness and Environmental Parameters
The pattern of Salicaceae species richness was obtained from 215 species potential distribution model layers.Botanical species richness ranged from 0 to 115 in each grid cell in China, with an average of 35 per cell.The southwest and central regions of China showed the highest Salicaceae species diversity, which included the eastern Tibet Plateau and the Hengduan Mountains, Qinling Mountains, Daba Mountains, Dalou Mountains, Wumeng Mountains, Wulian Mountains, Longmen Mountains, Daxue Mountains, Liupan Mountains, and Kunji Mountains.In short, high species diversity was mainly concentrated in alpine areas.The northeast region of China ranked second in species diversity.By contrast, species diversity of northwest China, south China, and the western Tibet Plateau were relatively low (Figure 4).periods), current and future (four scenarios) suitable areas, with the arrows depicting the magnitude and direction of predicted change through time.The map was plotted using ArcGIS v10.2.

Relationships between Species Richness and Environmental Parameters
The pattern of Salicaceae species richness was obtained from 215 species potential distribution model layers.Botanical species richness ranged from 0 to 115 in each grid cell in China, with an average of 35 per cell.The southwest and central regions of China showed the highest Salicaceae species diversity, which included the eastern Tibet Plateau and the Hengduan Mountains, Qinling Mountains, Daba Mountains, Dalou Mountains, Wumeng Mountains, Wulian Mountains, Longmen Mountains, Daxue Mountains, Liupan Mountains, and Kunji Mountains.In short, high species diversity was mainly concentrated in alpine areas.The northeast region of China ranked second in species diversity.By contrast, species diversity of northwest China, south China, and the western Tibet Plateau were relatively low (Figure 4).The GLMs showed that the variables related to contemporary energy had stronger explanatory ability than contemporary water availability, soil conditions, heterogeneity, and historical climate change.Specifically, UV-B radiation was the strongest variable, as it explained 35% (p < 0.05).Species richness decreased as UV-B factors and annual mean temperature increased, indicating that areas with high species diversity are warm with low UV-B radiation.Elevation and anomaly of precipitation seasonality both adversely affected species diversity patterning, while annual precipitation and nitrogen content of the topsoil were positively associated with it.Those places featuring a climate warmer and wetter than in the LGM tended to have more Salicaceae species than other regions (Table 3).The GLMs showed that the variables related to contemporary energy had stronger explanatory ability than contemporary water availability, soil conditions, heterogeneity, and historical climate change.Specifically, UV-B radiation was the strongest variable, as it explained 35% (p < 0.05).Species richness decreased as UV-B factors and annual mean temperature increased, indicating that areas with high species diversity are warm with low UV-B radiation.Elevation and anomaly of precipitation seasonality both adversely affected species diversity patterning, while annual precipitation and nitrogen content of the topsoil were positively associated with it.Those places featuring a climate warmer and wetter than in the LGM tended to have more Salicaceae species than other regions (Table 3).  1 for the abbreviations for environmental variables.SR: species richness, WE: weighted endemism, CWE: Corrected-weighted endemism, Bio1-Ano: Bio1 anomaly, Bio3-Ano: Bio3 anomaly, Bio7-Ano: Bio7 anomaly, Bio12-Ano: Bio12 anomaly, Bio15-Ano: Bio15 anomaly.

Relationships between Weighted Endemism (WE), Corrected-Weighted Endemism (CWE), and Environmental Factors
Plant-weighted endemism ranged from 0 to 0.0295 in each grid cell in China, with an average value of 0.0068, and 1.61% of the total grids had more values >0.0200 (Figure 5A).Weighted-endemism followed a similar pattern as did that found for species richness.The Himalayas, and the mountains of Hengduan, Qinling, Daba, Longmen, Wumeng, Wulian, Yunling, and Kunji all featured high values of weighted endemism.The corrected-weighted endemism ranged from 0 to 0.0004 (Figure 5B), with the high values concentrated in northwest China, including the Himalayas and Hengduan Mountains.The latter is where endemic species of Salicaceae had their ranges most narrowly distributed.
The GLMs showed that both historical climate change and contemporary water availability were the main factors shaping the patterns of weighted endemism (Table 3).The anomalies of precipitation seasonality and annual precipitation explained 22% of its variation.While nitrogen content of topsoil and UV-B radiation also had significant positive effects on the corrected-weighted endemism pattern.Elevation was a significant factor in the three types of pattern, but it showed weaker explanatory power for the weighted endemism and corrected-weighted endemism patterns.The relationships between nitrogen content of topsoil and patterns of weighted endemism and corrected-weighted endemism were significant, showing strong positive correlations (Table 3).The GLMs showed that both historical climate change and contemporary water availability were the main factors shaping the patterns of weighted endemism (Table 3).The anomalies of precipitation

Changes in the Potential Range of Salicaceae in China
Using the georeferenced occurrence data of Salicaceae, our results showed that its predicted potential distribution area covers most regions of China.The predicted results are consistent with the actual distribution of Salicaceae plants in China, which features strong regularity and forms part of the temperate deciduous forest community.The similarity of all the plant species-they are deciduous, water demanding, and require a temperate bioclimate, to name a few key traits-and the model's effective prediction accuracy (>0.7) support the generalization of our results.Second, by simulating Salicaceae species' ranges and patterns in different periods we could infer their historical and future population dynamics (Figure S1).Salicaceae may have experienced a significant and extensive expansion from south to north during the Last Inter Glacial period (LIG).This strongly suggests Salicaceae can track suitable niches during the glacial and interglacial cycles, in that its distribution shifted southward in the glacial period but shifted back northward after the glacial period.Salicaceae is a group of plants whose ecological amplitude is large and they can adapt to many different ecological environments, from those ranging from high (~20 • C) to low (~3 • C) temperatures, and from swamp to sand habitats [66,67].The temperature (average: ~5 • C) and rainfall (average: ~557 mm) in the LIG were within the tolerance levels of Salicaceae species.However, with significant reductions in temperature and rainfall, the climate became drier and colder in the Last Glacial Maximum (LGM) [68,69].In response, because those areas with suitable habitat were generally greatly reduced, Salicaceae's distribution range retreated to areas low in elevation [70].Therefore, from these results we infer this plant family was able to spread rapidly in the glacial-interglacial cycle [11].Most of its species (i.e., members of Salix subgenera Chamaetia and Vetrix, represent ca.75% species of Salicaceae) are Arctic-alpine taxa well adapted to cold and hostile environments [71].Therefore, the considerable northward expansion in China of the Salicaceae distribution area is both plausible and reasonable.In sum, Quaternary climate change profoundly influenced the historical spatial population dynamics of ancient Chinese Salicaceae species [72].
As the climate warms, extreme weather events are becoming more frequent and intense [73][74][75].Climate research has made much progress in recent years, but most of it is based on qualitative rather than quantitative analyses [76][77][78][79].In our study, we conducted a detailed quantitative analysis of the potential distribution of Salicaceae in the future.Nevertheless, simulation results are expected to vary between different climate models [80,81].Salicaceae are widely distributed, occurring in Africa, Europe, Asia, South America, and North America from ca. 82 • N to 52 • S [82], but their habitats are often fragmented.East Asia is the present distribution center of modern species of Salicaceae, and China harbors the most species numerically.
Our modeling study predicted ca.361.33 × 10 4 km 2 for Salicaceae species under the present climate scenario with its core distributional area being Hubei Province, which is broadly within the distribution area of extant Salicaceae.In our study's projections, the centroid of suitable distribution shifted northerly to higher latitude areas under the future climate scenarios, and this trend agrees with other works [83][84][85][86].Our results also show that if global warming becomes more intense, the suitable range of Salicaceae faces the risk of reduction under the medium and high carbon dioxide emission scenarios (RCP4.5 and RCP8.5), a result that is also in line with previous studies [79,87].Species of Salicaceae are, therefore, clearly vulnerable to climate change effects, with its range contractions mainly concentrated at low latitudes [88][89][90].Nevertheless, since the model involved all species of the family these results are not applicable to endangered species or any other particular species.Earlier studies by Bomhard et al. [91] and Midgley et al. [92] suggested that future climate change and land-use transformations would have synergistic effects on species habitats, giving rise to unsuitable habitats for their persistence.In this context, we found that Salicaceae species do display strong adaptability under the present and future climate scenarios; however, some of the suitable range may become unsuitable due to intense disturbance by human activities and land-use changes.For example, much of the suitable range for Salicaceae plants in China has been converted into farmland or construction land [93].Clearly, the effects from future land-use changes should also be taken into account in species distribution models [94].We suggest that local protection measures (Figure S3) can be adopted to reduce anthropocentric-driven damage to suitable habitat.Meanwhile, we can also try to implement ex situ protection projects, such as transplanting species from a threatened area to a site with non-or less-degraded habitat, or by planting them in multiple botanical gardens.

Species Richness and Endemism Patterns of Salicaceae in China
Using the distributions of 215 Salicaceae species, we explored their large-scale variation in species richness and endemism and the environment variables likely determining these patterns.China is famous for its high species diversity; better documenting and understanding of its diversity and endemic patterns will inform and strengthen biological conservation, by providing baselines for biodiversity management and relevant policy-making.China's southwest, lying south of the Hengduan Mountains, is well known for its high suitability for and diversity of Salicaceae species (Figures 2,  4 and 5).Our variance partition analyses revealed pure effects of contemporary environments and historical climate change that respectively accounted for 47.90% and 18.10% of the variance in Salicaceae diversity patterns across China (Figure S2).Consistent with other research [14,95], our study showed the climate exerts a important influence on Salicaceae species richness, weighted endemism, and corrected-weighted endemism (Table 3 and Figure 2 and Figure S2).For all three, the most important factor was the contemporary environment (energy, water availability and soil conditions), since it had the strongest explanatory power (Table 3).This indicates high species diversity of Salicaceae mostly occurs in areas with relatively warm and high humidity [96].High temperatures at low latitudes may curtail breeding times and accelerate species formation, favoring the emergence of new plant species [97].This is consistent with evidence that imposed warming considerably increased the root, stem, and leaf biomass of Salicaceae plants [98] with drought strongly impeding their growth [99,100].Meanwhile, considering that Salicaceae species originated in the temperate zone [101] and most of them prefer wet environments [102], the positive association between precipitation, temperature and species diversity lends supports to the water-energy dynamic equilibrium hypothesis as a species richness formation mechanism.It has been shown that water-energy dynamics (precipitation and coldest monthly evapotranspiration) were the main factors affecting the diversity of South African trees [103][104][105].Other studies have found biodiversity patterns that were also related to water availability [106,107].Temperature and liquid water are not only vital for plants to absorb and transport nutrients but also mediate key physiological activities (e.g., photosynthesis) [108,109].
Besides water and temperature, our results suggested UV-B radiation also contributes to the biodiversity pattern of Salicaceae plants.The negative associations indicate that where UV-B radiation is intense, both species richness and endemism are lower.UV-B radiation is a factor known to limit the distribution range of terrestrial and marine life: its negative relationship with biodiversity patterns arises from its adverse effects on the protective mechanisms and organs of plants [110].Related research has also found that high-intensity UV-B radiation can reduce the photosynthetic performance of Salicaceae species [43,111,112], in addition to diminishing the biomass accumulation [16], seed germination [113] and membrane structure integrity [44] of plants.Therefore, we propose that, generally, intensive UV-B radiation reduces Salicaceae plant biodiversity.Compared with other environmental variables, the association between soil fertility and biodiversity pattern was relatively weak, which is in line with other research [114].Nevertheless, it worth emphasizing that topsoil nitrogen content makes a non-trivial contribution to the diversity patterns of Salicaceae species.This may be linked to how N can markedly influence gender differences and competitiveness between the two sexes in these plants [47].Some studies do show that temperature stability [115] and historical climatic stability [116] also affect species diversity patterns.We found historical geological processes also played an essential role by altering the range (Figure 2A,B), diversity and endemism (Table 3, Figure S2) of Salicaceae species during the glacial period.Thus, while contemporary water-energy and historical climate change both have significant impacts on Salicaceae, other factors (such as disturbance history) warrant quantitative study in future research [117].It is well known that the ozone layer protects the forests and other life forms on the Earth's surface from the sun's harmful ultraviolet radiation [118].During the last few decades, the stratospheric ozone layer has been depleted by anthropogenic pollutants such as chlorofluorocarbons (CFCS), leading to an increase in surface-level ultraviolet radiation.As a result of the 1987 Montreal protocol and its amendments, the anthropogenic loading of ozone-depleting substances in the atmosphere is being reduced.Nevertheless, future climate change could cause extreme dynamic events [119], volcanic eruptions [120] and irregular changes in solar flux [121], all of which are likely to damage the ozone layer.Caldwell's results showed that for every 1% reduction in ozone, the surface's UV-B (280-320 mm) radiation increases by 2% [122].With greater UV-B radiation possible in the future, the geographical distribution and species diversity of Salicaceae plants may become threatened.

Conclusions
Based on niche model, we propose the Salicaceae family underwent an extensive population expansion of its species during the LIG, but retreated to low latitudes during the LGM.As climate warming intensifies, suitable habitats will shift to higher latitudes while those at lower latitudes will contract in area and abundance.Habitats of Salicaceae in China are most likely to be lost in the future, especially in Jiangxi, Fujian, and eastern Hunan.The western and central regions of China currently show the highest diversity of Salicaceae species.Because they are typical pioneer plants, with a strong ability to disperse and colonize, they could fill new or open habitats easily; these traits make it difficult to resort to a single hypothesis to reveal their pattern of species diversity.Nevertheless, our results point to the water-energy dynamic equilibrium and historical geological processes as being key drivers of Salicaceae diversity patterns in China.Our study could thus assist in the conservation and sustainable use of Salicaceae, and our results provide a solid baseline for biodiversity management and relevant policy-making for these important woody plants.

Figure 1 .
Figure 1.Present potential distribution pattern of Salicaceae plants in China, using the Albers projection.The map was plotted in ArcGIS v10.2.

Figure 1 .
Figure 1.Present potential distribution pattern of Salicaceae plants in China, using the Albers projection.The map was plotted in ArcGIS v10.2.

Figure 2 .
Figure 2. Spatial changes of Salicaceae plants in China under different combinations of climate scenario/years (Albers projection).Comparisons between the species distribution model (SDM) for Salicaceae in the present and (A) the SDM under the Last Inter Glacial ("lig"); (B) the SDM under the Last Glacial Maximum ("lgm"); (C) the SDM under future climate scenario RCP4.5 in 2050; (D) the SDM under future climate scenario RCP4.5 in 2070; (E) the SDM under future climate scenario RCP8.5 in 2050; (F) the SDM under future climate scenario RCP8.5 in 2070.The map was plotted in ArcGIS v10.2.

Table 1 .
Environmental variables used to predict the geographical distributions of Salicaceae plants.

Table 2 .
Dynamic changes in the suitable area for Salicaceae plants under different combination of climate scenario/years.Negative values indicate contractions in suitable habitat area.Comparisons between the suitable area for Salicaceae plants in the present and a the suitable area under the Last Inter Glacial ("LIG"); b the suitable area under the Last Glacial Maximum ("LGM"); c the suitable area under future climate scenario RCP4.5 in 2050; d the suitable area under future climate scenario RCP4.5 in 2070; e the suitable area under future climate scenario RCP8.5 in 2050; f the suitable area under future climate scenario RCP8.5 in 2070.

Table 3 .
Explanatory power (R 2 ) of the environmental variables for Salicaceae species diversity patterns evaluated by the generalized linear models (GLMs).Modified t-tests were used to examine the significance.