Distribution Pattern of Endangered Plant Semiliquidambar cathayensis (Hamamelidaceae) in Response to Climate Change after the Last Interglacial Period

Semiliquidambar cathayensis is a special and endangered plant in China, used for traditional Chinese medicine and in landscape applications. Predicting the impact of climate change on the distribution of S. cathayensis is crucial for its protection and the sustainable use of resources. We used the maximum entropy (MaxEnt) model optimized by the ENMeval data packet to analyze the potential geographic distribution changes of S. cathayensis in 12 provinces of Southern China for the different periods since the last interglacial period (LIG, 120–140 ka). Considering the potential geographic distribution changes in the province, and based on the two climate scenarios of Representative Concentration Pathways (RCP) 2.6 and RCP 8.5, the distribution range of S. cathayensis was analyzed and we predicted the range for the 2050s (average for 2041–2060) and 2070s (average for 2061–2080). The area under AUC (Area under the receiver operating characteristic (ROC) curve) is 0.9388 under these parameters, which indicates that the model is very accurate. We speculate that the glacial period refugia were the Nanling and Wuyi Mountains for S. cathayensis, and central and Western Fujian and Taiwan are likely to be the future climate refugia. In the mid-Holocene (MH, 6 ka), the growth habitat was 32.41% larger than the modern habitat; in the 2050s and 2070s (except RCP2.6–2070s), the growth habitat will shrink to varying degrees, so efforts to support its in situ and ex situ conservation are urgently needed. The jackknife test showed that the main factors affecting the geographical distribution of S. cathayensis were annual precipitation, precipitation of the wettest month, and precipitation of the driest month. The annual precipitation may be the key factor restricting the northward distribution of S. cathayensis. In general, the centroid of the distribution of S. cathayensis will move northward. The centroid of the adaptive habitats will move northward with the highest degree of climate abnormality. We think that Hainan Island is the most likely origin of S. cathayensis. These findings provide a theoretical basis for the establishment of genetic resources protection measures, the construction of core germplasm resources, and the study of the formation and evolution of Hamamelidaceae.


