Prediction of the Suitable Area of the Chinese White Pines ( Pinus subsect. Strobus ) under Climate Changes and Implications for Their Conservation

: White pines ( Pinus subsect. Strobus ) play important roles in forest ecosystems in the Northern Hemisphere. Species of this group are narrowly distributed or endangered in China. In this study, we used a species distribution model (SDM) to project and predict the distribution patterns of the 12 species of Chinese white pine under a variety of paleoclimatic and future climate change scenarios based on 39 high-resolution environmental variables and 1459 distribution records. We also computed the centroid shift, range expansion / contraction, and suitability change of the current distribution area to assess the potential risk to each species in the future. The modeling results revealed that the suitable habitat of each species is consistent with but slightly larger than its actual distribution range and that temperature, precipitation, and UV radiation are important determining factors for the distribution of di ﬀ erent white pine species. The results indicate that the Last Glacial Maximum (LGM) greatly a ﬀ ected the current distribution of the Chinese white pine species. Additionally, it was predicted that under the future climate change scenarios, there will be a reduction in the area of habitats suitable for P. armandii , P. morrisonicola , and P. mastersiana . Furthermore, some of the current distribution sites of P. armandii , P. kwangtungensis , P. mastersiana , P. morrisonicola , P. sibirica , and P. wallichiana were predicted to become more unsuitable under these scenarios. These results indicate that some Chinese white pine species, such as P. armandii , P. morrisonicola , and P. mastersiana , may have a very high risk of population shrinkage in the future. Overall, this study provided relevant data for the long-term conservation (both in situ and ex situ) and sustainable management of Chinese white pine species. distribution model to dynamically project and quantify their suitable habitats under paleoclimatic, current, and future climate scenarios. We also explored the centroid shift, range expansion/contraction, and suitability change of the current distribution area to assess the potential risk to each species in the future. Our aims are to find out how climate fluctuation since the last


