The study area is located in northeastern China. It can be divided into three regions along the latitudinal gradient: south (44.10°–47.22° N, 120.28°–128.48° E), central (51.30°–51.62° N, 120.82°–123.53° E), and north (51.83°–52.59° N, 121.10°–123.52° E) (Table 1, Figure 1). From 2015 to 2018, we conducted field work at 12 sites, which we divided into three regions according to latitude: the southern sites (LBS, TS, AES), the central sites (DBS, ME, HM, YK, XBS), and the northern sites (AK, FKS, YA, ZL) (Table 1, Figure 1). The mean annual total precipitation (1950–2014) in the study area ranges from 415.6 to 636.2 mm (Figure 1A). The mean annual temperature is between −6.3 and 0.5°C. January and July are the coldest (−38.2°C, ZL) and hottest (25.3°C, TS) months, respectively (Table 2).

In the study area, P. pumila and other plant species form two main plant communities: the subalpine Empetrum nigrum-P. pumila community at elevations greater than 1,350 m a.s.l and the Alnus mandshurica-P. pumila community at elevations between 1,250 and 1,450 m a.s.l. The former is mainly distributed along flat mountain tops above 1,350 m, on open ridges, gravel heaps and the tops of mountains. The trunk of P. pumila creeps along the ground, with the crown rising obliquely to a height of about 1 meter. Growth is relatively weak and sparse, with 700–900 trees (clusters) per hectare. The main companion species include Juniperus sibirica, Betula rotundifolia, Rhododendron dauricum, Larix gmelinii, Empetrum nigrum var. japonicum, Vaccinium vitis-idaea, Cladonia alpestris, Sphagnum cymbifolium, and Hylocomium splendens. The latter grows under the canopy of boreal trees (Okuda et al., 2008). Alnus mandshurica grows vigorously and densely among P. pumila. The ground cover is composed mainly of Sphagnum girgensohnii, Vaccinium vitis-idaea, Ledum palustre var. angustum, Pyrola incarnata, etc. At high elevations, P. pumila usually grows between stones, while at lower altitudes, P. pumila grows in brown coniferous forest soil that contains a lot of sand and gravel. The sand content is 30–85%, and the gravel content greater than 2 mm is 3–35%. Although the organic matter content of brown coniferous forest soil is high, the content of nitrogen and phosphorus is low, and the ratio of carbon to nitrogen is high.

Because there are no weather stations nearby, we used the monthly and seasonal temperature and precipitation data of the CRU TS 4.04 0.5° × 0.5° grid to analyze the growth-climate relationships from 1950 to 2014. The data were extracted from the sampling area using the KNMI Climate Explorer web page.¹ The CRU database is formed by interpolating data from regional meteorological stations. Climate variables for growth-climate response analysis include monthly total precipitation (P) and mean (T_mean), minimum (T_min), and maximum temperature (T_max) (Figure 2). We define winter as from December of the previous year to February of the current year, spring as from March to May, summer as from June to August, and autumn as from September to November.

At each site, 10–15 discs of P. pumila were cut with a handsaw. These samples were obtained only from isolated, mature, and healthy individuals to avoid competition with other plants that may affect stem growth. A total of 140 discs from 12 sites (Figure 1, Table 1) were sampled along the latitudinal gradient from south to north. The cross-section of each disc was naturally air-dried and polished with increasingly fine grades of sandpaper (120–800 grit) until the ring boundary could be clearly distinguished under a microscope (Fritts, 1976; Cook and Kairiukstis, 1990). The rings of P. pumila were cross-dated and ring widths were measured in the laboratory using traditional dendrochronological methods (Fritts, 1976; Cook and Kairiukstis, 1990). Two to four radii in each disc were cross-dated, and the ring width at each radius direction was measured using the Velmex measurement system with a resolution of 0.001 mm. The cross-dating and measurement accuracy were statistically checked using the COFECHA computer program (Holmes, 1983), which determined the degree of synchronization between series according to the correlation with the main chronology.

Each ring-width series was detrended and standardized by fitting a negative exponential curve or linear line using the ARSTAN program to remove non-climate signals related to age or stand dynamic effects (Cook and Holmes, 1986). The ring-width index was obtained by dividing the ring width by the fitting value of each ring. Three kinds of ring-width chronologies (standard, residual, and autoregression chronologies) were obtained by averaging all detrended series with a bi-weight robust mean (Cook and Kairiukstis, 1990). Raw measured ring widths and the standard chronologies (STD) were used in the subsequent analyses.