Introduction
The evolutionary direction, ecological habits, and geographical distribution of species are affected, constrained, and driven by the changing climate and human activities, which in turn affects ecosystem structure and biodiversity, except for human influence [1][2][3]. Periodic reciprocating oscillations, such as the glacial and interglacial periods and other sudden events in geological history, have led to regional differences in biodiversity; the last interglacial (LIG, 120-140 ka) is most closely related to modern times [4]. In the last glacial maximum (LGM, 22 ka), due to the rapid deterioration of the climate, several species have become extinct, and the distribution of most surviving species has shrunk sharply. As such, many biological refugia have appeared in many regions around the world [5][6][7]. The surviving terrestrial plants spread from their refuge under the conditions of post-glacial climate warming, and then redistributed. Climate change directly affects the distribution pattern of global plants, especially that of endangered plants [8][9][10]. Climate warming will gradually spread the centroid of the species distribution toward the north or the south as the suitable habitats migrate to higher latitudes and higher altitudes [11]. Li et al. used a combination of molecular markers and maximum entropy (MaxEnt) to determine how severe climate change may affect the spatial pattern of species distribution [12]. Juniperus drupacea, which is distributed in the Middle East and Greek Mountains, still had a wide habitats species in the LGM. With increasing climate change and human activities, habitat fragmentation is serious, and the species is now on the verge of extinction [13]. The extinction of species on earth is continually increasing, and 112,432 species have been listed on the International Union for Conservation of Nature (IUCN) Red List of threatened species [14]. The latest research reported that with the continuous development of industrialization, climate warming is becoming increasingly intense, and the frequency of El Niño and other extreme weather is also increasing [15]. China will have more extreme weather, and as endangered plants have more stringent requirements on habitat, studying the impact of climate change on the distribution of endangered plants is essential [16,17].
Species distribution models (SDMs), based on niche theory, are the main method used to evaluate the impact of climate change on species distribution using the known population distribution points of species and their associated climate factor data [18]. One of the various SDMs is maximum entropy (MaxEnt) [19][20][21]. However, each model has its own advantages and disadvantages. Through the research and comparison of various models, MaxEnt model was found to have better prediction ability, so is also one of the most widely used models [22,23]. At present, MaxEnt has been applied to a variety of aquatic and terrestrial animals to predict the change trend of potential distribution habitats [24], especially for invasive organisms, infectious diseases, and endangered species, such as Spartina alterniflora [25], Hex khasiana [26], and Thuja sutchuenensis [27]. MaxEnt can also use climate data to project to the past (such as the LIG and LGM) and the future to help determine the location of ice age shelters and the path of species migration and diffusion. MaxEnt can be used to provide auxiliary information for understanding species formation and differentiation, palynology, genetics, and species protection strategy formulation [28], which is important theoretical research and production practice.
Semiliquidambar cathayensis (Figure 1a-d) is a grade-2 protected plant in China. Study of its distribution pattern is essential for the phylogeny of Hamamelidaceae [29]. The whole plant is treasured in China. The roots, flowers, and leaves can be used in traditional medicine, with 85 kinds of active chemicals in the roots [30]. The extract and active components of the roots has inhibitory effects on cholinesterase and viral hepatitis antigen, respectively [31]. The plant is also one of the few colorful-leaved trees in South China and is used in landscape application. However, S. cathayensisis is on the verge of extinction, and most of whatever remains can be found in the south mountain area of China, with a fragmented distribution, so it needs to be protected. We found that the requirements for water, light, and heat of S. cathayensisis significantly differed between growth stages during the observation period in Tiantai mountain, Shaxian County (26 • 39'38" N, 117 • 49'32"E). The species likes shade and is sun-intolerant when young; it likes light and heat, and humid, fertile, and loose soil environment when adult. However, a wet winter and spring causes the seeds of S. cathayensisis to rot and deteriorate, resulting in few surrounding seedlings, indicating that climate change has considerably impacted the growth of S. cathayensisis. At present, reports are lacking on the response of Hamamelidaceae plants to climate and environment. The findings reported here can provide a reference for the protection and use of phylogeny of S. cathayensisis and Hamamelidaceae. In this study, based on the EMNeval-data-package-assisted selection of MaxEnt optimal parameter settings, the potential distribution habitats of S. cathayensisis in the LIG, LGM, middle Holocene (6 ka), modern, Representative Concentration Pathways (RCP)2.6-2050s, RCP2.6-2070, RCP8.5-2050s, and RCP8.5-2080s were simulated and inferred [32]. Applying percent contribution and permutation importance comparison, we synthetically used the jackknife test to evaluate the main factors restricting modern geographical distribution and used response curves to determine the appropriate interval of environmental variables. The objectives of this study were to: (1) reconstruct the history of the change in the geographical distribution pattern of S. cathayensisis since the LIG, (2) understand the restrictive mechanism of environmental factors on the potential geographic distribution of S. cathayensisis, (3) provide a scientific basis for the germplasm resources protection and management of S. cathayensisis, and (4) lay a theoretical foundation for the study on the formation and evolution of South China flora and Hamamelidaceae.

Study Area
S. cathayensisis is currently on the verge of extinction and remains only in the mountainous areas of Southern China [33]. The north boundary of its distribution is the Yangtze River, and the south boundary is Jianfengling National Nature Reserve of Hainan Island. The east boundary ranges from Wuyanling National Nature Reserve in Taishun County, Zhejiang, and Zhouning counties in Fujian and Daiyun Mountains Nature Reserve; the west boundary range from Cangyuan County in Yunnan-Guizhou Plateau to Gulin County in Sichuan [34,35]. The geographical coordinates of S. cathayensisis natural distribution area range from 18 • 41'38" to 29 • 7'44" N and 99 • 11'39" to 119 • 40'12" E.

Collection and Screening of Sample Data
In 2017-2019, our research team conducted field surveys of S. cathayensisis natural populations in 7 provinces of Fujian, Zhejiang, Jiangxi, Hunan, Guangxi, Guangdong, and Hainan, and collected 68 distribution records. We collected 17 distribution records of S. cathayensisis by consulting the published literature; we searched global species diversity information base [36], Specimen Resources Sharing Platform for Education [37], Chinese herbarium of nature [38], and Chinese Virtual Herbarium [39], collecting 13, Chinese herbarium of nature, 816, and 184 S. cathayensis specimen records, respectively, for a total of 1136. According to the method reported by Zhang et al. [40], all the above records were filtered, non-natural records were removed, and redundant data were deleted. To reduce the error caused by the cluster effect, each grid (10 × 10 km) retained only one distribution point, and finally 120 valid samples were obtained (Table S1; Figure 1e).

Environmental Variable Screening and Data Processing
The main determinant of species niche is climate variables, which are often used for plant niche modeling [41,42]. The LIG, LGM, MH, current and future climate data were downloaded from the WorldClim database [43]. Modern climate data are based on monthly meteorological data from different weather stations around the world from 1950 to 2000 and are generated by interpolation within a grid with a spatial resolution of 30 arc-second (about 1 km 2 on the ground).The climate data of the LGM, MH, and the year 2070 were generated by the Community Climate System Model version 4 (CCSM4), in which the future greenhouse gas emission scenario is a typical concentration target (RCP2.6 and RCP8.5) [44,45]. The selection of environmental factors is based on the essential variables of growth and development of S. cathayensis. First, bio5, bio6, and bio12 are selected. Then, the correlation analysis is carried out for all factors, and the most important environmental factor with a correlation greater than 0.8 is selected [46]. Then, the variance iflation factor (VIF) analysis is carried out for all factors, and the factor with a VIF value less than 10 is selected, and seven environmental variables were finally selected (Table 1) [47].

Model Establishment, Optimization, and Evaluation
In this study, MaxEnt3.4.1 [48] software was used to simulate and predict the potential geographical distribution probability of S. cathayensis under the LIG, the LGM, MH, current, and the 4 future scenarios (RCP2.6-2050s, RCP2.6-2070s, RCP8.5-2050s, and RCP8.5-2070s). To ensure the probability of S. cathayensis distribution appear close to normal, we selected 75% of the data for model training and the remaining data for model testing. Other values were set to default [49]. We used the package ENMeval in R v. 3.6.1 [50] to optimize the MaxEnt model, set the regulatory multiplier (RM) to 0.5-8, and each interval was 0.5, for a total of 16 regulatory multipliers [48,51]. We used 9 feature combinations (FCs): L, LQ, H, LQH, LQHP, LQHPT, QHP, QHPT, and HPT(The MaxEnt model provides 5 features, which are linear features (L), quadratic features (Q), segmented features (H), product features (P), and Threshold features (T)) [52,53]. The ENMeval data package was used to test the above 144 parameter combinations, and we finally used the Akaike Information Criterion (AIC) model of the Akaike information criterion, and used 10% training omission rate (OR 10 ) and the difference between the AUC values (AUC DIFF ) to check the fit and complexity of the model (Table S2) [50].
After the optimization was completed, the optimized parameters were used to simulate and predict the suitable growth habitat of S. cathayensis in different periods. The area under AUC (Area under the receiver operating characteristic (ROC) curve) was used to assess the accuracy of MaxEnt predictions [54]. The range of AUC values is (0, 1). The higher the AUC value, the higher the credibility of distinguishing between suitable and unsuitable habitats. An AUC less than 0.8 indicates that the reliability is low; 0.9-1.0 is extremely accurate [54]. We divided the suitability into four grades (I, 0%-20%; II, 20%-50%; III, 50%-70%; IV, 70%-100%), which respectively indicate unsuitable, less suitable, middle suitable, and highly suitable habitats, and the number of grids at each level was calculated to determine the percentage of area change.

Multi-Environment Similarity Surface (MESS) and the Most Dissimilar (MoD) Variable Analysis
Based on China's environmental variables as reference layers, the MESS and MoD were used to analyze the abnormal climate areas in the past and future situations and the key factors that caused the change of potential geographical distribution [32]. The multivariate similarity of each distribution point is the minimum of the similarity of each variable, and this variable is the MoD [55]. A negative multiple similarity value means that at least one variable value of the point exceeds the range of the reference layer, which is called a climate abnormal point and is expressed in blue. If it is positive, as the value becomes larger, the probability of a climate anomaly is lower, it is expressed in light blue, turquoise, and orange. The higher the score, the darker the color, the more normal the climate of the point. When the score is 100, the climate of the point is completely normal [56].

Analysis of Spatial Pattern Changes
We redefined spatial units with a species existence probability value of ≥0.5 as a suitable S. cathayensis habitat, and the spatial units with a species distribution probability value <0.5 as unsuitable habitats [32,57] so as to establish the existence / nonexistence (0,1) matrix of S. cathayensis under the past, current and future climate change scenarios. Based on the (0,1) matrix table, the pattern changes of suitable distribution areas of S. cathayensis under past and future climate change scenarios were further analyzed. The past and future area changes were calculated based on the current area of S. cathayensis suitable habitats. A matrix value 0 → 1 indicates the newly added suitable habitats, 1 → 0 indicates the lost suitable habitats, 1 → 1 is the reserved suitable habitats, and 0 → 0 is the unsuitable habitats. Finally, the matrix change values were loaded into Arc GIS 10.4.1 to create visualizations of the spatial pattern changes in suitable S. cathayensis habitats.

The core Distributional Shifts
The trends in the changes in suitable habitats were also calculated and the centroids were compared for the LGM, MH, and current and future suitable habitats using ZonalGeometry tool in Arc GIS 10.4.1, and the center points of the suitable regions in regions were compared [57]. We considered the suitable S. cathayensis habitat as a whole and reduced it to a vector particle and used the change of the position of the centroid to reflect the size and direction of the species' suitable habitat. Finally, we tracked the centroid with different SDMs to examine the centroid of S. cathayensis in different periods and under different climatic conditions to evaluate the migration distance of the suitable S. cathayensis zone in latitude and longitude coordinates [40].

Analysis of the Accuracy of the Model
Based on 120 distribution points and seven climate variables of S. cathayensis, we used the MaxEnt niche model to predict the potential geographical distribution of S. cathayensis. The default settings were RM = 1, FC = LQPHT, for which the minimum information criterion AICc value was 53.545. The difference between the training set and the test set AUC was 0.020149. When Feature Combination (FC) = LQ, Regularization Multiplier (RM) = 1.5, and ∆AICc = 0, the model was the optimal, and the difference between the model training set AUC and the test set AUC was 0.0198 (Table 2), which is 17.04% lower than with the default settings. The optimized ∆AICc and avg. diff. AUC (The average difference of AUC) values are lower than the default settings. The results showed that the optimized parameters reduce the fitting degree and complexity of the model. Therefore, FC = LQ and RM = 1.5 were selected as the parameter settings of the model. Under these parameter settings, we used MaxEnt to predict the modern potential distribution area. The average AUC value of 10 times of repeated training was 0.9388 and the average AUC value of the test data was 0.931, which shows that this model is accurate for S. cathayensis prediction.

Contribution Analysis of Environmental Variables
Using the response curve ( Figure 2), we obtained the thresholds of the seven main bioclimatic parameters (existence probability n > 0.2). These seven response curves were all single peak and nearly normal in distribution, which indicated that the prediction of S. cathayensis was successful. Annual mean temperature (bio1) ranged from 129 to 297 (12.9-29.7 • C), mean diurnal range (bio2) ranged from 39 to 98 (3.9-9.8 • C), mean temperature of wettest quarter (bio8) ranged from 149 to 346 (14.9-34.6 • C), annual precipitation (bio12) ranged from 1239 to 2940 mm, precipitation of the wettest month (bio13) ranged from 204 to 677 mm, precipitation of driest month (bio14) ranged from 20 to 61, and precipitation of warmest quarter (bio14) ranged from 431 to 1669 mm.

Environmental Variables Analysis
Among the seven environmental variables screened, the contributions of annual mean temperature (bio1, 37.3% of variation), mean diurnal range (bio2, 23% of variation), and annual precipitation (bio12, 17.8% of variation) were much larger than those of other environmental factors, and their cumulative contribution rate was as high as 78.1% (Table 1). The cumulative contribution rate of the temperature-dependent annual mean temperature and mean diurnal range was 60.3%. The top three important variables were annual mean temperature (bio1), annual precipitation (bio12), and mean temperature of wettest quarter (bio8), with a cumulative value of 64.4%. Figure 3 shows that when only individual variables were used, the regularization training gain, test gain, and the three variables with the highest AUC values were ranked the same. Annual mean temperature (bio12) and precipitation of warmest quarter (bio13) were the two variables that declined the most in the regularized training gain, test gain, and AUC values. Based on the above analysis, the main factors affecting the modern geographical distribution of S. cathayensis are precipitation factors (annual precipitation, precipitation of wettest month and precipitation of driest month) and temperature factors (annual mean temperature, mean diurnal range, and mean temperature of wettest quarter).

Potential Distribution Area in the Past
Our models predicted, during the LIG, the area of the middle-highly suitable habitats increased slightly compared with modern times, which was 369,965 km 2 , which is a growth of 7.38%. The original highly suitable habitats area of Hainan, Northeast Guangdong, and midwest Fujian became moderately suitable habitats or even unsuitable; the number of highly suitable areas in Hunan, Taiwan, and Jiangxi increased significantly ( Figure 4A). During the LGM, the middle to highly suitable habitats of S. cathayensis decreased significantly, leaving only 118,66 km 2 , which is only 3.44% of the modern medium-high suitable potential distribution area; the suitability of the suitable habitats at that time was significantly lower, and the original highly suitable habitats became moderately suitable or even unsuitable. The highly suitable habitats only included central and Southern Taiwan ( Figure 4B). In the MH, the area of suitable S. cathayensis habitat increased significantly compared with modern times, which was 456,208 km 2 , with is a growth of 32.41%. The highly suitable habitats in Hainan almost all became moderately suitable habitats, whereas the highly suitable habitats in Guangdong increased significantly. Highly suitable habitats appeared at the junction of Anhui, Zhejiang, and Jiangxi, and the area of low-grade suitable habitats in Yunnan increased slightly ( Figure 4C).
By comparing the four maps in Figure 5A-D, the areas that four historical periods were ecological stable the eastern part of Hainan Island and central Taiwan Island. Since the LIG, central and Western Fujian, Eastern Jiangxi, and Northern Guangdong have been stable ecological regions of S. cathayensis, and these regions may be the climate refugia S. cathayensis. In Northern Guangxi, Western Hunan, and Southern Hunan, the LGM and MH were middle-highly suitable S. cathayensis habitats, whereas in modern times, these habitats are low suitable. Combining Figures 2 and 4, we find that these regions, especially Guangxi and Guizhou, will gradually become unsuitable for S. cathayensis.

Future Changes in Suitable Habitat Area
Our model predicted the suitable habitat of S. cathayensis in four different climate scenarios in the future, and then conducts comparative analysis (Table 3, Figures 4 and 5). The results are as follows: Under the RCP2.6-2070s emission scenario, the habitats increase the most, by 31.64%, and the new habitats are mainly concentrated in the northwest of central and Western Hunan; the size of the middle-highly suitable habitats is about 419,908 km 2 , which is an increase of 21.88%, which indicates that many low-growth habitats under this scenario might be replaced by the middle-highly suitable habitats. This is the climate scenario where the future loss of habitat is the least, only 31,110.1 km 2 , and the loss rate is about 9.04%.
Under the RCP2.6-2050s emission scenario, the least amount of new habitat is added, with an increase rate of about 5.48%. The newly added habitats are mainly concentrated in Hainan Island, and the suitable habitats in other areas decrease, with an area of about 220,149 km 2 of middle-highly suitable habitat, about 63.90% of the area in the modern era. The area of suitable habitat loss is about 143,719 km 2 , the loss rate is about 41.76%, and the main loss habitats are Southern Hunan, Southern Guizhou, and Western, central, and Eastern Guangxi.
The distribution of newly added habitats is also larger under the RCP8.5-2050s emission scenario, but slightly smaller than under the RCP2.6-2070s emission scenario, which is an increase of 24.27%. The newly increased habitats are mainly concentrated in central Hunan and the border area between Hunan and Jiangxi; the middle to highly suitable habitats occupy about 321,938 km 2 , which is about 93.44% of the modern extent, and the suitable habitats shrink slightly compared with modern times. The habitat loss is large, concentrated in the border area of three provinces, Zhejiang, Anhui, and Jiangxi, and the border area of the two provinces of Guangdong and Jiangxi, about 107,810 km 2 , which is a loss rate of approximately 31.32%.
Under the RCP8.5-2070s scenario, the loss of habitat is the most severe, about 183,657 km 2 , which is a loss of 53.36%. The new suitable habitats under this scenario are relatively small, mainly concentrated in Southwestern Hunan, the border region of Hubei and Chongqing, and the border region of Hunan and Jiangxi. The increase is about 9.04%. The suitable habitats in Guangdong and Hainan are almost exhausted, as have the low-suitable habitats in Guangxi and Guizhou. This indicates that S. cathayensis suitable habitats will shrink extremely, and suitable habitat will only remain in most of Fujian and Northeastern Jiangxi under this scenario.

MESS and MoD Variable Analysis
In LIG, LGM, MH, RCP8.5-2050s, and RCP8.5-2070s, the average multivariate similarities of the 120 modern distribution points of S. cathayensis were 7. 63, 10.22, 5.19, 4.21, and 3.45, respectively. Only in the RCP8.5-2070s scenario did a distribution point with a negative multivariate similarity value indicate that the degree of climate anomaly is not optimistic. The most dissimilar variables were precipitation of wettest month, annual precipitation, and annual mean temperature ( Figure 6).

Shift in the Distribution Center of the Suitable Habitat
The current distribution center of S. cathayensis is located at 26 • 3'32.63" N and 114 • 50'19.35" E ( Figure 7B). The distribution centers of S. cathayensis during the LIG, LGM, and MH were located in Wuping County (24 • 58'20.

Bioclimatic Predictors and Model Performance
In this study, the results showed that the AUC given FC = LQ, RM = 1.5 for the parameter ROC curve was 0.9388, which is higher than that of Atrium brevifrons (0.901) [58], Tags lucida (0.92) [54], Homoia riparia (0.899) [59], and Phenacoccus madeirensis (0.9177) [60], which indicates that the simulation effect of the MaxEnt model of the potential geographical distribution area of S. cathayensis is accurate and reliable. In order to reduce errors, 144 combinations of features in ENMeval data packets were called in R software. In previous similar studies, the ENMeval data package was rarely used for optimization, and even the relevant studies that used ENMeval packets generally considered less than 96 parameter combinations. Therefore, after comparing 144 optimized combinations, the conclusions are relatively accurate [61].

Changes in Potential Distribution Areas
MaxEnt predicted that the highly suitable S. cathayensis habitat in the LIG was located in the central part of Taiwan Island and the eastern part of Hainan Island, but the suitable habitat was much smaller than that in contemporary times. This may have been caused by weather instability events, such as intense cold events, during the LIG. The average temperature of the extreme coldest month that can be tolerated by S. cathayensis is from 3 to 5 • C. Extreme cold events led to the rapid reduction of the suitable habitats of S. cathayensis. The LIG (120-140 ka) experienced cold events, and the average precipitation of the LIG was 557 mm, which was far lower than the requirements of S. cathayensis, which directly led to the prediction of the suitable habitats of S. cathayensis in the LIG being significantly smaller than currently [32]. Existing plants of the Hamamelis family originated in the eastern subtropical region, with an origin time in about the middle Cretaceous, and the northward diffusion occurred during the Upper Cretaceous, Paleozoic Eocene, and Oligocene. Molecular markers, related leaf fossils, and inflorescences have proven that S. cathayensis originated from the cross between Liquidambar and Altingia. Hainan Island is also an important distribution area for Liquidambar and Altingia [34]. Therefore, we speculate that Hainan Island, and even Southeast Asia, are likely to be important origin centers of S. cathayensis [62,63]. In the last glacial maximum, the Wuyi, Nanling, Hengshan, and Central Taiwan Mountain ranges were the most suitable habitats. It is generally thought that the large-scale climate shock in the last glacial maximum had a huge impact on the distribution of species, rapidly shrinking the distribution of species [64,65]. However, we found that LGM climate had little impact on the distribution of S. cathayensis, which is similar to Larrix laricina [66,67]. The results showed that the simulation error of LGM precipitation is larger than that of temperature, and the distribution of evergreen broad-leaved trees may be more sensitive to low precipitation than to low temperature due to the limitation of precipitation [68,69]. In the south of China, the LGM climate was much drier than the present, and the temperature, solar radiation, and other conditions were more unsuitable for the growth and development of S. cathayensis [70]. In the middle Holocene (MH, 6 ka), the suitable habitat of S. cathayensis was larger than during the other predicted periods, which may be related to the period (8.6-4.4 ka) when the East Asia hydrology created the most suitable habitat during the Holocene (13.1-2.5 ka), whereas the other periods of the Holocene were relatively dry [25]. The average annual temperature in the middle Holocene was slightly higher than that in the modern age, which led to the best survival period of S. cathayensis in the MH. The most suitable habitat for the growth of S. cathayensis obviously expanded to low altitude areas, which may be due to the increase in monsoon areas and precipitation by 10.7% and 18.7% respectively, during the MH compared with that in modern China [71].
The current middle-highly suitable habitat of S. cathayensis is 344,529.9 km 2 , and the highly suitable habitat is mainly concentrated in the middle and east of Hainan, mountainous areas in the middle and west of Fujian, and the Nanling Mountains in Guangdong. However, many large populations of S. cathayensis can be found in the Shiwandashan Nature Reserve of Guangxi, which is not in line with our expectations [33].
The findings showed that if the climate continues to warm, the potential distribution areas of S. cathayensis will shrink to varying degrees. We simulated the effects of RCP2.6 and RCP8.5 on S. cathayensis for two typical periods: 2050s and 2070s. Among them, under RCP2.6, S. cathayensis will first decrease and then increase; under the RCP8.5, the suitable habitat of S. cathayensis will continue to decrease and will be smaller than current. This result is consistent with the conclusion that, for the endangered plants of the same family Disanthus, the suitable habitats will shrink to different degrees [72]. The most suitable habitat, Nanling Mountain, will gradually shrink, especially in RCP8.5-2070s, with areas in Guangdong and Hainan, and the suitable habitat of S. cathayensis obviously moves northward. S. cathayensis is an excellent landscape tree species. As long as S. cathayensis is located in suitable habitats at present and in the future, many cities in Southern China can introduce S. cathayensis as urban street and ornamental trees.

Refugia from Climate Change
We classified the climate refugia of S. cathayensis as areas likely to have been or will be suitable from the LIG through to current and until at least the 2050s and 2070s. The Quaternary climate change had a profound impact on the historical space and population dynamics of the ancient S. cathayensis [73]. Although some errors may exist in the distribution during three historical periods since the LIG, we learned that the eastern part of Hainan Island, the central mountains of Taiwan, and Wuyi Mountains provide important biological shelters [74]. By comparing the changes under the five scenarios, the current, RCP2.6-2050s, RCP2.6-2070s, RCP8.5-2050s, and RCP8.5-2070s, we found that the central and western regions of Fujian (especially the Wuyi Mountains), the northeastern part of Jiangxi, and the middle part of Taiwan are middle to highly suitable habitats of S. cathayensis. We think that the above areas will be the future climate refugia for S. cathayensis.
However, we did not identify any distribution of S. cathayensis on Taiwan Island. The flora in Taiwan and mainland China are separated by the Taiwan Strait (200-410 km wide). With no artificial introduction, it is difficult to establish future climate refugia for S. cathayensis in Taiwan. We only predict possible climate refugia. Some of the above areas may not have S. cathayensis because other ecological factors will limit the distribution in these areas [73,75].

Restriction of Climate Variables on Geographical Distribution
Both the contribution rate and the permutation importance value indicate that the annual mean temperature and the annual precision are the most important variables affecting the distribution of S. cathayensis. The regularized training gain, test gain, and AUC value all indicates that annual precipitation (bio12), precision of wettest month (bio13), and precision of driest month (bio14) are the most important factors restricting the distribution of S. cathayensis. This is consistent with the results of the field surveys, where S. cathayensis populations were found to be mostly distributed on sunny slopes or semi-sunny slopes; with good water and fertilizer conditions; warm, humid, and loose soil; and in woodland environments, indicating that ambient temperature and humidity have a limiting effect on the distribution of S. cathayensis. Therefore, we think that the geographical distribution pattern of S. cathayensis is mainly restricted by two environmental factors: annual mean temperature (bio1, 12.9-29.7 • C) and annual precipitation (bio12, 1239-2940 mm). The results of MaxEnt analysis showed that the temperature and precipitation factors jointly limit the potential distribution pattern of S. cathayensis. The precipitation of driest month (bio14, 20-61 mm) is also one of the important factors affecting the distribution of S. cathayensis, both in terms of actual surveys and in the results of this simulation, as its seeds may decay in over-humid conditions per information from actual surveys or the simulation results.
However, climate factors do not affect the suitable S. cathayensis habitat independently. According to the MEES and MoD analyses, the climate change in RCP8.5-2070s is the most severe, and the suitability S. cathayensis habitat obviously reduces. The most dissimilar variable is annual mean temperature. The climate anomalies of RCP8.5-2050s, MH, LIG, and LGM decreased successively, with the most dissimilar variables being precipitation of wet month and annual precipitation. However, the change in the distribution of S. cathayensis is the result of multiple factors (especially during the LIG), and temperature and precipitation factors restrict the change in geographical distribution.

Protection of Genetic Resources
For the newly suitable habitats, as they are mostly distributed in the open forests on sunny slopes or semi-sunny slopes, a reasonable and sustainable land use plan should be formulated accordingly, and sufficient space should be reserved for its relocation. The newly suitable habitats are scattered in the marginal areas of the suitable habitats, which should reduce the human interference activities and increase the possibility of the immigration of S. cathayensis. The fruit of S. cathayensis is a capsule and its seeds are flat, light, and less viability, so artificial migration is required with its diffusion and colonization [76].
For the loss of suitable habitats, ex situ conservation measures should be taken, botanical gardens and core germplasm banks should be established, and S. cathayensis should be transplanted into an artificial environment for cultivation, conservation, and preservation. In addition, in situ protection is also important for the populations in the loss of suitable habitats, especially the small populations such as Leigong Mountain, Taiping Mountain, Yueliang Mountain, Le'an county, and Fanjing Mountain [77].
For the reserved suitable habitats can be protected to provide a safe place and refuge for S. cathayensis to cope with future climate change. Therefore, more attention should be paid to the protection and management of reserved areas. Establishing a nature reserve is the most effective method to protect rare and endangered wildlife resources in situ. At present, S. cathayensis protection forests have been established in Wuyi Mountain Nature Reserve, Jianfengling Nature Reserve, Dutianhe Nature Reserve, and Taiping Mountain Nature Reserve [33,34]. However, no protection measures have been taken in Liancheng, Shiwandashan, or Shaxian counties, where S. cathayensis is concentrated. The above areas can be prioritized as protected habitats for S. cathayensis.
In our team's field survey, S. cathayensis was rarely found to be concentrated; most of them were scattered or in isolated islands scattered in each mountain area, which complicates finding the distribution point of S. cathayensis, building its core germplasm resource bank, and protecting this important endangered plant. Our findings can provide basic reference data for these tasks.

Conclusions
In this study, the MaxEnt model was used to predict the response of the distribution pattern of S. cathayensis to climate change since the LIG, and to simulate the impact of climate change on its future distribution range under two representative concentration paths (RCP2.6 and RCP8.5). The results showed that S. cathayensis may have been formed by the differentiation of Liquidambar formosana and Altingia chinensis from the hybrid origin in Hainan Island to the present flora in South China. In the LIG, the most suitable habitats for the growth of S. cathayensis included the central part of Taiwan Island and the eastern part of Hainan Island, but the distribution range during the LGM was larger than that during the LIG. In the future, except for RCP2.6-2070s, the middle to highly potential suitable of S. cathayensis are predicted to shrink to different degrees, especially under the high emission scenario of RCP8.5-2070s, which is predicted to shrink by nearly 50%. In addition, we identified refugia for S. cathayensis under the future climate. These findings provide an important reference for the development of climate change adaptation strategies for S. cathayensis.