Introduction
There are 17 so-called "megadiverse" countries, predominantly located in tropical or subtropical regions, which harbor the majority of earth's species and have a high proportion of endemism [1]. China is one of the few megadiverse countries in the temperate zone of the Northern Hemisphere [2]. Therefore, the conservation of China's biodiversity will make a great contribution to the natural heritage of temperate zones [3]. However, biodiversity loss is continuing, and some researchers consider the earth to be undergoing the sixth mass extinction in its history [4][5][6]. The distribution of organisms has been greatly influenced by climate changes and habitat fragmentation, which bring about severe threats to biodiversity [3,7]. In the last few decades, rising temperatures have shifted the latitudinal and altitudinal ranges of many species, especially those distributed in alpine areas [7,8]. Thus, determining the effects of climate changes on the distribution of organisms is fundamental for understanding the conservation of biodiversity [9][10][11]. Due to its importance in ecosystems, P. subsect. Strobus has attracted much attention in recent years. Previous studies have studied white pines in terms of their phylogeny, ecology, population structures of endangered species, etc. [30,[35][36][37][38][39][40][41][42][43][44][45]. However, it is still unclear what are the key ecological factors determining the current distribution of P. subsect. Strobus, how the distribution pattern reacted to changes in paleoclimate, and how the distribution pattern will react under future climate changes. This knowledge will help us to identify response mechanisms, which will contribute to the development of effective conservation strategies in the future.
The maximum entropy method (MaxEnt) is a species distribution model (SDM) that is used to simulate the habitats that are suitable for species based on presence data and information of environmental variables [46]. Due to the model's effectiveness, it has been widely used in biodiversity (especially for endangered species) conservation [47][48][49][50][51], invasive species control [52,53], as well as nature reserve planning and management [54][55][56]. It has also been applied to reconstruct habitat shifts of individual species or vegetation during the last glacial period, and simulate how paleoclimatic fluctuations have influenced the present distribution [57][58][59], predict the potential distribution in the future, find conservation gaps, and guide protection management [60][61][62].
Niche conservatism is a key assumption underlying the application of SDMs [63]. According to a phylogenetic study of P. subsect. Strobus species, the current species diversity and distribution of this group was greatly influenced by climate changes during the late Neogene [30]. Therefore, we hypothesize that the current distribution and endangered status of the Chinese white pine species may have been greatly but differently influenced by the climate variation during the last glacial period due to the different ecological niches of the species. Under projected future global warming, the suitable habitat of each species is expected to change differently, which may raise serious conservation issues for some species, especially those with restricted distribution.
In this study, we used comprehensive occurrence data of the Chinese white pine species and a species distribution model to dynamically project and quantify their suitable habitats under paleoclimatic, current, and future climate scenarios. We also explored the centroid shift, range expansion/contraction, and suitability change of the current distribution area to assess the potential risk to each species in the future. Our aims are to find out how climate fluctuation since the last Due to its importance in ecosystems, P. subsect. Strobus has attracted much attention in recent years. Previous studies have studied white pines in terms of their phylogeny, ecology, population structures of endangered species, etc. [30,[35][36][37][38][39][40][41][42][43][44][45]. However, it is still unclear what are the key ecological factors determining the current distribution of P. subsect. Strobus, how the distribution pattern reacted to changes in paleoclimate, and how the distribution pattern will react under future climate changes. This knowledge will help us to identify response mechanisms, which will contribute to the development of effective conservation strategies in the future.
The maximum entropy method (MaxEnt) is a species distribution model (SDM) that is used to simulate the habitats that are suitable for species based on presence data and information of environmental variables [46]. Due to the model's effectiveness, it has been widely used in biodiversity (especially for endangered species) conservation [47][48][49][50][51], invasive species control [52,53], as well as nature reserve planning and management [54][55][56]. It has also been applied to reconstruct habitat shifts of individual species or vegetation during the last glacial period, and simulate how paleoclimatic fluctuations have influenced the present distribution [57][58][59], predict the potential distribution in the future, find conservation gaps, and guide protection management [60][61][62].
Niche conservatism is a key assumption underlying the application of SDMs [63]. According to a phylogenetic study of P. subsect. Strobus species, the current species diversity and distribution of this group was greatly influenced by climate changes during the late Neogene [30]. Therefore, we hypothesize that the current distribution and endangered status of the Chinese white pine species may have been greatly but differently influenced by the climate variation during the last glacial period due to the different ecological niches of the species. Under projected future global warming, the suitable habitat of each species is expected to change differently, which may raise serious conservation issues for some species, especially those with restricted distribution.
In this study, we used comprehensive occurrence data of the Chinese white pine species and a species distribution model to dynamically project and quantify their suitable habitats under paleoclimatic, current, and future climate scenarios. We also explored the centroid shift, range expansion/contraction, and suitability change of the current distribution area to assess the potential risk to each species in the future. Our aims are to find out how climate fluctuation since the last interglacial period have impacted the present distribution pattern and endangerment of Chinese white pines in order to determine the vulnerability of the studied species to future climate change, identify the habitats at risk of extinction, and provide guidance for species conservation.

Study Species Field Survey and Occurrence Data
According to Flora of China [32], ten species of P. subsect. Strobus are naturally distributed in China. The authors treated P. dabeshanensis as a variety of P. fenzeliana, and treated P. mastersiana as a variety of P. armandii (with the acceptance of the treatment by Hayata) [64]. However, the distribution areas of both P. dabeshanensis and P. mastersiana are far from the distribution areas of P. fenzeliana and P. armandii, and additionally the former two taxa are genetically differentiated from P. fenzeliana and P. armandii [32,36]. Hence, we kept the species status of P. dabeshanensis and P. mastersiana and recognized 12 natural species of white pines in China.
From 2016 to 2018, we conducted extensive field surveys throughout the major distribution areas of these 12 species in mainland China and Taiwan Province. The geographical coordinates of the surveyed locations were recorded with an accuracy of 10 m. Additionally, online databases were consulted to obtain comprehensive occurrence data for the white pines, including the Chinese Virtual Herbarium database (CVH, http://www.cvh.ac.cn) and the Global Biodiversity Information Facility (GBIF, http://www.gbif.org). For CVH data, the specimens were checked, and misidentifications were corrected based on the taxonomic treatments mentioned above. When no exact geographical coordinates were recorded for a specimen, the longitude and latitude coordinates of the specimens were determined using Google Earth 7.0 (https://www.google.com/earth) [65,66]. We used the R package "dismo" [67] to retrieve the occurrence data from GBIF. The specimens were checked according to previous taxonomic studies to avoid the risk of misidentification [32,68,69]. For those species that are not endemic to China, entire natural distribution records were used for niche modeling. In order to avoid sampling bias [70], in the case of multiple occurrence sites, only one was retained in each 10 × 10 km square using the "spatially rarify occurrence data" tool in the SDM toolbox for the ArcGis software (Esri, Redlands, CA, USA) (http://www.sdmtoolbox.org) [71].
Finally, a total of 1459 records of the 12 white pine species were kept for analysis. Among them, there were only four distribution records for P. wangii, only 10 distribution records for P. dabeshanensis, and only 17 distribution records for P. bhutanica. For all of the other species, a much larger number of records were available (Table S2).

Environmental Variables
The distribution of white pines is profoundly influenced by environmental factors since their growth and reproduction requires specific combinations of climatic and geographic conditions. For example, P. koraiensis is a photophilous species which needs a cold moist climate and deep fertile soil, whereas P. pumila is shade-tolerant and can grow in shallow soil and P. fenzeliana usually grows in warm, wet mountainous areas with weakly acidic soil [31]. Therefore, we utilized 39 environmental variables for niche modeling, including temperature, precipitation, soil, ultraviolet radiation, terrain, and land cover. Nineteen widely used bioclimatic variables were accessed from the WorldClim online database (http://worldclim.org/data/worldclim21.html) [72]. Four climatic variables-namely wet days, ground frost frequency, water vapor, and cloud cover-were obtained from the IPCC data distribution center (http://ipcc-data.org/observ/clim/get_30yr_means.html). We obtained another three variables (growing degree days, soil pH, and soil organic carbon) from the Center for Sustainability and the Global Environment, Nelson Institute for Environmental Studies, University of Wisconsin-Madison (UW) (http://nelson.wisc.edu/sage/data-and-models/atlas/maps.php). Additionally, we downloaded six UVB variables from glUV, an online global UVB radiation dataset (https://www.ufz.de/gluv/index. php?en=32435) [73]. The other five geographic variables were derived from a dataset provided by the Food and Agriculture Organization of the United Nations (FAO) (http://www.fao.org/soils-portal/ soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/zh/) [74]. All of the 39 environmental variables are listed in Table S3.
Collinearity may exist among the environmental variables, which may exaggerate the importance of highly correlated variables and hence may lead to over-fitting of the simulation [75,76]. Therefore, we determined the Pearson correlation coefficients among the 39 variables using the "banding collection statistics" tool in ArcGis 10.7 [77] and conducted the Jackknife test using Maxent to evaluate the contributions of each variable. All of the 39 variables scored 0 contribution values for all of the 12 species, and the variables with less ecological importance in highly correlated (r ≥ 0.9) variable pairs were excluded [77]. Finally, 18 variables were used in the subsequent analysis (Table 1).

Assessment of Important Environmental Factors
Three methods were used to estimate the importance of the environmental variables, namely percentage contribution, permutation importance, and Jackknife of regularized training gain. The percentage contribution value was an average value derived from each of the 10 training replicates and depended on the specific path that the model selects to find the optimal solution [78]. The permutation importance value, which depends more on the model's final performance instead of an individual path, was obtained by switching the variable values between occurrence and background spots and observing the change leading to the area under the curve (AUC) value [60]. The Jackknife test was conducted by creating a series of new models with exclusion or inclusion of a variable in turn to reveal the change in the information content of the model that is caused by each variable [78].

Paleoclimate and Future Climate Scenarios
Aside from the current scenario, three paleoclimatic scenarios were modeled: the last interglacial period (LIG, 120 to 140 ka), the LGM (22 ka), and the Middle Holocene (MH, 6 ka). These scenarios were modeled in the Community Climate System Model version 4 (CCSM4), which is an efficient climate model developed by the National Center for Atmosphere Research (Boulder, CO, USA) [79] and is widely used in studies on the response of species distribution to past and future climate changes [65,66,80].
We also used the CCSM4 model under two typical CO 2 representative concentration pathways (RCPs), namely RCP2.6 and RCP8.5, which represent two future climate scenarios depending on different levels of greenhouse gas (GHG) emissions, with a projected average temperature increase of 1.0 • C and 3.7 • C respectively, in 2081-2100 relative to the situation in 1986-2005 [8]. Then, the future climate datasets were applied to four combined scenarios: RCP2.6-2050 (average for the period 2041-2060), RCP2.6-2070 (average for the period 2061-2080), RCP8.5-2050, and RCP8.5-2070. All of the paleoclimatic and future data were downloaded from the WorldClim online database.

Model Simulation and Model Performance Evaluation
We used Maxent 3.4.0 (http://www.cs.princeton.edu/~schapire/maxent/) to construct the ecological niche and predict the potential distribution of each species under the current climate, three paleoclimatic climate scenarios, and four future climate scenarios [81]. All of the environmental variables were converted into ASCII files with the same resolution (2.5 × 2.5 , approximately 4.6 km × 4.6 km) in ArcGis 10.7 as required by the software. The model was parameterized with a maximum of 10,000 background points, a convergence threshold of 0.00001, and a maximum of 500 interactions. To optimize the model, we ran 10 bootstrap replicates, in which 75% of the presence data was randomly selected for model training and the other 25% was used for model testing. During the modeling, response curves were generated for each individual variable in order to demonstrate how the occurrence probability was influenced by their variations [48,65,66,70].
Additionally, two assessment indicators were adopted to measure the accuracy of the SDM, namely the AUC of the ROC [82] and the TSS [83]. ROC curves were generated by plotting sensitivity (true positive rate) against false positive rate. The AUC is one of the most widely used threshold-independent accuracy measures for SDM [84][85][86]. TSS is defined as the sum of sensitivity and specificity (true negative rate) minus 1, with a range of -1 to 1. As TSS is a threshold-dependent index used for binary prediction data, the continuous predictions from the Maxent model were transformed into binary ones by setting the threshold with the maxTSS method [87,88]. It is usually regarded that AUC values of 0.5-0.7 are an indicator of poor performance, values of 0.7-0.9 are an indicator of medium performance, and values of 0.9-1.0 are an indicator of good performance [89]; meanwhile, it is usually considered that TSS values from 0.7-0.85 indicate good performance and values of 0.85-1.0 indicate excellent performance [90].

Potential Distribution Areas Classifying, Suitable Habitat Changes in Different Scenarios
The output maps of Maxent were ASCII grid maps with the value of occurrence probability (ranging from 0 to 1) of the species for each grid cell (2.5 × 2.5 ). According to occurrence probabilities, the simulated potential distribution areas were categorized into four classes of suitability: not suitable (<0.2), low suitability (0.2-0.4), medium suitability (0.4-0.6), and high suitability (>0.6); areas ranked in the latter three classes were considered as suitable habitats [11,66,91]. The suitable areas for each species within China were calculated. A map of China (1:4,000,000) showing the administrative regions was downloaded from the Resource and Environment Data Cloud Platform (http://www.resdc. cn/Default.aspx) provided by the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, and spatial data for other countries were obtained from the GADM online database (https://www.gadm.org/).
Distribution changes for the past and the future climate scenarios were analyzed by calculating the interactions and complements between the distribution areas of each period and the current scenario. In order to examine how the suitable habitats change altitudinally and latitudinally in response to climate changes under the multiple scenarios in detail, we extracted the center point of each grid cell in the suitable areas (in Asia) of each species with ArcGis 10.7 and plotted the elevation and latitude of each point in a scatter diagram [47,92].
To understand how the suitability of the current actual distribution sites (presence spots) of each species in China change in the future scenarios, we extracted the probability values of the grid cells where the current presence points were located and calculated the proportions of points with the four different suitability ranks (not suitable, and low, medium, and high suitability) under the current climate and the four future scenarios. The proportion was calculated using the following equation: where Ni is the number of presence points within China with a certain suitability rank, and S denotes the total number of actual presence points within China of each species. We estimated the species' distribution centroid of each scenario using the "spatial statistics" toolbox in ArcGis 10.7 by calculating the occurrence probability-weighted Euclidean distances among potential occurrence grids [52,93].

Model Accuracy and Key Environmental Variables
The mean AUC and TSS values of the 12 white pine species for the current scenario were all above 0.9 and above 0.7, respectively ( Table 2), indicating that the models had a high level of accuracy and validity. The percentage contribution, permutation importance, and Jackknife of the regularized training gain results are shown in Table 3. According to the response curves of the environmental variables ( Figure S1), we obtained the range (occurrence probability > 0.2) of each variable for the 12 species (Table S4). The top four most contributing variables and permutation importance values of each species cumulatively account for more than 70% of the total contributions, except for P. wangii. Our results demonstrated that the precipitation in the warmest quarter (Bio18) is a dominant factor for all the tested species except P. wallichiana and P. sibirica. According to response curves of Bio18 ( Figure S1), the two Taiwanese species are the most moist adaptive (with suitable rainfall of over 800 mm in wettest season, Table S4); P. fenzeliana, P. kwangtungensis, P. armandii, P. koraiensis, and P. pumila can all grow under a wide range of rainfall levels; however, the first two species cannot endure rainfall of under 500 mm in the rainy season. The other three precipitation variables (Bio12, Bio13, and Bio16) were also found to be important for many species, such as P. bhutanica, P. dabeshanensis, P. fenzeliana, P. kwangtungensis, P. mastersiana, P. morrisonicola, and P. pumila (Table 3 ) are also important for many white pine species. Among them, temperature seasonality was found to be important for all species except P. bhutanica and P. dabeshanensis. Two UVB factors (UVB1 and UVB3) contributed a large amount for the two Himalayan species (P. wallichiana and P. bhutanica), two northern species (P. koraiensis and P. pumila), and the endemic P. armandii and P. dabeshanensis.

Suitable Habitats in Current Condition
The area of suitable habitats (high, medium, and low suitability) within China of the 12 species under the current climate, three paleoclimatic climate scenarios, and four future climate scenarios are shown in Table 4, Figure 2, Table S5 and Figure S2. For each species, the predicted suitable habitat was generally larger than the actual distribution range. Of the three northern white pines, P. sibirica is predicted to occupy the widest suitable habitat in the Eurasian continent, mainly in Russia (approximately 50 •~6 0 • N, 40 •~1 10 • E), with only small marginal areas extending into China (the Altai Mountains and Tianshan Mountains in Xinjiang, and the eastern margin of Northeast China). The suitable habitat of P. koraiensis spreads from the Changbai Mountains of Northeast China to the Korean Peninsula, far-eastern Russia, and Japan. The suitable habitat of P. pumila is similar to that of P. koraiensis, but extends further north to the Greater Khingan Mountains, the Lesser Khingan Mountains, and the Kamchatka Peninsula. The largest predicted suitable habitat (1,419,288 km 2 ) of the 12 species was obtained for the most widely distributed endemic species, P. armandii, which was predicted to cover a vast area of China, including the Qinling-Daba Mountains, the Yungui Plateau, and the Hengduan Mountains to the eastern Himalaya (Table 4). The suitable habitats of the two endangered species, P. dabeshanensis and P. wangii, are relatively limited and fragmented in China. Those of P. dabeshanensis are scattered either along the Yangtze River or in coastal areas of South and Southeast China. For P. wangii, the suitable areas occur in small patches in southwestern and southern mainland China and Southwestern Taiwan. For these two species, the predicted suitable areas are much larger than their actual distribution areas.

Suitable Habitats in Current Condition
The area of suitable habitats (high, medium, and low suitability) within China of the 12 species under the current climate, three paleoclimatic climate scenarios, and four future climate scenarios are shown in Tables 4 and S5 and Figures 2 and S2. For each species, the predicted suitable habitat was generally larger than the actual distribution range. Of the three northern white pines, P. sibirica is predicted to occupy the widest suitable habitat in the Eurasian continent, mainly in Russia (approximately 50°~60° N, 40°~110° E), with only small marginal areas extending into China (the Altai Mountains and Tianshan Mountains in Xinjiang, and the eastern margin of Northeast China). The suitable habitat of P. koraiensis spreads from the Changbai Mountains of Northeast China to the Korean Peninsula, far-eastern Russia, and Japan. The suitable habitat of P. pumila is similar to that of P. koraiensis, but extends further north to the Greater Khingan Mountains, the Lesser Khingan Mountains, and the Kamchatka Peninsula. The largest predicted suitable habitat (1,419,288 km 2 ) of the 12 species was obtained for the most widely distributed endemic species, P. armandii, which was predicted to cover a vast area of China, including the Qinling-Daba Mountains, the Yungui Plateau, and the Hengduan Mountains to the eastern Himalaya ( Table 4). The suitable habitats of the two endangered species, P. dabeshanensis and P. wangii, are relatively limited and fragmented in China. Those of P. dabeshanensis are scattered either along the Yangtze River or in coastal areas of South and Southeast China. For P. wangii, the suitable areas occur in small patches in southwestern and southern mainland China and Southwestern Taiwan. For these two species, the predicted suitable areas are much larger than their actual distribution areas.  Compared to the other 10 species, P. fenzeliana and P. kwangtungensis have wider and more continuous suitable habitats in vast areas south of the Qinling Mountains. The predicted suitable habitats of these two species are highly similar, except that there is less suitable habitat for P. fenzeliana in Hunan and Jiangxi Provinces. The two Himalayan species, P. wallichiana and P. bhutanica, also have similar suitable habitats in China; their suitable habitats mainly occur in the Hengduan Mountains and along the Himalayas. The suitable areas of the two Taiwanese species, P. mastersiana and P. morrisonicola, are also similar and mainly occur in the Central Mountain Range of Taiwan with a small area along the coast of South China and Northeast India. However, the suitable area of P. mastersiana is smaller than that of P. morrisonicola.

Distribution Shifts in Different Climate Scenarios
Predicted changes in suitable habitat for different scenarios are presented in Figure 3 and Table 4, Figures S3-S9 and Tables S5 and S6. The suitable area of each species underwent great change during the LIG and LGM periods. In the LIG period, compared to the current scenario, the suitable habitats in China of P. dabeshanensis, P. bhutanica, and P. koraiensis significantly expanded, whereas the two Taiwanese species and P. wangii lost much of their suitable habitat (Figure 3, Table 4). Additionally, the results showed that most of the white pine species experienced great changes in their suitable habitat area during the last glacial maximum. In this period, the suitable areas of seven species contracted notably (the two Taiwanese species as well as P. bhutanica, P. dabeshanensis, P. fenzeliana, P. koraiensis, and P. kwangtungensis). On the contrary, the suitable area of P. sibirica increased in China during the LGM (Figure 3). A large range expansion was predicted during the MH for the two Taiwanese species, P. dabeshanensis, P. pumila, P. fenzeliana, P. koraiensis, and P. kwangtungensis. For the future scenarios, the suitable areas of the two Taiwanese species and P. armandii were predicted to decrease. Meanwhile, the suitable areas of the nine other species were predicted to increase in most of the scenarios except P. sibirica and P. wangii for RCP2.6-2050, P. dabeshanensis for RCP2.6-2070, and P. sibirica for RCP8.5-2050. (Figure 3, Table 4).  The shifts in the distribution of the suitable habitats were further simulated by calculating the centroids of each suitable habitat for different climate scenarios (Table S7, Figure S10) as well as the altitude and latitude change of the suitable habitats in each grid cell under climate changes ( Figures  4 and 5, and Figures S11-S15). The suitable habitats of the three northern species displayed obvious changes in latitude during the last interglacial/glacial periods. In the LIG scenario, the suitable habitats of P. sibirica, P. koraiensis, and P. pumila expanded to higher latitudes compared to the current condition. Their habitats also tended to spread to lower latitudes at higher elevations ( Figure S11). In the LGM, at higher latitudes, the habitats of the three northern species contracted (Figure 4). Meanwhile, in the MH, at higher latitudes, the habitats of these three species expanded compared to the LGM scenario. (Figure S12). The shifts in the distribution of the suitable habitats were further simulated by calculating the centroids of each suitable habitat for different climate scenarios (Table S7, Figure S10) as well as the altitude and latitude change of the suitable habitats in each grid cell under climate changes (Figures 4 and 5, and Figures S11-S15). The suitable habitats of the three northern species displayed obvious changes in latitude during the last interglacial/glacial periods. In the LIG scenario, the suitable habitats of P. sibirica, P. koraiensis, and P. pumila expanded to higher latitudes compared to the current condition. Their habitats also tended to spread to lower latitudes at higher elevations ( Figure S11). In the LGM, at higher latitudes, the habitats of the three northern species contracted (Figure 4). Meanwhile, in the MH, at higher latitudes, the habitats of these three species expanded compared to the LGM scenario. (Figure S12).   During the last interglacial/glacial periods, two southern species, namely P. fenzeliana and P. kwangtungensis, showed similar habitat changes to the northern species in response to climate change. In the LIG, these two species were mainly distributed between 18-30° N ( Figure S11). Meanwhile, in the LGM, the suitable habitats of these two species were greatly reduced ( Figure 4). Subsequently, during the MH, their suitable habitats expanded greatly at latitudes of 18-34° N ( Figure S12). The suitable habitat of the Himalayan species P. bhutanica showed altitudinal shifts during the last glacial period. In the LIG, the suitable habitat of P. bhutanica extended to elevations of up to 5600 m, while in the LGM, the suitable habitat retreated to lower elevations. The elevation of the suitable habitat did not increase in the MH (Figures S11 and S12). The other Himalayan species, P. wallichiana, did not show clear altitudinal and latitudinal habitat shifts during the last glacial period. In the LGM, the During the last interglacial/glacial periods, two southern species, namely P. fenzeliana and P. kwangtungensis, showed similar habitat changes to the northern species in response to climate change. In the LIG, these two species were mainly distributed between 18-30 • N ( Figure S11). Meanwhile, in the LGM, the suitable habitats of these two species were greatly reduced ( Figure 4). Subsequently, during the MH, their suitable habitats expanded greatly at latitudes of 18-34 • N ( Figure S12). The suitable habitat of the Himalayan species P. bhutanica showed altitudinal shifts during the last glacial period. In the LIG, the suitable habitat of P. bhutanica extended to elevations of up to 5600 m, while in the LGM, the suitable habitat retreated to lower elevations. The elevation of the suitable habitat did not increase in the MH (Figures S11 and S12). The other Himalayan species, P. wallichiana, did not show clear altitudinal and latitudinal habitat shifts during the last glacial period. In the LGM, the suitable habitat of P. wallichiana covered a much narrower range of both altitude and latitude compared to the LIG and MH. Similar results were obtained for the two Taiwanese species and the endangered P. dabeshanensis. An interesting result was obtained for the endangered species, P. wangii, whose habitat was predicted to expand to higher latitudes and higher elevations during the LGM relative to the LIG (Figure 4). The predicted habitat change of the most widely distributed endemic species, P. armandii, was unique in that its suitable habitat was predicted to keep moving northward after the LIG, even in the LGM.
In the four future climate scenarios, the predicted habitat changes of the 12 species were largely similar. For three of the 12 Chinese white pine species, namely P. armandii, P. mastersiana, and P. morrisonicola, the suitable habitat was predicted to reduce in all the future scenarios except P. morrisonicola for RCP2.6-2050, whereas the suitable habitats of the nine other species were predicted to expand (Table 4).
In all the future scenarios, suitable habitats of P. armandii were predicted to spread to the north and to higher elevations ( Figure 5 and Figures S13-S15). This was also the case for the two Himalayan species P. bhutanica and P. wallichiana, as well as a northern species P. pumila and the endangered P. wangii. On the other hand, in all of the future scenarios, P. koraiensis, P. fenzeliana and P. kwangtungensis were predicted to spread to higher latitudes but not to higher elevations ( Figure 5 and Figures S13-S15).
In all of the future scenarios, the suitable habitat of P. mastersiana was predicted to reduce to a very narrow latitudinal range ( Figure 5 and Figures S13-S15). Meanwhile, for the other Taiwanese species, P. morrisonicola, the suitable habitat was predicted to move to the north in the RCP2.6-2050, RCP2.6-2070, and RCP8.5-2050 scenarios; however, in the RCP8.5-2070 scenario, the range was predicted to reduce to a very narrow latitude range, similar to P. mastersiana. Pinus sibirica was also predicted to have obvious range shifts in all of the future scenarios: at low latitudes, its suitable habitat was predicted to move to higher elevations, while in suitable habitats at low elevation, it was predicted to move to higher latitudes ( Figure 5). For the vulnerable species, P. dabeshanensis, no change in its suitable habitat (either altitudinally or latitudinally) was predicted in any of the future scenarios.
In this study, we also checked the suitability change of the present distribution sites in the future scenarios to predict the potential risk to each species. The proportions of presence points with the four different suitability ranks under the current and the four future scenarios are shown in Table 5. There are six species whose proportion of not suitable presence points was predicted to increase in the future, namely P. armandii, P. kwangtungensis, P. mastersiana, P. morrisonicola, P. sibirica, and P. wallichiana. This proportion was predicted to decrease in only one species, P. koraiensis, for which all of the actual sites were predicted to become suitable in all four of the future scenarios.

Discussion
This study predicted the habitat suitability of the 12 Chinese white pine species in the current climate and under paleoclimatic climate scenarios and future climate scenarios using Maxent. The AUC and TSS values produced by the models adequately described the distribution patterns of the 12 species. The results showed high precision, indicating that our modeling was reliable and accurate for the prediction of the suitable habitats of the white pine species. This study revealed dynamic range shifts of the Chinese white pine species under climate changes, which may be useful for conservation strategies.
The predicted current suitable habitat was basically consistent with, albeit wider than, the actual distribution of each species in China (Figure 2). This indicates that there may be more areas in potentially suitable habitats for the ex situ conservation of these species. However, small and discontinuous suitable habitats (such as for P. sibirica and P. dabeshanensis) may be one of the key factors for the endangered situation of these species [94,95]. For example, although the endangered species, P. dabeshanensis and P. wangii, have a suitable habitat outside of their actual distribution areas, these suitable habitats are still very narrow and dispersed. This is also the case for the two Taiwanese species, which have some small patches of suitable area outside of Taiwan Island. Although the three northern species have large suitable habitats, the majority of this habitat is in countries neighboring China. Of the 12 Chinese white pine species, only three (P. armandii, P. fenzeliana, and P. kwangtungensis) were predicted to have relatively large and continuous habitats in China.
The adaption of species to different environmental variables is a determinant factor for their distribution patterns [96]. The 12 species of Chinese white pine respond to environmental variables in different ways (Table S4), reflecting their differentiation in the ecological niche. According to the contribution indexes obtained from Maxent modeling and the Jackknife test, it was found that the precipitation of the warmest quarter (Bio18) and temperature seasonality are the factors that contribute the most to the species distribution. The importance of Bio18 indicates that most of the Chinese species are adapted to the East Asian monsoon climate with warm and moist summers [97,98]. Previous studies have suggested that Bio18 may be the most important factor for the growth or biomass accumulation of the Chinese white pines [99,100].
The four important temperature variables (Bio2, Bio4, Bio6, and Bio9) highlighted the contribution of minimum temperature and temperature variations on the distribution of the Chinese white pines, probably due to the fact that, in China, the summer is commonly warm and the minimum temperature and temperature variation in winter differ greatly with latitude [97]. Temperature seasonality generally increases with increasing latitude [101], and strong seasonal variation of temperature may inhibit vegetative growth [102]. As a result, temperature seasonality may determine the upper latitude limit for the distribution of each species. Bio6 may determine whether the temperature meets the need of vernalization for the northern species such as P. sibirica [103] and whether the pines from low-latitude areas can survive the harsh weather during the winter [104]. The response curves ( Figure S1) of these temperature variables show that, compared with the northern species, the two southern species (P. fenzeliana and P. kwangtungensis) and the two Taiwanese species are more adapted to moisture and high temperature and prefer stable thermal conditions with low diurnal and seasonal temperature fluctuations (Table S4, Figure S1).
Another dominant factor affecting the distribution of Chinese white pines is UVB radiation, which is a damaging wavelength of radiation that can cause problems to plants such as DNA damage and growth rate reduction [105]. Plants' responses to UVB radiation vary with species [106], and some alpine plants exhibit adaption to the long-term enhanced UVB radiation at high elevation [107]. In this study, the UVB radiation levels that the two Himalayan species (P. bhutanica and P. wallichiana) can tolerate were found to be remarkably higher than those of the other studied species ( Figure S1). This may be an adaptation to the strong UV radiation in the alpine habitat of the Himalayas and the Hengduan Mountains [108].
Climate changes exert a significant influence on the geographical distribution of organisms, and consequently, the distribution patterns of species can reflect climatic conditions [14,109,110]. Global climate has fluctuated greatly since the Late Pleistocene, especially the last glacial-interglacial cycle, which has deeply influenced the current distribution of most living organisms [13,14,19]. During the glacial periods, the global temperature was 5-12 • C lower than today, and the environment of East Asia became much drier due to the weakness of the Asian summer monsoon [17]. Under these conditions, forest biomes in the Northern Hemisphere migrated southward and the distribution of many temperate species contracted remarkably [14,17].
The results of this study suggest that the distribution areas of most of the Chinese white pine species underwent evident fluctuations between the LIG and the present. According to our results, the suitable habitat of the northern species (P. sibirica, P. pumila, and P. koraiensis) fluctuated in response to glacier movement since the LIG (Figures 4 and 5, Figures S11 and S12). During the LGM, the polar ice sheets spread considerably [14]. Probably forced by this spread, the distribution of the boreal species P. sibirica contracted to a large degree in Northern Eurasia and retreated to lower-latitude regions in Europe (not shown in this study), Northeast China, the Korean Peninsula, and Japan ( Figure S4). As the last glaciers receded, its habitat expanded to the north, which may have affected the current distribution pattern in China; namely, a few limited populations in the Altai Mountains and the northern Greater Khingan Mountains (Figure 1). A similar scenario occurred for the other two northern white pine species, indicating that the last glacial period greatly affected the current distribution of white pine species in Northeastern Asia [58].
The four species that are most adapted to warm and moist conditions-namely P. mastersiana, P. morrisonicola, P. fenzeliana, and P. kwangtungensis-lost a tremendous amount of habitat during the LGM, when the habitat area decreased to less than half of that during the LIG (Table 4). Although the direct impact of the last interglacial/glacial period was smaller in South China [111], the habitats of the above four species also greatly reduced in this region during this period, indicating that they were vulnerable to the cool and dry climate [50,63]. Pinus fenzeliana and P. kwangtungensis inhabit extremely narrow and fragmented habitats in the mountainous areas of South China, which are characterized by high biodiversity and abundant endemic plant species and are considered to be typical refugia during glacial-interglacial cycles [112,113]. Notably, in the present study, the coastal areas of Southeast China were estimated to be suitable habitats for the above two species. This is supported by evidence from pollen fossils, which indicates that many upland conifer species of South China, including Pinus species, were common in coastal areas in South China during the LGM [114]. These regions are still predicted to be suitable for the two southern white pine species for the current condition; however, no actual distribution is recorded there. This indicates that the distribution of these species in the coastal regions of South China may have been inhibited before the glacial period and contracted during the LGM but that the two species failed to recover in these regions after glacial retreat. Mountains in Guangxi Province and the Nanling Mountains may be refugia for P. kwangtungensis, and these regions have been found to be genetic diversity hotspots of this species via mitochondrial DNA analysis [38]. The LGM climate and the species' response to this climate, especially the resulting refugia pattern, may have been an important cause of the current highly narrow and disjunctive distribution of the species.
Except for the abovementioned species, the climate of the LGM influenced the Chinese white pines in different ways. The suitable habitat of one Himalayan species, P. bhutanica, showed altitudinal shifts during the climate fluctuation of the last glacial period, whereas that of the other Himalayan species, P. wallichiana, narrowed both in altitude and in latitude (Figure 4, Figures S11 and S12). In the LIG, the vulnerable P. dabeshanensis had a large suitable habitat in China; however, this habitat shrank strikingly to a very limited area in the LGM. On the other hand, in the LIG, the endangered P. wangii had a very small habitat and then expanded its range in the LGM, although this expanded area was still limited. This result implies that these two species had different endangerment mechanisms. For P. dabeshanensis, the last glacial maximum may be the main determining factor for its current endangerment. A recent population genetic study showed that populations of P. dabeshanensis underwent a reduction about 10,000 years ago [43]. Whereas, for P. wangii, its current endangerment may have been influenced by factors which occurred prior to the LIG. Pinus armandii showed a significantly different response to climate changes since the LIG compared to other species. Between the LIG and the present time, its suitable habitat was relatively stable but continuously reduced. On the whole, its habitat continued to move northward even during the LGM and reached the Qinling Mountains at this time. Its suitable habitat expanded to the Loess Plateau during and after the LGM. This result is also supported by pollen fossil records [115]. The different response pattern of P. armandii may be due to the fact that this species occupied different ecological niches and adapted in a different way compared with other southern species such as P. fenzeliana, P. kwangtungensis, and P. mastersiana, and was more cold-tolerant than these species (Table S4). This also indicates that the Qinling Mountains, which run from west to east in China, may have been a geographical barrier against the ice sheet [116][117][118].
Predicting distribution changes that occurred under paleoclimate conditions could provide an effective way to predict distribution areas under future climate change [119]. The Fifth Assessment Report of the IPCC estimated that average global temperatures will have increased by 0.3-4.8 • C by the end of the 21st century [8]. This may profoundly influence the species distribution in temperate forests [10]. Shifts in suitable habitats may provide chances for species to adapt to predicted future climate change. However, the results of this study suggest that the habitats of three endemic species, namely P. armandii, P. mastersiana, and P. morrisonicola, will reduce in the future (Table 4). Among them, the two Taiwanese species were predicted to be in a critically endangered situation. The suitable habitats of both of these species were predicted to keep contracting to higher latitudes in the Central Mountains in Taiwan, where few new suitable regions are predicted to arise (Table 4, Figure 5 and Figures S13-S15). Meanwhile, regions at the southern edge and/or at low altitudes of the suitable habitats of P. armandii may disappear, and new suitable habitats may appear at higher latitudes and altitudes ( Figure 5).
The results of this study suggest that, under most of the future climate scenarios, the suitable habitats of the eight species other than P. armandii, P. mastersiana, and P. morrisonicola will generally expand to higher latitudes and higher elevations. However, this does not mean that these species will actually expand their distribution ranges in response to climate change, since this may be limited by the dispersal ability of the species and interspecific competition [120,121]. Furthermore, in the future, the actual distribution area of a species may contract when its current distribution habitats become unsuitable, even if the total suitable area has expanded. In the present study, we analyzed the predicted future changes in the suitability of the presence sites for each species within China. The results showed that, for the three species (P. armandii, P. mastersiana, and P. morrisonicola) whose suitable habitats are predicted to contract in the future, the suitability of their actual distribution sites is also predicted to reduce (Table 5). Therefore, P. armandii and the two Taiwanese species are predicted to have a high risk of population shrinkage or even local extinction in the future. In view of their restricted current distributions, P. mastersiana and P. morrisonicola require urgent protection in their entire distribution range. Additionally, for P. armandii, attention should be paid to locations at low altitude (below 2500 m) in the southern margin of the Hengduan Mountains, the southeast margin of the Yungui Plateau, the edge of the Sichuan Basin, and the southern margin of the Dabieshan Mountains ( Figures S6-S9).
Furthermore, the results of this study suggest that P. sibirica is also at high risk under future climate change, despite the fact that its suitable habitat within China is predicted to expand in three of the four future climate scenarios (Table 4). This species is extremely cold-adapted [42] and its distribution sites in China are at the southern and southeastern edges of its distribution range. Thus, it is likely that a high proportion of its current presence sites are not currently suitable and will continue to be unsuitable in the future. Furthermore, the results of this study suggest that, in the future, the suitable habitats of P. sibirica in China will be discontinuous ( Figures S6-S9), which will hinder the natural dispersal of this species.
The results of this study indicate that the suitable area of P. kwangtungensis will increase (Table 4) under all the future climate scenarios. However, a large area of its habitat in South China is expected to disappear (Figures S6-S9). Our results supported the suggestion of ex situ conservation of P. kwangtungensis propose by a previous study [39]. A similar situation was predicted for P. fenzeliana. Thus, it is reasonable that the IUCN has listed these two species as near-threatened [33]. For two species in Northeast China, P. pumila and P. koraiensis, no contraction of their suitable habitat in China was predicted and the suitability of their current presence sites in China was not predicted to decrease. P. koraiensis is the most important conifer species with high ecological and economic value in Northeast China [40]. Although P. koraiensis is currently listed in the second conservation class in China [34], this study suggests that its future situation will be better than those of P. dabeshanensis, P. wangii, and P. kwangtungensis, which are also listed in the second conservation class. Populations genetic study showed that high genetic diversity is present within the natural population of P. koraiensis [41], indicating that an in situ conservation strategy for this species should be taken into consideration.
The two Himalayan species, P. wallichiana and P. bhutanica, are listed as Least Concern by the IUCN. However, some of the presence sites of P. wallichiana are predicted to be unsuitable in the future. Given that the two Himalayan species are restricted, special care should be paid to the conservation of their natural populations. Additionally, more attention should be paid to the little-known species P. bhutanica, which has been the subject of only a few taxonomic studies to date [29,68,122]. Moreover, two endemic species, P. wangii and P. dabeshanensis, have a restricted distribution in China [35,45,123,124]. Although the suitable habitats of P. dabeshanensis are predicted to increase in the future, they are predicted to be discontinuous. Furthermore, for P. wangii, it is predicted that some of its few presence sites will disappear in the future. Therefore, both in situ and ex situ conservation strategies should be considered for these endangered species.
Combining paleoclimatic reconstruction with future prediction is of great importance for species conservation strategy in changing climates [63]. However, the uncertainties of our modeling should be kept in mind. For instance, it is difficult to implement some ecological and evolutionary processes which greatly affect the distribution of species in the modeling process, such as species competition and human activities [125]. Despite this, SDMs can still provide a reference which can be used to inform the choice of conservation strategies.

Conclusions
In this study, a species distribution model (SDM) not only revealed the possible paleoclimatic roots of the distribution patterns of Chinese white pine species but also shed light on appropriate directions for their conservation and sustainable management. The results of this study suggest that most of the Chinese white pine species responded sensitively to the last interglacial-glacial cycles and the ongoing global warming. The results suggest that P. armandii and the two Taiwanese species may be at high risk in China in the future, which indicates that both in situ and near situ protection strategies should be considered. Pinus kwangtungensis, P. dabeshanensis, and P. wangii are considered as threatened or endangered species, and therefore special care should be taken to enlarge their natural populations, especially for the latter two species. Additionally, this study suggests that P. sibirica is also at a high risk of losing its current populations. Given that the suitable habitat of this species is predicted to expand to several unconnected mountains outside of its current range, ex situ conservation should be considered for this species. Although P. koraiensis is listed as a second-class species for conservation, this study suggests that this species is not at risk under future climate change scenarios. Similar findings were also obtained for P. pumila and P. fenzeliana, for which the suitable area in China is predicted to increase in response to future global warming. The two Himalayan species are also rare within China and have small distribution areas. The SDMs revealed that, under future climate change, the suitable habitats of these species will expand in the Hengduan Mountains. However, the high mountains may impede the dispersal of white pine species. Efforts should be made to protect the current distribution ranges of these rare and fragile species and transplant them to wider suitable habitats.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/9/996/s1, Figure S1: Response curves of the key environmental variables; Figure S2: Bar chart of suitable habitat areas (of high, medium and low suitability, within China) under the current, three past and four future climate scenarios; Figure S3: Suitable habitat changes under the last inter-glacial (LIG) scenario compared with the current condition; Figure S4: Suitable habitat changes under the last glacial maximum (LGM) scenario compared with the current condition; Figure S5: Suitable habitat changes under the middle Holocene (MH) scenario compared with the current condition; Figure S6: Suitable habitat changes under the RCP2.6-2050 scenario compared with the current condition; Figure S7: Suitable habitat changes under the RCP2.6-2070 scenario compared with the current condition; Figure S8: Suitable habitat changes under the RCP8.5-2050 scenario compared with the current condition; Figure S9: Suitable habitat changes under the RCP8.5-2050 scenario compared with the current condition; Figure S10: Shifts of centroids for each species under different scenarios; Figure S11: Altitude and latitude change of the suitable habitats under the LIG scenario compared with the current condition; Figure S12: Altitude and latitude change of the suitable habitats under the MH scenario compared with the current condition; Figure S13: Altitude and latitude change of the suitable habitats under the RCP2.6-2050 scenario compared with the current condition; Figure S14: Altitude and latitude change of the suitable habitats under the RCP2.6-2070 scenario compared with the current; Figure S15: Altitude and latitude change of the suitable habitats under the RCP8.5-2050 scenario compared with the current condition; Table S1: Species of Chinese white pines and their conservation status; Table S2: Presence data of the 12 Chines white pine species; Table S3: Information of the 39 environmental variables; Table S4: Suitable ranges of the 18 environmental variables for the 12 Chines white pine species; Table S5: Area of suitable habitats (of high, medium and low suitability) within China under the current condition, three past and four future climate scenarios; Table S6: Suitable habitat area changes (within China) under the three past and four future climate scenarios compared with the current condition; Table S7: Geographic coordinates of the distribution centroids under the eight climate scenarios.