Chronology parameters are statistical using the ARSTAN program. The expressed population signal (EPS), defined as the proportion of each series signal of the total series variance, was used to quantify the reliability of the chronology (Wigley et al., 1984). The mean sensitivity (MS) and the first order autocorrelation (AC1) were calculated on the individual detrended index series, then averaged to measure the year-to-year variability and the degree to which current-year growth was influenced by previous-year climate factors. The mean series correlation between trees (R_bar) allows us to evaluate the strength of the common growth signal over time. The standard deviation (SD) of inter-annual ring-width variability was calculated as a proportion of mean ring-width (Chen et al., 2011). Variance in the first eigenvector (VF1) of all series identifies the common growth variability among all trees at each site. The signal-to noise ratio (SNR) is a measure of the strength of the common high-frequency signal in the ring-width indices of trees from the same site.

Pearson correlation was used to analyze the relationship between the chronology and monthly and seasonal climate variables to compare the climate factors limiting radial growth in different regions. The climatic factors from May to October of the previous year were selected for analysis because the climate of the previous and current year affects the radial growth of trees in the current year (Fritts, 1976). We also carried out correlation analyses for two periods (1950–1980 and 1981–2014) to investigate the impact of recent rapid warming (after 1980) on the growth of P. pumila.

The responses of P. pumila individuals to temperature were categorized as one of four response patterns: positive (more than 67% of individuals at a site are significantly positively correlated with temperature), negative (less than 33% are significantly positively correlated with temperature), mixed (between 33 and 67% are significantly positively correlated with temperature), or none (no significant correlation with temperature) (Lloyd et al., 2011). The response patterns of P. pumila to precipitation were similarly classified, and the proportion of shrubs exhibiting each response pattern was tallied for each site. Then, linear mixed-effects models (LMEs) were used to identify the effects of eight climate variables (T_wi: temperature in winter, T_sp: temperature in spring, T_su: temperature in summer, T_au: temperature in autumn, P_wi: precipitation in winter, P_sp: precipitation in spring, P_su: precipitation in summer, and P_au: precipitation in autumn) on the residuals of the previous function. Climate variables during 1950–2014 were used as fixed factors, while individual shrubs were considered random factors. Fitted models follow the equation:

where RW_i represents the ring width at year i; a is the vector of fixed factors (seasonal climate variables), b is the vector of random factors (shrub identity), X and Z are the fixed and random effects regressor matrices, respectively, and ei is the within-group error vector (Camarero et al., 2017). We ranked all the potential models that could be generated with the different explanatory variables according to the Akaike information criterion (AIC). We chose the most parsimonious models; that is, the models with the lowest AICs (Burnham and Anderson, 2002). For each site, the model with the lowest number of variables among those with the lowest AIC was selected as the final model (Burnham and Anderson, 2002). In addition, we used the Akaike weights (Wi) of each model to measure the conditional probability of the candidate model, assuming it was the best model. By using seasonal climate averages instead of monthly data, we can create more parsimonious models while maintaining a reliable representation of climate trends (Matías et al., 2017). The linear mixed model analysis was performed using the LME4 package in R 4.0.3 (R Core Team, 2015). Hierarchical cluster analysis and correlation analysis were performed using the SPSS 18.0 software package (IBM SPSS Inc., Chicago, IL, USA). Finally, the ring-width growth trend of P. pumila under four emission scenarios (RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5) was simulated through the linear regression prediction under the CMIP5 scenarios of 2020–2100 predicted by HadGEM2-ES.

The mean age of shrubs at the 12 sites ranges from 69 to 151 years, and the mean ring width varies from 0.27 to 0.65 mm (Table 2). ME has the highest MS, AK has the highest SD, VF1 and R_bar, YK has the highest AC1, and HM has the highest EPS and SNR (Table 2). In short, the chronological statistics of P. pumila at different sites show significant trends. All analyses indicate that the 12 chronologies have abundant climate signals and are suitable for analyzing growth-climate relationships. The Pearson correlations analysis result shows a high correlation among sites in the same region (Figure 1C). From 1950 to 1980, the ring-width index of P. pumila increased significantly in the south and central regions, but did not increase significantly in the northern region. After 1980, the chronology trends diverge in different regions, with the southern and central regions showing a slight downward trend, and the northern region showing an upward trend (Figure 2A). There is no significant difference in ring width between different regions (Figure 2B).

