Study area
The study area is located in northeast China, with latitudinal gradients from 44.06° N to 52.36° N and 120.17° E to 128.29° E (Table 1, Fig. 1). From 2015 to 2018, we conducted a field investigation at 12 sites. We divided the 12 sampling sites into three regions according to latitude: South (LBS, TS, AES), Central (DBS, ME, HM, YK, XBS), and North (AK, FKS, YA, ZL) (Table 1, Fig. 1). The mean annual precipitation (1950–2014) in the study area ranges from 415.6 mm to 636.2 mm, and more than 68% occurs from June to August (Fig. 1a). The mean annual temperature is between −6.3 ℃ to 0.5 ℃. January and July are the coldest (−38.2 ℃, ZL) and hottest (25.3 ℃, TS) months, respectively (Table 2). The annual frost-free period is 80-100 days, with early and late frost occurring in September and May.
P. pumila is a multi-trunk, ground-creeping shrub which grows obliquely. In the symbiotic community with other plant species, P. pumila forms two types of plant communities: high-altitude subalpine shrub community and low-altitude high canopy tree forest community. The former grows in subalpine areas with fruticulose and herbs. The latter grows under the canopy of boreal trees (e.g., Larix gmelinii, Pinus sylvestris var. mongolica, Betula ermanii) (Okuda et al. 2008). The dominant P. pumila in the subalpine dwarf forest is mainly distributed in the altitude range of 800-1700 m in northeast China. Pinus sylvestris var. mongolica and Larix gmelinii are the main trees associated with P. pumila at low altitude, and they are also the main constructive species of boreal forest in the Daxing’an Mountains, China (Zhang et al. 2019).
Table 1 Characteristics of the 12 study sites with the name of the mountain, the sample code and sample number, longitude and latitude of site, altitude of plots above sea level, and divided region.
Site
|
Code
|
Longitude
(E)
|
Latitude (N)
|
Altitude (m a.s.l.)
|
Sample number
|
Region
|
Laobai Mountain
|
LBS
|
128° 03′
|
44° 06′
|
1685
|
35
|
South
|
Tao Mountain
|
TS
|
128° 29′
|
46° 38′
|
1369
|
31
|
South
|
A’er Mountain
|
AES
|
120° 17′
|
47° 13′
|
1164
|
52
|
South
|
Dabai Mountain
|
DBS
|
123° 08′
|
51° 18′
|
1431
|
28
|
Central
|
Mo’erdaoga
|
ME
|
120° 49′
|
51° 22′
|
1072
|
33
|
Central
|
Hanma
|
HM
|
122° 24′
|
51° 31′
|
1000
|
58
|
Central
|
Yikesama
|
YK
|
121° 14′
|
51° 32′
|
1262
|
36
|
Central
|
Xiaobai Mountain
|
XBS
|
123° 32′
|
51° 37′
|
1400
|
33
|
Central
|
Aokelidui Mountain
|
AK
|
122° 03′
|
51° 50′
|
1104
|
30
|
North
|
Fuke Mountain
|
FKS
|
121° 40′
|
52° 28′
|
1096
|
30
|
North
|
Yong’an
|
YA
|
121° 06′
|
52° 33′
|
1196
|
60
|
North
|
Zhalinku’er
|
ZL
|
123° 31′
|
52° 35′
|
1135
|
58
|
North
|
Dendrochronological sampling, treatment and analysis
At each site, 10-15 discs of P. pumila were sliced with a hand saw near the base of the trunk. These samples were obtained only from isolated, mature and healthy individuals to avoid the fact that competition with other plants may affect stem growth. A total of 140 discs from 12 sites (Fig. 1, Table 1) were sampled on the latitude gradient from south to north. The cross-section of each disc was naturally air-dried and polished with mesh gradually finer sandpaper (120-800 grit) until the ring boundary could be clearly distinguished under the microscope (Fritts, 1976; Cook and Kairiukstis, 1990). Shrub rings were cross-dated and ring width was measured in the laboratory using traditional dendrochronological methods (Fritts, 1976; Cook and Kairiukstis, 1990). Two to four radii in the 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 determines the degree of synchronization between series according to the correlation with the main series of 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 derived 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). The standard chronologies (STD) were used in the subsequent analyses.
Climate data
We used CRU TS 4.04 0.5° × 0.5° gridded monthly and seasonal temperature and precipitation data to analyze growth-climate relationships for 1950–2014 because no nearby weather stations exist. The data were extracted from the sample area using the KNMI Climate Explorer web page (http://climexp.knmi.nl). The CRU database is formed by interpolated values from regional meteorological stations. Climate variables including monthly total precipitation (P), mean (Tmean), minimum (Tmin) and maximum temperature (Tmax) (Fig. 2) were used for growth-climate response analysis. We define winter as December of last year to February of that year, spring for March-May, spring for June-August, and autumn for September-November.
Climate–Growth Relationship Analysis
To determine the main climate factors limiting the radial growth of each chronology at each site. Pearson correlation was used to determine the relationship between the tree-ring index and monthly and seasonal climate variables. Radial growth is affected by the current and previous year’s climate (Fritts, 1976). Therefore, climate variables over 17 months, from May of the previous year to October of the current year, were used for the correlation analysis. Meanwhile, we carried out two periods (1950-1980, 1981-2014) correlation analyses to investigate the temporal stability of the growth-climate response.
The response of each shrub to temperature was categorized into four response types: positive (67% of significant correlations with temperature were positive), negative (33% of significant correlations with temperature were positive), mixed (between 33% and 67% of significant correlations with temperature were positive), or none (no significant correlations with temperature) (Lloyd et al. 2011). Shrubs were similarly categorized concerning their response pattern to precipitation, and the proportion of shrubs exhibiting each response type was tailed for each site. After that, linear mixed-effects models were used to identify the effects of 8 climatic variables (Twi, Tsp, Tsu, Tau, Pwi, Psp, Psu and Pau) on the residuals of the previous function, using climatic variables during 1950–2014 as fixed factors and shrub as a random factor. Fitted models followed the equation:
where RWi represents the ring-width in year i; a is the vectors of fixed effects (seasonal climate variables), b is the vector of random effects (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 selected the most parsimonious models, that is the ones with the lowest AIC (Burnham & Anderson, 2002). The final model was selected for each site as the one with the lowest number of variables among those with the lowest AIC (Burnham & 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. The use of seasonal climatic averages instead of monthly data allowed the creation of more parsimonious models while maintaining a reliable representation of climatic 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 22.0 software package (IBM SPSS Inc., Chicago, IL, USA).
Finally, the growth trend of shrub ring width under four emission scenarios of RCP 2.6, RCP 4.5, RCP 6.0 and RCP 8.5 was simulated by a linear regression forecast under the CMIP5 scenario predicted by HadGEM2-ES for the period 2020–2100.