Predicting growth and habitat responses of Ginkgo biloba L. to climate change

We developed a climatic response function using 20-year tree height observed from 45Ginkgo bilobaplantations in China and used it to predict the growth and habitat responses to anticipated climate change. We projected northward and upward shifts in the species habitat and productive areas, but a dramatic contraction of the species distribution is unlikely to occur at least during the present century. Ginkgo biloba is the only living species in the division Ginkgophyta. The species exists in small natural populations in southeastern China but is cultivated across China and the world. The species’ future under climate change is of concern. This study was initiated to model the species’ growth response to climate change and to predict its range of suitable habitat under future climates. Using height data from 45 20 years old plantations growing under a wide range of climatic conditions across China, we developed univariate and bivariate climatic response functions to identify the climate requirements of the species. According to the amount of variance explained (> 70%) and the high level of agreement (> 99%) with independent species occurrence coordinates, the developed climate response function was highly accurate and credible. Projections for future periods (2011–2040, 2041–2070, and 2071–2100) under the Representative Concentration Pathway 4.5 (RCP4.5) scenario indicated that the areas of potential suitable habitat would increase (25–67 million hectares). It would also be associated with northward (0.21–0.62° in latitude) and elevational (24–75 m) shifts. Global climate change is projected to increase the area of potential suitable habitats for Ginkgo and shift its spatial distributions northward and upward.


Introduction
Ginkgo biloba L., commonly known as ginkgo or gingko, is the only living species of the division Ginkgophyta (Leistner and Drewke 2010). The origin of ginkgo can be traced back to approximately 270 million YBP and is believed to have spread worldwide (Hsieh 1992). However, its range has been gradually reduced to small areas in southeastern China during the Pleistocene (Gong et al. 2008;Tang et al. 2012). In the refuge area, the majority of ginkgo trees exhibit multiple secondary trunks and show a strong tendency to asexual reproduction. Ginkgo's morphology is highly conservative and this species exhibits a remarkable genetic stability (Bainian et al. 2003;Major 1967;Royer et al. 2003;Tredici et al. 2010). Due to its high economic, medicinal, ecological, ornamental, and scientific values, the species has been the subject of