The chronology of P. pumila in the north is significantly positively correlated with temperature in January and May (Figures 3A–C). Warm winter conditions enhance the growth of P. pumila in the north and south (Figures 3A–C). In addition, wet winter conditions are positively correlated with the growth of P. pumila in the north. The ring widths of P. pumila are significantly correlated with minimum temperature from April to June, especially in the north. Wet conditions in June promote the growth of P. pumila in all three regions. It should be noted that there is no significant correlation between the growth of P. pumila and monthly maximum temperature in summer. However, warm autumn conditions are conducive to the growth of P. pumila in the north and south. Temperature during the previous growing season has a positive effect on the radial growth of P. pumila. The results show that the growth response of P. pumila varies with latitude, and that each region has at least three different response classes (Figure 3).

Northern P. pumila individuals respond more negatively to temperature than their southern and central counterparts, but many individuals exhibit a positive relationship with temperature (Figure 4). Moreover, P. pumila individuals that are positively correlated with temperature show more growth in colder years (Figure 4A).

During the period 1950–1980, the growth of P. pumila was negatively correlated with summer temperature, positively correlated with winter temperature, and positively correlated with summer and winter precipitation (Figure 5). However, from 1981 to 2014, the growth of P. pumila was positively correlated with summer temperature and negatively correlated with winter temperature with increasing latitude (Figure 5). This may be due to the higher winter temperature growth rates at high latitudes from 1981–2014. In the more recent 30-year period, the correlation between the growth of P. pumila and summer precipitation decreased with increasing latitude, while the correlation between the growth of P. pumila and winter precipitation increased (Figures 4, 5).

We developed a prediction model for the relationship between the growth of P. pumila and the main climate factors. For the climate change predicted by the HadGEM2-ES RCP 2.6 scenario, our model predicts that the growth of P. pumila in the central region will stabilize before 2100 (Table 3, Figure 6). In contrast, the growth of P. pumila increases slightly in the south and decreases slightly in the north. Under the RCP 4.5 and RCP 6.0 scenarios, the growth of P. pumila in the south increases rapidly. Under the RCP 8.5 scenario, the growth of P. pumila in the south increases rapidly until 2100 (Figure 6). In contrast, the growth of P. pumila in the northern and central regions shows a downward trend in the RCP 4.5, RCP 6.0, and RCP 8.5 scenarios (Figure 6).

The growth of P. pumila in the south is limited by temperature during the growing season. This might be because the average elevation of the southern sampling sites (1,310 m) is higher than that of the central and northern sites; average annual precipitation is also greater (Li et al., 2020). Therefore, moisture is not the main factor limiting the growth of P. pumila in the southern part of the study area. The increase in winter temperature of the previous year promotes the growth of P. pumila in the southern and northern regions but has little impact on the growth of P. pumila in the central region. High temperatures in winter accelerate snow melting, providing sufficient water supply during the subsequent growing season and enhancing radial growth rates (Duan et al., 2017; Lyu et al., 2017). In addition, warm winters help protect needles and fine roots from frost damage (Myers-Smith et al., 2011; Duan et al., 2017). The timing of soil thaw plays a crucial role in initiating the physiological activity of shrubs (Vaganov et al., 1999; Kirdyanov et al., 2003; Lloyd et al., 2011). Therefore, winter warming may be beneficial to the growth of high-altitude shrubs.

The sensitivity of P. pumila growth to temperature varies along the latitudinal gradient. The positive correlations between P. pumila growth and minimum temperature highlight the importance of winter and spring temperatures in the north (Figure 4). Frost may affect shrub growth at high latitudes before the growing season (Li et al., 2020). In a cold environment, spring temperatures and the length of the growing season are important factors limiting the production and differentiation of xylem cells (Morales et al., 2012; Lenz et al., 2013; Yang et al., 2020). Higher temperatures in early spring may favor coniferous species because their evergreen leaves enable them to take advantage of early warmth (Jiang et al., 2016). Myers-Smith et al. (2011) found a significant positive correlation between the ring widths of some of the most common shrub species in the tundra ecosystem and the mean growing season temperature. However, winter and spring temperatures seem less critical to the growth of P. pumila in the central region as compared to the northern and southern regions. The impact of maximum summer temperature on the growth of P. pumila is not significant, suggesting that the initial climate conditions of the growing season are more important for a positive growth response (Figure 4). However, we found that high temperature and insufficient water in summer still has a negative impact on the growth of P. pumila (Yang et al., 2020).

