Radial growth response of two oaks to climate at their disparate distribution limits in semiarid areas, Beijing, China

Impending climate warming is expected to influence plant growth and distribution, and the distribution range limit of tree species is extremely sensitive to climate change. However, synchronous comparisons of responses of different tree species with overlapping ecological niches to climate at their distribution range limit in the same region have received little attention. In the present study, we assessed the discrepancy in radial growth responses of two dominant oak species to climate at their distribution range limit (southern range limit for Quercus mongolica, northern range limit for Quercus variabilis) in Beijing, China. Furthermore, growth–climate relationships were examined using linear mixed-effects models, and growth trends up to the year 2100 were forecast based on future climate scenarios (RCP2.6 and RCP8.5). Our results indicated that there were no significant differences in the growth response to climate between Q. variabilis at its northern limit and Q. mongolica at its southern limit, and the growth of the two tree species was positively correlated with the minimum temperature and negatively correlated with the climatic moisture deficit of autumn. However, the growth response to climate variables varied at different sites. Additionally, the model forecast showed an increase in radial growth of Q. variabilis at its northern limit and Q. mongolica at its southern limit up to the end of the present century.


INTRODUCTION
Climate (temperature and precipitation) is the critical driver of plant distribution at the global scale (Woodward 1987, Walck 2011, because it strongly influences growth and survival at the population scale (Arft 1999, Diemer 2002, Walck 2011. Impending climate change is expected to influence plant growth and regeneration, with consequent effects on species range shifts, especially at their distribution range limit. Therefore, there is an urgent need to move beyond reporting changes already underway to increasing our predictive capacity, enabling us to better estimate and model the ecological consequences of climate change for species under future climate change scenarios (Steinkamp et al. 2015).
Previous studies demonstrated that the growth of tree species at the northern edge of their range was mainly affected by low temperature (Babst 2013, S anchez-Salguero et al. 2015, Gao et al. 2018), but was restricted by the high temperatures and warming-induced water shortage at the southern limit (S anchez-Salguero et al. 2015, Gao et al. 2018. Most of these studies were based on the northern range limits in alpine timberline/high-latitude areas (Gao et al. 2018) or southern limits in Mediterranean areas (Camarero 2015, Arzac 2016, S anchez-Salguero et al. 2017. However, in continental climatic regions (such as semiarid area of temperate), whether the growth of tree species was also limited by water shortage at their southern edges has been less investigated. In addition, the growth of tree species with different distribution ranges is likely to respond differently to climate change. For example, Huang (2010) examined four boreal species in eastern Canada and found inconsistent responsiveness for these species to recent warming; that is, the radial growth of Populus tremuloides was mainly affected by water conditions during the year, whereas temperature in January and precipitation during the growing season positively impacted the growth of Betula papyrifera, but the growth of Picea mariana and Pinus banksiana was positively affected by the temperatures during the current-year winter and spring or over the entire growing season. However, Tardif et al. (2010) found that growth of Quercus alba and Quercus rubra exhibited no significant differences in responses to climate at their northern range limit. Therefore, the growth response among different tree species to climate at distribution range limits is complex, and there is a lack of systematic research on the responses of different tree species with overlapping ecological niches to climate in the same region.
Climate warming can accelerate the decomposition rate of soil organic matter, thus increasing plant soil mineral nutrition uptake, especially that of the most common limiting nutrient, N (Luo 2004, Huang 2010. Furthermore, rising atmospheric CO 2 concentrations and associated climate warming can also enhance plant photosynthesis and water-use efficiency, thus increasing biomass accumulation of above-and belowground biomass (K€ orner 2006(K€ orner , Huang 2007. Together, these effects might lead to a shift in species growth and distribution (Huang 2010). Previous studies demonstrated that temperate and boreal tree species presented a growth increase toward the northern parts of the species range limits in the Northern Hemisphere (Jump 2007, S anchez-Salguero et al. 2015, Gao et al. 2017. However, there are still research gaps on how rising temperatures will affect tree growth under projected climatic conditions, especially the southern range limit in semiarid area. Therefore, within the context of sustainable forest management, prediction of the potential growth trends of forest under climate change is important to reveal the distribution range dynamics. The Mongolian oak (Quercus mongolica Fisch) and Chinese cork oak (Quercus variabilis Blume) are two of the most widespread and abundant deciduous broadleaf tree species in China, ranging from~34°to 53°N, 112°to 134°E and 19°to 42°N, 97°to 140°E, respectively. They provide not only supplies for raw material for cork and timber products, but ecological services (e.g., carbon sequestration, and water and soil conservation) as well. According to species distribution models and field investigations, our previous work (Gao et al. 2017(Gao et al. , 2018, as well as that of Li and Zhang (2015), revealed that the distribution range of the Chinese cork oak may shrink and retreat under the expected climate warming. However, the growth response at their range limit is not yet understood. Therefore, in the present study, we assessed the climate response of radial growth of these two tree species with different distribution range limits in the same area (Beijing, China, is the southern range limit for Q. mongolica and northern range limit for Q. variabilis). Furthermore, we predicted the future growth trends under projected climate warming. In particular, we considered the representative concentration pathways (RCPs), which are climatic projections consistent with a wide range of possible changes under future anthropogenic greenhouse gas emission scenarios (Collins 2014). Our specific objectives were to (1) characterize growth responses to climate at the southern range limit for Q. mongolica and northern range limit for Q. variabilis and (2) predict the future growth trends of forests under different climatic force scenarios.

Study site and stand selection
This work was conducted in Beijing, China, which is the southern range limit for Q. mongolica and northern range limit for Q. variabilis (Appendix S1: Fig. S1). The climate is the typical warm-temperate continental monsoon climate with hot and rainy summers and cold and dry winters. Based on climate data from 1950 to 2015, the annual mean temperature of the study region varies from 5.5°C to 10°C, and annual total precipitation ranges from 570 to 680 mm (Table 1). The mean annual temperature at all sites has increased significantly over time, whereas the annual total precipitation has significantly decreased (Appendix S1: Fig. S2). These climate changes have the potential to affect forest ecosystems. A total of 3 sites along an altitudinal gradient ranging from 200 to 500 m for Q. variabilis (sites: SPL-CP, SPL-PGA, and SPL-PGB) and from 300 to 1100 m for Q. mongolica (sites: MGL-CP, MGL-PG, and MGL-HR) were selected (Table 1).

Field investigation
Fieldwork was conducted from August to September in 2018. At each site, at least three temporary plots of 400 m 2 (20 9 20 m) were established. The criteria for plot selection were as follows: (1) Q. mongolica and Q. variabilis were the dominant species in the overstory; (2) the withinstand conditions were relatively homogeneous with little anthropogenic disturbance; and (3) the position of the plots was randomly arranged within a stand if the previous criteria were satisfied. In total, 30 plots were investigated. Topographic factors (latitude, longitude, elevation, aspect, and slope) were recorded (Table 1). Ten dominant trees per population in each plot were selected for dendrochronological study to remove the competition effect. Sampling was performed following standard dendrochronological methods (Fritts 1976). For each sampled tree, two cores were taken using a 4.3 mm increment borer and diameter at breast height was recorded (dbh; 1.3 m above the ground). A total of 302 cores were collected (Table 1). Samples were then stored and air-dried in paper straws. The wood core samples were mounted, polished with successively finer grit sandpaper until rings were clearly visible, scanned at 1200 dpi, and then measured to an accuracy of 0.01 mm using a WinD-ENDRO (Regent Instruments, Canada). All ringwidth series were corrected via cross-dating, performed by both visual check and COFECHA for quality control of dating (Holmes 1983). To control for the geometric trends of decreasing ring width with increasing tree size and age, the ring width was converted into basal area increments (BAIs) which could remove the impact of tree size and age on annual growth trends while keeping the high-and low-frequency signals from treering series (Biondi and Qeadan 2008): where R is the radius of the tree and n is the year of tree-ring formation.

Climate data
The climate data of the historical (1950-2015) and future (2016-2100) periods for all sites were obtained from ClimateAP v2.20 (http://asiapacif ic.forestry.ubc.ca/research-approaches/climate-mod eling/). ClimateAP is an application for dynamic local downscaling of historical and future climate data in the Asia-Pacific area, which generated scale-free (i.e., specific to point locations) climate Note: H, tree height; T, annual mean temperature; P, annual precipitation; and CMD, climatic moisture deficit. data for historical  and future (2016-2100) years and periods (Wang 2017

Statistical analysis
Climate-growth correlation analysis.-Mean BAI chronologies for each site were computed over all trees for the same year (Jump et al. 2006). The relationships between climate variables and the BAI were investigated based on Pearson's correlation coefficient. The BAI chronologies were correlated separately for monthly climate (T min , T max , T ave , and P), seasonal variables (T min , T max , T ave , P, PAS, and CMD), vegetation period (T min , T max , T ave , P, PAS, and CMD), and physiological year (T min , T max , T ave , P, PAS, and CMD) for every species and all sites. To determine spatial variation of growth sensitivity to climate among the sites for the two species, a principal component analysis (PCA) was performed based on the climate-growth correlation coefficients and then plot the sites in the PCA, following the method from Merian et al. (2011).
Tree radial growth prediction under future climate change. -Linear mixed-effects models (LMM) were used to analyze the effects of climate variables (including vegetation period (March-September) and physiological year (previous October-current September)) on growth (Jump et al. 2006, Matias et al. 2017. Climate was treated as a fixed effect, and plots were treated as random effects to show the nested nature of our experimental design. Accordingly, the model was specified as follows (Matias et al. 2017): where BAI ij is the mean annual basal area increment of tree j at plot i; b 0 is the estimated fixed intercept; b i0 represents the random effects for plots; b 1 is the model coefficient vector estimated for the climate variables; climate represents the candidate vegetation period and physiological year variables; and e ij is the residual error.
The models were fitted based on the maximum likelihood and were compared according to the corrected Akaike information criterion (AIC c ; Schwarz 1978), which is unbiased for small samples. We used the package lme4 (Bates et al. 2014) in R (R Core Team 2013) to fit the models. The significance of each source of variation included in these models was tested using F tests (Zuur et al. 2009). The total variance explained by the fixed effects alone (marginal R 2 ) and by the fixed and random structure of the models (conditional R 2 ) was calculated according to Nakagawa and Schielzeth (2013). The selected models were conducted to forecast annual BAI of tree species at each site for the period 2017-2100.

Observed past growth
We found a consistent growth trend for the populations of Q. variabilis and Q. mongolica during the growth period, but it varied among the different sites (Table 1; Fig. 1). All populations of Q. mongolica and Q. variabilis exhibit an increasing growth trend until 2015, and the slope was highest at PG and lowest at CP (MGL-CP R 2 = 0.06, P < 0.0001, slope = 4.14; MGL-PG R 2 = 0.34, P < 0.0001, slope = 47.29; MGL-HR R 2 = 0.200, P < 0.0001, slope = 16.61; SPL-CP R 2 = 0.06, P < 0.0001, slope = 4.29; SPL-PGA R 2 = 0.21, P < 0.0001, slope = 20.46; SPL-PGB R 2 = 0.40, P < 0.0001, slope = 34.00; Table 2 Responses of annual growth to monthly and seasonal climate For both tree species, similar patterns of growth response to monthly temperature and precipitation were identified, although there was some variation among different sites (Figs. 3, 4). The BAI of the two tree species was positively correlated with monthly temperature, especially the minimum temperature, but negatively correlated with monthly precipitation of the previous and current August (Figs. 3,4). Similar results were found for the seasonal climate response for growth of the two tree species, which were positively correlated with seasonal temperature, but negatively correlated with precipitation during summer (Figs. 5,6). In addition, the PCA for correlation coefficients between seasonal climate and growth among different populations showed that the first axis (PC1) represented 43.5% of the total inertia; the scores of the chronologies (e.g., in MGL-CP and SPL-CP) for the two species in the same site were similar, but varied among different sites, although the climate-growth correlation was not correlated with altitude (Fig. 7a). The main structuring climate variables were related to temperature, precipitation of summer, and CMD of autumn, with positive scores for temperature and summer precipitation and negative scores for autumn CMD (Fig. 7b). The second axis represented 31.6% of the total inertia and was also not correlated with altitude, and the primary distinguished climate variables were precipitation and temperature (Fig. 7a, b). The growth of Q. mongolica was negatively related to precipitation of summer, vegetation period and physiological year, and precipitation as snow between August in the previous year and July in the current year (PAS) of winter, and vegetation period in CP, but positively correlated with autumn precipitation and negatively correlated with PAS of summer in PG. The growth of Q. variabilis was negatively Note: ΔT, mean annual temperature change (°C); ΔP, mean annual precipitation change (mm); ΔCMD, climatic moisture deficit change (mm). CMD, climatic moisture deficit.   correlated with precipitation of summer, CMD of autumn, and PAS of spring and physiological year in PGA, but positively correlated with autumn precipitation in the three sites. The growth of MGL-HR and SPL-PGB populations was negatively correlated with CMD of autumn, vegetation period, and physiological year, but positively correlated with autumn precipitation.
The LMM results indicated that the growth of both Q. variabilis and Q. mongolica was mainly affected by the mean temperature of vegetation period (T ave _VP) and physiological year (T ave _PY), and secondly by the precipitation and climatic moisture deficit of vegetation period and physiological year (CMD_VP/PY; Table 3). The growth of both tree species was significantly positively related to T ave of physiological year/T ave of vegetation period and negatively related to precipitation and CMD_VP/PY, except in PGB for Q. variabilis. In addition, the growth Q. mongolica in HR exhibited a significantly positive relationship with PAS of vegetation period (Table 3). Fig. 5. Pearson's correlation coefficients between basal area increments and climate variables for Quercus mongolica. Dashed horizontal lines correspond to significance levels at 0.05. VP, vegetation period (March-September); PY, physiological year (previous October-current September); T max , max temperature; T min , min temperature; T ave , mean temperature; P, precipitation; PAS, precipitation as snow (mm) between August in the previous year and July in the current year; CMD, climatic moisture deficit).
The fixed effects of the LMM related to climate variables explained 18% in CP, 40% in PGA, and 15% in PGB of the total variance of tree growth for Q. variabilis (marginal R 2 ), and 8 % in CP, 33% in PG, and 33% in HR for Q. mongolica (Table 3). When both the fixed and random effects (sites) were included in the models, the total explained variance (conditional R 2 ) increased considerably, by 29%, 44%, and 23% for Q. variabilis and 47%, 44%, and 46% for Q. mongolica, respectively (Table 3).

Climatic factors controlling tree radial growth
Understanding the growth response to climate among tree species with various range limits is fundamental to the assessment of forest dynamics and responses to forecasted climate warming. Previous studies indicated that the growth of tree species at their northern limit was restricted by low temperature (Babst 2013, S anchez-Salguero et al. 2015, but was limited by heat at their southern limit, for example, Pinus sylvestris (S anchez-Salguero et al. 2015), Picea abies (M€ akinen et al. 2002), and Q. alba (Goldblum 2010). Other studies have detected a positive relationship between radial growth and higher moisture at both the northern and southern edges of the distribution of Pilgerodendron uviferum and the opposite relationship at the core of its distribution (Holz et al. 2018). However, our results indicated that the growth of both Q. variabilis at the northern limit and Q. mongolica at the southern limit was limited by low temperature and the autumn CMD; that is, BAIs of both tree species were positively correlated with the minimum temperature but negatively correlated with the autumn CMD. This was consistent with previous studies on boreal tree species, such as Picea mariana and Pinus banksiana in eastern Canada, which indicated a positive relationship between growth and temperatures of the current-year winter and spring or the entire growing season (Huang 2010). These studies suggest that restriction factors varied in different study areas but were not necessarily different in the same study areas for different distribution areas of tree species. Furthermore, our results indicated that the growth of Q. mongolica in all three sites was also limited Fig. 7. Scatter plots of the principal component analysis loaded on correlation coefficients. (a) Chronologies, increasing circle size with increasing altitude; (b) climatic regressors. Wi, winter; Sp, spring; Su, summer; Au, autumn; VP, vegetation period (March-September); Py, physiological year (previous Octobercurrent September); T max , maximum temperature; T min , minimum temperature; T ave , mean temperature; P, precipitation; PAS, precipitation as snow (mm) between August in the previous year and July in the current year; CMD, Hargreaves climatic moisture deficit. Results are presented on the two first components. by CMD (Table 3). This finding is consistent with numerous previous studies showing that the growth of tree species (e.g., Fagus sylvatica, Q. alba, Picea abies) at their southern limit was mainly affected by water shortage (M€ akinen et al. 2002, Jump et al. 2006, Goldblum 2010. The growth of northern populations of Q. variabilis in all three sites was also affected by autumn water shortages (positively correlated with autumn precipitation and negatively correlated with autumn CMD), especially at low altitudes (CP; Fig. 6), although the LMM analysis presented a negative relationship between growth and precipitation (SPL-CP and APL-PGA in Table 3). However, our results were contrary to previous studies because the growth of populations at their northern edge/high latitudes was limited by water shortages (Figs. 5, 6, Table 3). A possible explanation for this discrepancy may be that the study areas were located in a semiarid region with low altitude. Previous studies indicated that drought was a major constraint of forest development and tree growth in semiarid regions (Barber et al. 2000, Jiao 2019. Many studies also found that precipitation was strongly limiting for tree species growth at their lower limits (e.g., F. sylvatica [Jump 2007]; Pinus pinaster [Candel-P erez 2012]; Pinus sylvestris [Matias et al. 2017]). For example, Tardif et al. (2010) found that the radial growth of Q. alba and Q. rubra at the northern boundary of their distribution was restricted by low temperature and drought. Latreille et al. (2017) confirmed that according to the elevation-dependent response patterns, the radial growth of trees was limited by drought at lower elevations, which is Notes: The climate variables are included as fixed effects, whereas plots are random effects. Marginal R 2 (M) reflects the variance explained by the fixed effects alone; conditional R 2 (C) reflects the total variance explained by the fixed effects and random structure of models; d, random-effect variance. T ave , mean temperature; P, precipitation; PAS, precipitation as snow (mm) between August in the previous year and July in the current year; CMD, climatic moisture deficit; VP, vegetation period (March -September); PY, physiological year (previous October-current September). SE, standard error. consistent with our results in that lower growth at low altitudes was found not only in the southern populations (MGL-CP) but also in the northern populations (SPL-CP; Appendix S1: Table S1).

Forecasted radial growth for the 21st century
The populations of Q. mongolica in PG and HR at their southern edge are expected to increase in BAI up to the end of the present century, although not at the same rate. The current growth rate is highest in the PG for Q. mongolica, and it will continue rising until 2080 under the RCP2.6 scenario, and then be maintained until 2100, which agrees with the trend of increasing temperature (+1.8°C under the 2.6 scenario, +3.6~3.7°C under the 8.5 scenario) and precipitation (+76.9 AE 0.6~143.8 AE 1.6 mm under the 2.6 scenario, +160.7 AE 1.1~250.0 AE 3.3 under the 8.5 scenario) detected in these areas (Table 2; Appendix S1: Fig. S3), although the CMD also will increase (Table 2; Appendix S1: Fig. S3), but disagrees with the recent observations of growth decline at the southern edge of their distribution for Pinus nigra, Pinus sylvestris, and Pinus uncinata (Andreu et al. 2007, Lapointe-Garant et al. 2010, and Cathaya argyrophylla and Pinus sylvestris (S anchez-Salguero et al. 2017). We found evidences suggesting expected growth improvement in cold-limited Q. mongolica at its southern edge, which also agrees with previous findings that growth is expected to increase because of longer and warmer growing seasons under the RCP8.5 emission scenario in mountainous areas of NE Spain (S anchez-Salguero et al. 2017). This research suggests that the The populations of Q. variabilis in PGA and PGB at their northern edge were also projected to increase in BAI up to the end of the present century. The current growth is similar between the PGA and PGB, and it will continue to increase until 2040 under the RCP2.6 scenario and be maintained until 2100, which agrees with the trend of increasing temperature (+1.8°C under the 2.6 scenario, +2.99-4.63°C under the 8.5 scenario) and precipitation (+134.4 AE 2.4~139.6 AE 1.8 under the 2.6 scenario, +145.5 AE 61.4~231.4 AE 62.2 mm under the 8.5 scenario) detected in these areas (Table 2; Appendix S1: Fig. S4). This implies that northern populations have the potential to continue to increase their growth after 2100, except at lower altitudes. However, the populations are expected to present a higher response under the predicted climate alteration of the RCP8.5 scenario. The BAI of both species (Q. mongolica in PG and HR, and Q. variabilis in PGA and PGB) increases up to the year 2100 and beyond (Table 2), which is consistent with Matias et al. (2017) who forecast a general increase in Scots pine growth by predictive models across its latitudinal distribution (including northern, central, and southern areas). One possible explanation is that climate warming could accelerate the decomposition rate of soil organic matter, thus enhancing the utilization efficiency of mineralized nutrients to plants, especially the most common limiting nutrient, N (Luo 2004). Furthermore, increasing the concentration of CO 2 in the atmosphere might have a positive effect on the above-/belowground biomass by improving plant photosynthesis and water-use efficiency in a CO 2 -enriched world, that is, the direct CO 2 fertilization effect (Huang 2007(Huang , 2010. In addition, the highest responsiveness of growth of Q. mongolica and Q. variabilis to climate both appeared at central altitude populations (MGL-PG, SPL-PGA/PGB). However, no significant increasing trends were found for the populations of Q. mongolica and Q. variabilis under the future climate change in CP where growth is most sensitive to moisture deficits. These projections imply that some southern Q. mongolica populations and northern Q. variabilis populations in relatively dry and hot areas (MGL-CP, SPL-CP) may compensate for growth declines forecasted during the late 21st century under the RCP2.6 and RCP8.5 scenarios. One possible explanation for our results is that the site CP was located at a low altitude, where warming and water shortages would be more intense (Table 3). High temperatures in autumn may accelerate the evaporation of soil water, causing drought stress and subsequently limiting tree growth, although precipitation would also increase (Fang et al. 2010). Additionally, high temperatures during the growing season could exacerbate respiration and limit the uptake of nutrients that are needed for growth (Liang et al. 2011).
Although linear mixed-effects model was often used to examine climate-growth relation and future prediction in the literature and our study (Gonz alez-Muñoz et al. 2014, Matias et al. 2017, there are possible problems when it is extrapolated since nonlinear relationships are far more likely (S anchez-Salguero et al. 2017). On the other hand, nonlinear relationship could be complex with interactions and difficult to be implemented. The process-based model can be used to simulate radial growth based on the principle of limiting factors and nonlinear growth response functions in the future (Fritts 2001, S anchez-Salguero et al. 2017).

CONCLUSION
Our study is unique in that it simultaneously described the radial growth of Q. variabilis at its northern range limit and Q. mongolica at its southern range limit. No significant differences were found for growth responses to climate between Q. variabilis at its northern limit and Q. mongolica at its southern limit. The radial growth of the two species was positively correlated with temperature and precipitation (autumn) and negatively correlated with CMD (autumn), although the growth response to climate variables was varied among different sites for the two tree species. Our models forecast a clear increase in radial growth of Q. variabilis at its northern limit and Q. mongolica at its southern limit up to the end of the present century, although water availability also affected the growth patterns (except at low altitudes), which would likely lead to a northward (Q. variabilis) and southward (Q. mongolica) expansion of the species' distribution. However, these predictions should be interpreted with care, because other factors, such as soil conditions, extreme climatic events, or biotic factors (e.g., stand structure), could strongly affect demographic performance. In summary, our study advances our understanding of the growth of woody species response to climate at the different range limits, and the results suggest a potential range shift of the oak species studied, which may occur across horizontal geographical gradients under climate change in boreal forests, although further evidentiary support is required from future studies.