Current and future distribution of the deciduous shrub Hydrangea macrophylla in China estimated by MaxEnt

Abstract Climate change has a significant impact on the growth and distribution of vegetation worldwide. Hydrangea macrophylla is widely distributed and considered a model species for studying the distribution and responses of shrub plants under climate change. These results can inform decision‐making regarding shrub plant protection, management, and introduction of germplasm resources, and are of great importance for formulating ecological countermeasures to climate change in the future. We used the maximum entropy model to predict the change, scope expansion/reduction, centroid movement, and dominant climate factors that restrict the growth and distribution of H. macrophylla in China under current and future climate change scenarios. It was found that both precipitation and temperature affect the distribution of suitable habitat for H. macrophylla. Akaike information criterion (AICc) was used to select the feature combination (FC) and the regularization multiplier (RM). After the establishment of the optimal model (FC = QP, RM = 0.5), the complexity and over‐fitting degree of the model were low (delta AICc = 0, omission rate = 0.026, difference between training and testing area under the curve values = 0.0009), indicating that it had high accuracy in predicting the potential geographical distribution of H. macrophylla (area under the curve = 0.979). Overall, from the current period to future, the potential suitable habitat of this species in China expanded to the north. The greenhouse effect caused by an increase in CO2 emissions would not only increase the area of high‐suitability habitat in Central China, but also expand the area of total suitable habitat in the north. Under the maximum greenhouse gas emission scenario (RCP8.5), the migration distance of the centroid was the longest (e.g., By 2070s, the centroids of total and highly suitable areas have shifted 186.15 km and 89.84 km, respectively).


| INTRODUC TI ON
Global climate change affects many ecosystems and biota. It not only affects the growth and distribution of global vegetation, but also leads to the loss of species diversity and germplasm resources (Kozak et al., 2008). By 2100, global temperatures are predicted to increase by 1.4 to 5.8°C from 1900, two-to ten-fold higher than that in the 20th century, and the average temperature in China will rise 2.2°C by 2050 (Kumar, 2007). According to the fifth assessment report of the Intergovernmental Panel on Climate Change, global average temperature will increase 3.7°C by the end of the 21st century under the scenarios of Representative Concentration Paths (RCPs) (Lee et al., 2020). Climate variables affect vegetation composition and species distribution, and these, in turn, can directly reflect global climate change (Ni & Song, 1997). Many common species, including grasses, hybrids, and shrubs, are highly sensitive to climate change (Kane et al., 2017). The prediction of the potential geographical distribution of species based on the biological-climate relationship can not only lay the foundation for theoretical research on species origin, vegetation division, and floristic formation, but also play a significant role in genetic improvement, and species introduction and domestication . It is also very important for analyzing the potential distribution of species, and for planning species protection and sustainable resource use strategies. Understanding the spatial patterns of species and their dependence on the environment is one of the basic goals in ecology and evolution (Merow & Silander, 2014).
Climate models can determine the geographical distribution of biological ecosystems, and the species range is usually defined by their bioclimatic range (Christopher & Elsa, 2014). The change in species distribution range shows a polar variation in latitude and altitude to track the climate niche, which is the commonly expected result of global warming (Katherine & Monique, 2015;Laura et al., 2014), especially at high latitudes and elevations. The maximum entropy (MaxEnt) model is the most widely used tool to study species distribution and, in particular, to predict the potential future geographical distribution of species, need for only species occurrence data and environmental information (Jane et al., 2011;Phillips et al., 2006). Current research has attempted to predict the impact of future climate change on the distribution of shrub species (Bhandari et al., 2020;Xu et al., 2019). Shrubs play an important role in biogeochemical cycles, prevent soil and water erosion, provide forage for livestock, and are a source of food, wood, and non-wood products (Hageer et al., 2017). Hydrangea macrophylla (Figure 1a) is a semi-cold tolerant deciduous shrub found in the warm temperate zone (Wu et al., 2020). As a woody shrub plant, it is an important component of undergrowth in many areas and has been widely introduced globally as a horticultural plant. Presently, it is an important planting material for contemporary urban forestry construction (Galopin et al., 2010).
It is named after its inflorescence that is very similar to the embroidered ball in traditional Chinese culture; its inflorescence is large, rich in color, and has high ornamental value. It is known as one of the flowers with great development potential (Zeng et al., 2018). In summer, even growing in the shady environments of Nanjing, Jiangsu Province, the leaves of plants are prone to wilting due to long-term arid conditions (Figure 1b). In winter, the process of defoliation is not completed until the middle of December (Figure 1c). Its natural distribution areas, and introduction and cultivation areas have spread globally, including in the temperate zone, warm temperate zone, north subtropical zone, and south subtropical zone (Zeng et al., 2018). Therefore, H. macrophylla is an ideal model species for studying the distribution of plants and their response to climate change.
The Hydrangea genus in China is rich and diverse, and the country is considered the center of Hydrangea germplasm resources. Most of the parent species of the cultivated species are distributed in China, and the cultivation history is very long. According to literature, the cultivation research of H. macrophylla in China could be traced back to the Tang Dynasty, and was introduced to England in the 18th century (Ceng et al., 2018;Qiao et al., 2020). Records show that there are approximately 73 species of H. macrophylla globally, with 47 species and 11 varieties in China (Cai et al., 2018). At present, there has been some research conducted on the response of H. macrophylla to high and low temperature stress and drought (Cai et al., 2018;Zhang, Li, et al., 2018;Zhang, Meng-qi, et al., 2018), but its distribution, limiting factors, and possible response to future climate change in China have not been reported. This study was the first attempt investigating the spatial distribution of this plant species in China as a whole. These are important basic data for the development of H. macrophylla planting, scenic spot planning, and layout division, and therefore, it is important to collect such data considering the abundant germplasm resources of H. macrophylla in China.
In combination with the literature and sample collection records of H. macrophylla, the MaxEnt model was used to predict the potential distribution of this species under three different RCPs (RCP2.6, RCP4.5, andRCP8.5) under current (1960-1990) and future conditions: 2050s (2041-2060) and 2070s (2061-2080), and to explore the key environmental factors causing potential geographical distribution changes. The purpose of this study was to further understand the species distribution characteristics and ecological adaptability of small deciduous shrubs, and to provide theoretical reference for the protection and utilization of the wild resources of H. macrophylla, as well as future large-scale introduction and cultivation (Qin et al., 2020 (Coban et al., 2020) are used for future climate data. The RCP scenarios consist of four pathways, including RCP2.6, RCP4.5, RCP6.0, and RCP8.5 , Remya et al., 2015. RCP4.5 and RCP6.0 are both moderate greenhouse gas emission scenarios, and RCP4.5 has a higher priority than RCP6.0 .

| Selection of environmental variables
In order to avoid correlation between various climate variables, which leads to over-fitting and affects the accuracy of prediction, 19 bioclimatic variables data from 307 data points were first extracted using QGIS v.3.14, and then imported into R v.3.6.3. Pearson correlation analysis of data was carried out using the cor function in R Wang, Wu, et al., 2018). For highly correlated variables (correlation coefficient >0.8), only one of the more suitable variables for model interpretation was selected as the representative variable to participate in subsequent model prediction Radosavljevic & Anderson, 2014). Finally, 10 biological variables were selected from 19 biological variables to participate in the model prediction (Table 1).