At the individual level, the responses of P. pumila growth to precipitation are consistent, with precipitation more strongly influencing growth in the northern areas. The more northward P. pumila grows, the more negative its growth response to temperature is. This is inconsistent with the results of previous studies, which show a pattern of more positive responses to warming in northern region (Lloyd et al., 2011). However, we found that more than half of the northern P. pumila individuals have a negative response to temperature, which may be caused by the synchronous constraining effects of hydrothermal coupling on the growth of P. pumila (Yang et al., 2020). Individual shrubs differ in many respects, including size, age, microsite conditions, competition, and genome (Huebner and Bret-Harte, 2019). The population-based approach emphasizes the mean climate signal, but the cost is the loss of individual-level information (Galván et al., 2014). Thus, although population chronology can highlight regional climate information, it may bias subsequent inferences. This is especially true for dwarf shrubs growing in harsh environments (Sedlacek et al., 2016), as each shrub individual may experience different microenvironment conditions. The difference between individuals in the population is the reality of nature, whereas the mean value is a statistical abstraction (Carrer, 2011). However, Bär et al. (2008) found that while the topoclimate at different microsites leads to slight changes in the annual growth pattern of Emperum hermaloditum, the prominent regional climate signals are still reflected in shrub growth at all microsites.

The altitude and longitude that cause environmental heterogeneity may also affect the growth and physiology of shrubs (McDowell and Sevanto, 2010). High-altitude shrubs are more limited by low temperatures, whereas low-altitude shrubs are more limited by water (Lu et al., 2021). The precipitation difference caused by longitude also affects the relationship between shrub growth and climate factors and may further affect their responses to drought and warming (Li et al., 2016; Anderegg et al., 2019; Buchwal et al., 2020). Therefore, our results may include the combined effects of latitude, longitude, and elevation on shrub growth.

From 1950 to 1980, the radial growth of P. pumila is limited by moisture deficits in summer and low winter temperatures. However, the growth of P. pumila in the northern region is more susceptible to winter precipitation from 1981 to 2014 due to climate warming. The weakening of temperature seasonality in the Northern Hemisphere due to anthropogenic forcing, as well as future human impacts, may further weaken the annual temperature cycle, thereby affecting biological and ecological systems, especially at high latitudes (50°–70 °N) (Rosenzweig et al., 2008; Qian and Zhang, 2015; Duan et al., 2017). The degree of warming in northern China over the past 60 years has exceeded that of southern China (Xu et al., 2021). Seasonal warming during the 1958–2017 period was highest in winter and lowest in summer (Xu et al., 2021). Winter warming has exceeded summer warming since the end of 19th century, and the magnitude of the annual temperature cycle has decreased by 0.45°C/100 a (Duan et al., 2017). The rapid warming since 1980, particularly in winter, will affect the radial growth of P. pumila, with a more pronounced influence at higher latitudes (increased drought). If warming continues, water stress may become the main factor limiting the growth of P. pumila in the future.

The growth of P. pumila in the north shows an unexpected but clear downtrend under the climate warming scenario, whereas growth in the south shows an upward trend. The relationship between climate factors and P. pumila growth is not perfectly linear due to the influence of multiple climate factors (Shao et al., 2021; Zheng et al., 2021). In addition, climate also compensates the radial growth of trees or shrubs (Jiang et al., 2015; Cao et al., 2016). The model results may be affected by the linear relationship between climate factors and shrub growth (Bonan, 2008). Many studies have confirmed that shrub growth increases and expands to higher latitudes and altitudes (Myers-Smith et al., 2015b; Myers-Smith et al., 2020) as global warming leads to extended growing seasons (Huang et al., 2011; Rossi et al., 2016). Compared with shrubs, the proportion of trees (larch and pine) threatened by high summer temperatures will increase by 45.5 to 63.6% (Li et al., 2020). In our study area, the degradation of permafrost due to the increase in ground surface temperature has led to a rapid increase in the growth of larch in the past decade (Zhang et al., 2019). However, the impact of soil warming on the growth of P. pumila, which has a shallow root system, remains unclear. Under different scenarios of HadGEM2-ES, the growth of P. pumila in the south exhibits an increasing trend until 2100 (Figure 6). In contrast, the growth of P. pumila in the northern and central regions decreases slightly or shows a steady trend (Figure 6). Under the future warming scenario, the growth of P. pumila in the northern and central regions will be significantly reduced. This is mainly because water plays an important role in shrub growth in these regions. When there is sufficient moisture, shrub growth benefits from the continuous warming (Myers-Smith et al., 2015a). Therefore, in areas where trees are sensitive to warming and drought due to rising temperatures, shrubs may become the dominant plant life form. In these areas, even if shrubs gradually replace trees, it will take a long time. Conversely, changes in shrub distribution and coverage may also have a feedback impact on climate change (Myers-Smith et al., 2011).