Handling editor: Bruno Fady
Contribution of the co-authors Guibing Wang and Tongli Wang conceived and designed the experiments; Ying Guo, Yue Lu, and Lei Feng performed the investigation, analyzed the data, and wrote the manuscript; Tongli Wang and Yousry A. El-Kassaby revised the manuscript.
* Guibing Wang guibinwang99@163.com * Tongli Wang tongli.wang@ubc.ca 1 introduction and cultivation in 32 out of China's 34 provinces, autonomous regions and municipalities (Aziz et al. 2018;Murakami et al. 2012;Ping et al. 2007). Ginkgo has also been introduced to some Asian and European countries, mainly for ornamental and greening purposes (Zhao et al. 2010). However, the significant variation in performance of growth among different introductions has attracted our attention. The major environmental factors driving such a variation and the range of suitable habitat remain unclear, which may limit the success of introduction and cultivation of this species. Additionally, the potential impact of climate change on ginkgo's survival and productivity is also a major concern. A rapidly growing volume of literature shows evidence of climate change and its impact on forest ecosystem and tree species (Franklin et al. 2016). The impact of climate change may even cause a major transformation to terrestrial ecosystems worldwide (Nolan et al. 2018). Climate change often leads to a shift of tree habitat in latitude and elevation, and to some extent is associated with compromised growth performance (Moukrim et al. 2018;Río et al. 2018). The response to climate change varies among tree species. The geographic distribution of the suitable habitat is projected to expand and shift northward for Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco) (Wang et al. 2012) and western Larch (Larix occidentalis Nutt.) in North America (Rehfeldt and Jaquish 2010), while dramatic contractions in habitat would occur for Chinese fir (Cunninghamia lanceolata (Lamb.) Hook), Masson pine (Pinus massoniana Lamb.) , and yellowhorn (Xanthoceras sorbifolium Bunge) in China (Wang et al. 2017). Many of ginkgo's life history characteristics are not conducive to survival in highly disturbed habitats, such as limited shade intolerance, high growth rates, large seeds, long-life span, ability for clonal reproductions, and late sexual maturity (Royer et al. 2003). Thus, modeling the range of suitable habitat and predicting the impact of climate change are particularly critical for this economically and ecologically valuable species. To define the suitable habitat of this species, we focused on the effect of abiotic variables, i.e., climatic or environmental variables.
A suitable habitat can be predicted either by ecological niche models (Pearson and Dawson 2003;Wang et al. 2016) or by climatic response functions (Chakraborty et al. 2018;Hu et al. 2019;Rehfeldt et al. 1999;Wang et al. 2006Wang et al. , 2010. A niche model (a.k.a., ecological niche model (ENM) or bioclimatic envelope model (BEM)) is built based on the relationships between the species occurrences and climatic variables at the corresponding locations. A climatic response function is built based on the relationship between the performance of growth traits and climatic variables from a wide range of climatic conditions under controlled environments (such as common garden experiments), where interactions with other species are eliminated (Chakraborty et al. 2018;Dormann and Singer 2012). Response functions for Douglas-fir built on observations in Austria and part of Germany can predict the performance of this species in Europe (Chakraborty et al. 2016) and its ecological niche in North America (Chakraborty et al. 2018). As most climate response functions are built at population level (Rehfeldt et al. 1999;Wang et al. 2006), those population response functions need to be aggregated to predict at a species level, while this step is not required if a universal response function is developed (Chakraborty et al. 2018;Hu et al. 2019;Wang et al. 2010). The response function can be developed with a relatively small number of samples, but requires tree growth data observed at controlled environments (Campbell and Marco 2004;Fritts 1925). A clear advantage of the response function approach is its capacity to predict growth potential in addition to the suitable habitat of the species if it is built on a growth trait, which has been the case for all response functions up to date.
In recent decades, the large-scale plantation establishment and introductions of ginkgo in China resulted in the availability of a substantial number of ginkgo trees with the same age and growing under different climatic conditions. This provided suitable experimental material to investigate the response of ginkgo to different climatic conditions using the response function approach. Provided the relative uniformity in genetics among ginkgo varieties, we assumed that the variation in growth among genotypes was relatively small and could be treated as random noise in developing response functions. Thus, instead of using observations from common garden experiments as previously done, we used observations on tree height from a set of plantations to build a climatic response function for this particular species. Tree height is considered as one of the ideal growth trait for observation. Dominant height growth is a basic dimension of the ontogenetic development of trees (Falster and Westoby 2003;Ran et al. 2001). Tree height is also an important parameter in ecological genetics to represent adaptability or fitness of a tree species to the environment (Rehfeldt et al. 1999;Ying and Yanchuk 2006). Therefore, the relationship between 20year-old tree height and climate variables is mathematically specified as causal relations although the response function is correlative and empirical in nature (Dormann and Singer 2012). In addition, we obtained ginkgo occurrence samples from the Chinese Digital Herbarium and forest inventory to serve as an independent dataset to validate the range of suitable habitat of this species.
Through building the response functions, we can assess the impact of climate change on this species. Thus, the main objectives of this study were to (1) identify the major climatic factors affecting ginkgo tree height growth; (2) build response function for a major growth trait (i.e., tree height); and (3) predict the shift of the potential distribution of ginkgo in response to global climate change.

Selection of sample locations and data collection
We used ClimateAP software (Wang et al. 2017) to produce climate maps for two climatic variables, annual mean temperature (MAT) and mean annual rainfall (MAP) of China's mainland for the reference period 1961-1990. These two temperature and moisture variables were divided into 10 and 5 grades, respectively, resulting in 50 combinations, which were used as references for choosing sample locations. In addition, we also overlaid a 500 × 500-km latitude and longitude grid (50 grids in total) over the climate maps to consider the geographic distribution of the sample locations.
Through field research and literature review, we found that 20-year-old ginkgo plantations were widely available throughout the regions where this species is introduced. The majority of the fast growth in height of ginkgo is generally achieved in 20 years; and during this period, the vegetative growth of the plants is dominant (Cao 2002). Using the combination of climatic and geographic gradients with the availability of 20-year-old plantations, we identified 15 sample locations located in 14 provinces, autonomous regions and municipalities in China with latitude and longitude spans between 25-45°N and 88-121°E, respectively. We selected three sample plots around each of the 15 locations; thus in total, 45 sample plots were selected as shown in Fig. 1. From each sample plot, 30-50 trees were measured for tree height and the latitude, longitude, and altitude of each plot were recorded. The dataset is accessible at Mendeley Data Repository (Wang 2019).
In addition, we obtained 351 ginkgo occurrence samples from the Chinese Digital Herbarium and forest inventory (the distribution of red dots in Fig. 1). These known occurrence coordinates are located between 24.270°N and 44.174°N in latitude, and between 88.473°E and 123.610°E in longitude (Guo et al. 2019). We used these samples to validate the model predictions.

Climatic variables
ClimateAP software (Wang et al. 2017) was used to obtain climate data for this study. For model development, climate variables for the period 1981-2010 were generated for each specific sample locations based on the latitude, longitude, and altitude of each plot. For spatial mapping, we used 800 × 800 m digital elevation model (DEM) to prepare the input file for ClimateAP to obtain gridded climate data for the current (1981-2010) and future periods. For future climates, we used the ensemble of 15 atmospheric circulation models in the IPCC Fifth Assessment Report (http://www.ipcc.ch/) provided in ClimateAP and selected the Representative Concentration Pathway 4.5 (RCP4.5) scenario, which is closest to the target of the Paris Agreement (Thomson et al. 2011). Climate data were generated for three future periods including 2011-2040 (2020s), 2041-2070 (2050s), and 2071-2100 (2080s) throughout China. The climate variables applied to this study included 14 commonly used climate variables (Table 1).

Univariate quadratic polynomial function
To identify important climate variables that drive the variation in performance of the species across different climatic conditions, we applied a quadratic polynomial regression of tree height growth on each of the 14 climatic variables. The strength of the regression models was compared and the top 10 important climatic variables were identified according to the amount of variance explained by each of the climatic variables. The quadratic polynomial is as follows: where, the dependent variable Y is the average tree height; the independent variable X is a climate variable; and b 0 , b 1 , and b 2 are the parameters to be determined. Results of the regression on the 14 climate variables were compared and an optimal climatic variable was selected to establish the best univariate climate response function.

Bivariate quadratic polynomial function
The following bivariate quadratic polynomial equation was used to build bivariate response models: where Y is the tree height, X 1 and X 2 are the two climatic variables, X 1 X 2 is the interaction term, and b 0 , b 1 , b 2 , b 3 , b 4 , and b 5 are the parameters to be determined. Through comparisons of the regression models on each combination of climate variables, we selected the best combination of two climatic variables to establish the optimal bivariate climate response equation. Due to the limited number of samples, we did not introduce more climatic variables to avoid over fitting.
Both the univariate and bivariate climatic response models were built at the species level. Thus, the genetic variation among all investigated trees (if existing) was considered as a random error in the modeling process. We used optimized bivariate climate response functions to predict potential 20year height across China's mainland for the reference and the three future periods with the raster climate data generated by ClimateAP. If the predicted value was zero or negative, the grid point was considered not suitable for G. biloba growth.

Model evaluation and validation
Models were evaluated based on R 2 value (i.e., the proportion of variance explained), root mean squared error (RMSE) (i.e., the average prediction error), and Akaike information criterion (AIC) (an estimator of the relative quality of models for a given set of data) of the regression model. The model was chosen for its largest R 2 , smallest RMSE, and AIC values, which were consistent in our cases. In addition, we used the 351 ginkgo occurrence samples as an independent dataset to validate the model predictions. The proportion of these occurrence locations in agreement with the model predictions was used as a measure of the reliability of the model predictions.

Univariate climate response models
Using each of the 14 climate variables to build a univariate climate response function for growth trait (20-year-old tree height), we identified the top 10 important climate variables for growth based on their contribution to the model ( Table 2). The model built with mean annual temperature (MAT) was the strongest and explained 62.5% of height total variation of trees growing under different climatic conditions. The second "best model" was the one built with extreme minimum temperature (EMT), explaining 61.3% of the total variation in tree height growth. The parameters and the climate response curve of the best model for tree height are shown in Table 3 and Fig. 2, respectively. As shown in Fig. 2, with the increase of MAT, predicted height showed an increase at the lower temperature end and a decrease at the higher temperature end. The peak values of the curves were around 16°C for tree height. The optimal height growth conditions for ginkgo trees would be between 13 and 18°C (MAT). Annual average temperatures below 11°C or above 20°C would limit its growth. It would be too cold for ginkgo trees to achieve any long-term growth or survive in the case of MAT < 8°C.

Bivariate climate response models
Using a combination of every two of the 14 climatic variables, 91 bivariate climate response models were built for tree height. We excluded the interaction term from the model as it generated predictions far beyond the range of observations. Through comparisons of the model fit, we identified a combination of two climatic variables, DD > 5 (degree-days above 5°C , growing degree-days) and PAS (precipitation as snow (mm) between August in previous year and July in current year), for building the "best bivariate climate response models" for tree height growth. This model displays a bellshaped response surface (Fig. 3, left panel) and explained 70% (R 2 = 0.70, P < 0.0001) of the total variation in tree growth with prediction error (RMSE) of 1.2 m and AIC of 152.4, which is a much better model fit than that of the best univariate model (Table 3); thus, it was used for current and future predictions. The model parameters are listed in Table 4. The response surface for tree height growth peaked at the climatic conditions between 4200~5600°C-days and 10~25 mm for DD > 5°C and in PAS, respectively. Under these climatic conditions, tree height was predicted to reach over 10 m at the age of 20.

Predicted growth potential in China
After the response surface was mapped over the country (Fig.  3, right panel), we found that except for parts of the Northeast (such as Heilongjiang, Jilin, and Liaoning Province) and Southwest China (such as Qinghai Province and Tibet Autonomous Region), the climate conditions in most areas of China were suitable for ginkgo growth (20-year-old tree height > 0 m). The climatic conditions in the central and   eastern regions of China (such as Shandong, Jiangsu, Anhui, and Hubei Province) were predicted to be more conducive to tree height growth. We found a strong agreement between the distributions of the species occurrence samples and the areas of potential distribution (Fig. 3, right panel). All but three of the 351 ginkgo occurrence samples (one in Sichuan province and two in Northeast provinces) fell into areas where the species was predicted to have a height growth potential, and the coincidence was 99.2%.

Future potential distribution of ginkgo
Tree height was projected using the "best" bivariate climate response model in Table 4 for three future periods 2020s, 2050s, and 2080s under the RCP4.5 scenario (Fig. 4). Under the future climatic conditions, the range of the suitable habitat showed a trend of shifting northward (0.21-0.62°in latitude) and elevational upward (24-75 m) with increases in total areas (25-67 million hectares) (Table 5). Geographic expansions of the suitable habitat were projected to occur mainly in the central part of Xinjiang Uygur Autonomous Region, the western part of Inner Mongolia Autonomous Region, and parts of Ningxia. Expansions would also occur in the southern part of Liaoning province, Heilongjiang Province, and northeastern Inner Mongolia Autonomous Region.

Discussion
It is of great importance to define the range of suitable climatic conditions and to assess the impact of climate change for an iconic tree species such as Ginkgo biloba. In this study, we identified the top 10 most important climatic variables affecting its height at age 20 using univariate climatic response functions and defined the range of suitable climatic conditions (Fig. 4) using bivariate climatic response functions.

Prediction of suitable habitat using a response function
In this study, we built a climatic response function for 20-year tree height based on observations from a set of commercial  The high accuracy of our model indicates a small random error without considering genetic effect, confirming our assumption of a relatively small genetic variation in growth among varieties in this species. With such a model, we could not only predict the range of the suitable habitat, but also the height growth potential within its suitable habitat, suggesting a clear advantage of the climatic response model over the conventional ecological niche models.

Climate variables contributing to the fundamental niche models
It is worth noting that the climatic factor with the largest contribution to growth revealed by the univariate models is MAT, whereas the two climatic variables in the "best bivariate models" were DD > 5 and PAS. In general, DD > 5 is strongly correlated with MAT and specifically in the present study (r 2 = 0.98, P < 0.0001). Still, in combination with PAS, DD > 5 improved the model over MAT (r 2 = 0.70 vs. r 2 = 0.65). Nevertheless, both MAT and DD > 5 reflect heat or temperature status of a site. Temperature is considered a key factor in controlling tree growth in general (Körner 1994). It has been found that for every 1°C increase in MAT, the length of the growing season of ginkgo may increase by 10 days (Matsumoto et al. 2010). Our finding is in agreement with other studies suggesting that temperature plays the most important role in affecting tree growth of many forest tree species (Rehfeldt et al. 1999;Wang et al. 2010). The second climatic variable in the model, PAS, was not among the top 10 important climatic variables revealed by the univariate models. PAS is an integral variable reflecting a combination between temperature and precipitation because temperature affects the proportion of precipitation as snow. The optimal PAS values for ginkgo were between 15 and 25 mm, suggesting that that the regions with relatively cool and dry (small amount of snowfall) weather favor ginkgo growth, which is in agreement with productive areas of this species. Thus, the combination of DD > 5 and PAS could be somehow considered as a proxy of temperaturemoisture combination.

Future projections and management implications
Under future climatic conditions, we projected that the range of the ginkgo would shift northward and upward. Many studies support this finding, speculating that plants will move to higher altitudes and high latitudes due to global warming (Costa et al. 2017;Innes 1991;Winkel et al. 2015). Studies have also shown that temperate and  cold temperate trees may grow faster in moderate warming conditions (Davis and Shaw 2001;Norby and Luo 2004;Peltola et al. 2002;Ruiz-Labourdette et al. 2015). For example, the seed resources of Pinus banksiana Lamb. in Canada and northern Europe are growing below the optimal temperature range and will benefit from future climate increases (Savva et al. 2007). Our findings entail a common concern that deciduous species may grow in the north as the climate warms, and will gradually develop into a true northern coniferous forest in future climatic conditions. On the other hand, local population extirpation of ginkgo may occur in the south at low elevations under climate change, and monitoring the sign of maladaptation should be considered in these areas. However, ex situ conservation may not be necessary as all populations originate from small refuge areas. Overall, we believe that the fate of this species is not under drastic threat and a dramatic contraction in its distribution caused by climate change is not anticipated at least within this century based on moderate climate change scenarios.

Conclusion
In this study, we found that temperature is the dominant climate variable affecting tree height at age 20 and the range of suitable habitat of ginkgo. The most favorable climatic conditions for this species are located in the central and eastern regions of China. We also found that areas with suitable climatic conditions for this species would increase, and shift northward and upward. We believe that the climatic response model and its outputs presented in this study can serve as an effective tool for assessing the growth potential and future distributions of ginkgo under current or future climate conditions.