| Optimization of model parameters and model building
The predictive performance of the model is influenced by two parameters that require optimization: (i) the regularization multiplier (RM) and (ii) the feature combination (FC). Optimization was achieved using the kuenm package (Cobos et al., 2019) in R v3.6.3.
In the optimization process, 307 distribution records were divided into four equal parts by block method, including three for training and one for testing. Then, we set RM constants 0.5 to 4.0 with steps of 0.5 . For the selection of FC, we selected 29 combinations of five feature classes (linear = l, quadratic = q, product = p, threshold = t, and hinge = h) for testing, including L, LQ, LQH, and LQHPT, etc (Appendix 1). Finally, the above 232 candidate models were tested by kuenm package. Best models were selected according to the following criteria: (1) significant models with (2) omission rates ≤5% (Cobos et al., 2019). Then, the model with the minimum AICc value (i.e., delta AICc = 0) was considered as the final model, and all models with delta AICc <2 had high reliability Wang, Wu, et al., 2018).
Then, we saved 307 records in.csv format initially, and then reflecting the sensitivity and specificity of a model (Gebrewahid et al., 2020). In practical applications, the AUC index derived from ROC is commonly used to evaluate the model performance index . Since the AUC value is unaffected by the threshold value, it is a more reasonable basis for comparing different models. When the AUC value exceeds 0.8, the simulation model is accepted. For a better simulation model, this value should exceed 0.9 Sun et al., 2020;Zhang, Li, et al., 2018;Zhang, Meng-qi, et al., 2018).

| Classification of suitable habitat
The data generated by MaxEnt software were imported into ArcGIS v.10.6. The conversion tool was used to convert the data into raster data, and then the classification of suitable H. macrophylla habitat was carried out using the reclassify tool (Hu & He, 2020). According to the method of artificial classification, the distribution areas of H. macrophylla were classified as follows: unsuitable area (0 ≤ p ≤ .1), low-suitability area (.1 < p ≤ .3), medium-suitability area (.3 < p ≤ .5), and high-suitability area (.5 < p ≤ 1) (Qiu et al., 2020). Finally, we obtained the potential distribution of suitable habitat for H. macrophylla in China (p > .1).

| Changes in suitable habitat area and centroids
To further examine trends, the centroids of the current and future climate distribution areas were calculated using the Python-based GIS toolkit, SDMtoolbox (Brown et al., 2017). We imported.asc

| Optimal model and accuracy evaluation
The Maxent model was used to simulate the potential distribution area of H. macrophylla in China. When the model was set with the default parameters, delta AICc was 126.56, but when it was set with the optimization parameters (FC = QP and RM = 0.5), delta AICc was 0 (

| Future changes in suitable habitat area
Under different future climate scenarios, four grades of the potential geographical distribution areas vary significantly, but the change in the total suitable area is less than that of current times. As shown in Figures 3 and 4, in the 2050s, under three different emission scenarios (RCP2.6, RCP4.5, and RCP8.5), the area of unsuitable habitat will reduce by 2%, 3.67%, and 6.33%, respectively, compared with that of current times. In the 2070s, the model predicts that the area of unsuitable habitat of H. macrophylla will be 3.04% (RCP2.6), 4.86% (RCP4.5), and 7.41% (RCP8.5), less than it is currently. This indicates that many unsuitable areas will be replaced by low-, medium-, or high-suitability areas. Therefore, the total suitable area of H. macrophylla will increase in the future. The area of low-suitability habitat increases by 0.87% (RCP2.6), 1.70% (RCP4.5), and 3.32% (RCP8.5) in the 2050s. This further increases by 0.58% (RCP2.6), 1.11% (RCP4.6), and 0.72% (RCP8.5) in the 2070s. From current times to the 2070s, the low-suitability habitat of H. macrophylla gradually increases, and the range increases gradually with rising greenhouse gas concentration. Under the three different emission scenarios, the area of medium-suitability habitat decreases by 2.27%, 3.83%, and 3.02%, respectively, from current times to the 2050s. From current times to 2070s, it decreases by 3.11%, 3.46%, 2.15%, respectively.
Some of these suitable areas will be transformed into low-suitability habitat and some into high-suitability habitat. The high-suitability habitat area of H. macrophylla also shows a trend of expansion in the future. By the 2050s, the high-suitability habitat area increases by 3.40%, 5.79%, and 6.03%, respectively. Compared with the current area, the area of these regions will increase by 4.70%, 5.50%, and 5.51%, respectively, in the 2070s. In general, the change of suitable habitat of H. macrophylla in the future shows a trend of decreasing non-suitable habitat and increasing suitable habitat increasing under three different greenhouse gas emission scenarios. At the same time, low-and high-suitability habitat will increase, and mediumsuitability habitat will begin to differentiate.

| Analysis of spatial pattern changes
We

| Environmental variables analysis
The results of jackknife test showed that the top three climate variables with the largest influence on the normalized training gain, test gain, and AUC value were as follows: precipitation of the warmest quarter (Bio18), min temperature of the coldest month (Bio6), and annual precipitation (Bio12), when modelling with a single environmental factor (Appendix 4). Furthermore, the contribution rate of environmental variables can also be used to evaluate the importance of the environment. The contribution rate of each environmental variable to the distribution of H. macrophylla was combined (Table 3)

| Predictive capabilities of the MaxEnt model
The MaxEnt model has been widely used in ecology, conservation biology, evolutionary biology, and invasive species management (Phillips et al., 2006). It has many advantages, such as wide application range, high precision, simple operation, low sample number requirement, and stable operation results .
Hydrangea macrophylla is the most popular flowering shrub that is widely used in landscape gardening and has crucial commercial value (Lisa, 2019). As a globally introduced plant species, a less studied but important characteristic is its adaptability to the environment. The key environmental factors affecting or limiting its growth should be studied by assessing its species distribution patterns.
Based on the results of the jackknife test and the contribution rate of climate variables calculated by the model, we determined the main environmental factors affecting the variation in H. macrophylla distribution. The predicted results showed that the variables involved in the regional changes of suitable habitat for H. macrophylla were as follows: min temperature of the coldest month (Bio6), precipitation of the warmest quarter (Bio18), annual precipitation (Bio12), and temperature seasonality (Bio4). Therefore, its distribution was determined by both temperature and precipitation factors.
Both high temperature and low temperature stress could affect its growth. When H. macrophylla was stressed at high temperature, a series of changes occurred in its internal physiological mechanism, which mainly showed chlorophyll degradation, increased relative permeability of cell membrane, and increased MDA content (Xin & Shi, 2009 (Pagter et al., 2008). H. macrophylla is a very popular flowering plant, low temperature or frozen soil will seriously damage the buds, thus considerably reducing the attractiveness of H. macrophylla to consumers (Ren et al., 2020). In addition, precipitation is one of the most important environmental factors affecting plant growth, and it is an important factor in the survival and growth of seedlings (Liu, 2020). Abundant rainfall leads to a significant increase in soil water content near the planting site, providing an excellent growth environment for H. macrophylla, which has a large water demand. However, if the rainfall is too much, the waterlogging stress will lead to the closure of stomata, the decrease of photosynthetic capacity and the increase of respiratory energy consumption, which is not conducive to the accumulation of organic matter. But some varieties of H. macrophylla had strong adaptability to flooding stress, and could be planted in waterfront, low-lying wetland, and other sites .
In the arid environment, the leaf transpiration of H. macrophylla was high, and even a short period of wilting would lead to the visual phenomena of leaf edge drying and flower necrosis, which would lead to the decrease of yield, quality, or ornamental value (Ren et al., 2020). The harm degree of environmental factors to H. macrophylla was also related to the length of time it was under stress. When H. macrophylla is under high temperature or drought stress in the early-life stages, it can resist stress damage by adjusting its osmotic regulation substances and protective enzyme activities; however, when the stress duration increases by a certain extent, the damage caused is irreversible. The damage of double stress is far greater F I G U R E 5 Comparison of changes in the spatial pattern of total suitable Hydrangea macrophylla habitat under three different climate scenarios (RCP2.6, RCP4.5, and RCP8.5) in the 2050s and 2070s. Total suitable area change: spruce green represents expansion area, red represents contraction area, gray represents unchanged area, and white represents non-suitable area than that of single stress . All of these stress phe-

| Changes in spatial pattern and centroid of potential Hydrangea macrophylla distribution
The overall change in suitable habitat area of H. macrophylla was as follows: from the current time to 2070s, the total suitable habitat area would continue to expand with the increase of greenhouse gas emission. It was consistent with the results of a study on the Chinese Sea Buckthorn (Hippophae rhamnoides subsp. sinensis) suitable area, which the planting potential of this shrub in the high emission scenario was greater than that in the low emission scenario (Huang et al., 2018). A study on the distribution of the alpine vegetation distribution in the Greater Himalaya also showed that under the high concentration emission scenario RCP 8.5, the low-and high-suitability habitat for shrub species would be increased in the middle of the 21st century with the increase of CO 2 concentration (Bashir, 2021). A similar expansion was also expected in the western Himalayas under future global warming (Changjun et al., 2021). The reason why the distribution area would expand in the future is that the climate in China would transition to warm and humid under the influence of the greenhouse effect induced by the increase of CO 2 concentration. However, the change of the medium-suitability areas F I G U R E 6 Spatial pattern of the Hydrangea macrophylla high-suitability habitat area in the 2050s and 2070s under three different climate scenarios (RCP2.6, RCP4.5, and RCP8.5). High-suitability area change: grass green represents expansion area, purple represents contraction area, gray represents unchanged area, and white represents non-suitable area of H. macrophylla was not completely consistent with the change of the above area. to high-altitude and high-latitude areas (Ashraf et al., 2016;Wardle & Coleman, 2011;Wu et al., 2021). Only under RCP2.6-2070s and RCP4.5-2050s, the centroid of the highly suitability habitat moved to the southeast. But under the scenarios of other periods, it was consistent with that of the total suitable area and moved to the north. Combined the trajectory distance and direction of centroid in different periods, we speculated that the centroid would not always move in the same direction, and the smaller the time span, the more complex its trajectory is, which was similar to a previous study on other plant Wang, Wu, et al., 2018).

| Study limitations
The first limitation of this study was that the 10 environmental variables selected cannot reflect all the factors affecting the geographical distribution of H. macrophylla. In addition to the selected factors, other abiotic factors, such as light, soil, and air; biological factors, such as the interaction between species; and anthropogenic impacts all have a significant impact on the species distribution (Qin et al., 2020). Secondly, most species distribution models speculate the environmental conditions of species based on the observed species occurrence data. Therefore, we cannot simply equate the habitat of species occurrence data with the suitable habitat for a species, that is, the suitable habitat area predicted by the model does not necessarily contain the species (Li et al., 2013).
Uncertainties such as the internal assumptions of species distribution models need to be considered and further studied.

ACK N OWLED G M ENTS
The authors would like to thank the editors, as well as the anonymous reviewers, for their valuable comments. We thank International Science Editing (http://www.inter natio nalsc ience editi ng.com) for editing this manuscript.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

A PPE N D I X 4
Relative predictive power of different environmental variables based on the jackknife of AUC for Hydrangea macrophylla. The gain (%) with all predictors = Red bar; the gain (%) with only predictor = Blue bar; the gain (%) without the corresponding predictor = Green bar