Ecological contingency in species shifts: downslope shifts of woody species under warming climate and land-use change

A predicted impact of a warming climate is an upslope shift of montane plant species. These upslope shifts may be amplified by land-use changes or attenuated by forest recoveries at low elevations where historical disturbances were ceased allowing for plant regrowth. Consequently, species may shift downslope back to low elevations where they had been previously harvested. The cessation-driven downslope shifts are hypothesized to dampen or even reverse climate-driven upslope shifts. We tested this hypothesis by a 20 year (1989–2009) forest inventory dataset from five mountainous areas in eastern China. In our study region, intense deforestation occurred mostly at low elevations until 1970, but was then ceased to facilitate natural forest recovery. Based on the analyses of 30 216 woody plants in 609 plots, we found that: (1) forest recovery occurred over the 20 year survey period, and increment rates of both recruitment and basal area increased up to 2004. However, in the last period (2004–2009), increment rates of basal area leveled off and recruitment was close to zero; (2) forest recovery was faster at lower elevations, as indicated by the higher increment rates there; (3) despite rising regional temperatures, the mean elevations of study species showed a downslope shift over the 20 years; and (4) the contribution of forest recovery to elevational shifts was supported by the fact that the species shifts were positively related to elevational changes in the recruitment increment, e.g. the negative (or downslope) shifts occurred in association with higher increments at lower elevations. These results suggest that, the cessation of disturbances and consequent lowland forest recovery had greater effects on the species distributions than did warming climate. In mountain systems that are being allowed to recover from historical disturbances, the effects of forest recovery on species distributions should be explicitly accounted for when assessing and predicting climate change impacts.


Introduction
Plant species in mountain systems are generally predicted to shift their distributions upslope in response to warming climate (Breshears et al 2008, Feeley et al 2012, Kuhn and Gégout 2019, but downscaling the impacts to regional scale with various local conditions could lead to ecological contingencies (Harrison et al 2010), e.g. downslope species shifts caused by changes in highland permafrost (Cannone et al 2007) or climatic water balance (Crimmins et al 2011). In addition, downslope shifts could be potentially caused by land-use changes (Lenoir and Svenning 2015), such as cessation of human disturbances (e.g. logging and grazing) which contributes growth space to forest recovery (Bodin et al 2013, Chisholm et al 2014). Location of the forest recovery along elevational gradients potentially determine the direction of species shift, and lowland forest recovery may lead to downslope shifts since species can expand their realized niches into lower elevations (Lenoir et al 2010a). However, the net effects of warming climate and disturbance cessation on elevational species shifts have been poorly investigated according to a recent review (39 studies in Guo et al 2018). Understanding the net effects of both processes will improve predictions of climate change impacts on species distributions at regional scale, which is more relevant to local management of natural resources. Here, we tested the hypothesis that the disturbance cessation may obscure climate-driven upslope shifts, and even drive downslope shifts of plant species.
The cessation of human disturbances should be viewed as another form of global change that can have marked impacts on species distributions (Jackson and Baker 2010). The disturbance cessation is associated with two causes. The first one is local economic transformations, e.g. abandonments of agricultural lands in Europe, particularly in highland mountain areas, have been motivated by their reducing productivity and increasing focus of agricultural activities on more fertile Notably, these highland abandonments in Europe are contrasting with our study sites where the disturbance cessation occurred at low elevations (Gehrig-Fasel et al 2007). The other cause is increasing concerns about the negative influences which human disturbances impose on ecosystems, such as the degradation of ecosystem services resulting from deforestation (Chazdon and Brancalion 2019). Consequently, there is an increasing number and size of protected areas or protection programs. For example, several programs have been launched to protect and increase forest coverage in China (Nüchel and Svenning 2017), and over half of the world's tropical forests are naturally recovering from different forms of deforestation (Food and Agriculture Organization of the United Nations 2010).
Investigating the net effects of these two aforementioned global drivers (warming versus disturbance cessation) requires a long-term dataset of forest surveys across elevations in an area where disturbances have been ceased. In this study, we used a 20 year (1989-2009) forest survey dataset sampled from four mountainous areas of eastern China. In their low elevations, intense deforestation occurred until 1970 and were then ceased to facilitate natural forest recovery. Concurrently, temperatures in our study sites have increased steadily while precipitation has remained relatively stable. Capitalizing on the data from these areas with the unique conditions, we examined the process of forest recovery by analyzing the temporal patterns of forest growth (e.g. increment rates of recruitment and basal area). During forest recovery, forest growth rates generally increase over time. This increasing pattern contrasts with the relative stable and low growth rates of mature forests (Bormann 1979).
Additionally, elevational patterns of forest growth were analyzed by the changes in growth rates with elevation. We tested these hypotheses: (1) the increment rates of both recruitment and basal area are increasing through time, and the increment rates are higher at low versus high elevations; (2) plants would exhibit primarily downslope shifts; (3) the elevational species shifts are associated with elevational changes in recruitment rates.

Study area and human disturbance history
The study was conducted in the mountain regions of Anhui province in eastern China (figure 1). Based on a 30 year climate dataset (1983-2013) from three meteorological stations covering our study region (Hefei, Tongling and Huangshan in figure 1), the local climate is characterized by warm humid summers and cool dry winters, with most precipitation occurring from April to August. Figure 2 shows both the mean annual temperature and annual sum of precipitation over the past 30 years.
Prior to 1970, the dominant form of local disturbance in this study region was deforestation for timber and fuel woods. The disturbance levels were estimated by both the number of individuals cut and conditions of disturbed soil in 1989. The lowland plots suffered more from disturbances than ones located at high elevations (figure A1(A)), since lowland areas can be accessed easier than highland ones. In 1970, a project named 'Closing Hillsides to Facilitate Afforestation' was initiated under the Chinese conservation policy of 'Conversion of Degraded Farm Land into Forest and Grass Land' and related conservation law of forest resources. Therefore, both artificial deforestation and natural disturbances (e.g. fire) are strictly limited by local government agencies (e.g. Forest Service), which monitored these forests closely. The disturbance level was only estimated in 1989. There have been no subsequent largescale disturbances and the small-scale disturbances may occur, but their effects (e.g. mortality) are not detected from our forest survey data.
The natural forest recovery proceeded from surviving small individuals (or seedlings) and soil seed bank. Both the low-and highlands showed an increasing number of individuals per plot from 1989 to 2009 (figure A1(B)). Given this strict cessation of disturbances, our study sites provide an ideal opportunity to investigate the effects of disturbance cessation on elevational species shifts.

Forest dynamics: increment rates of recruitment and basal area
At a given elevation (e) from a given study site ( j), increment rates of recruitment ( + R e j y y , , , 1 ) and basal area ( + BA e j y y , , , 1 ) between the yth and yth +1 years were calculated as: where N p e j y , , , was the total number of living individual in a given plot (p) , and n e j , was the total number of plots at a given elevation (e) from a given study site ( j) during a given survey year (y).
2.4. Species mean elevation (SME) and species elevational shift (SES) SME was used to measure the elevational center of 'gravity' or optimum elevation of the study species. Following the protocols of Chen et al (2009) and Feeley et al (2011), we calculated SME as the mean elevation weighted by number of individuals per elevation in each census. We then calculated the 20 year change of SME from 1989 to 2009 as a measure of SES (equation (3)). ,2009 , ,1989 where SES i was the species range shift of a given species MSE i j , ,1989 were the mean species elevation of a given species (i) at a given study site ( j) in 2009 and 1989, respectively. J was total number of study sites (J=5 in this study). The change of SME was used to estimate the SES instead of elevational margins, since the margin-based estimates are sensitive to sampling efforts and short-term disturbances (e.g. hurricane (Zhai et al 2019)) that can obscure long-term species shifts (Lenoir et al 2008). For example, a study of mangroves showed that the margin-based northward expansions could be a result of short-term recovery from hurricanes, instead of a long-term response to warming climate (Giri and Long 2014).
2.5. Statistical analysis 2.5.1. Temperature and precipitation Linear mixed-effects models were used to analyze temporal (1983-2013) patterns of temperature (equation (4)) and precipitation (equation (5) where Temperature y j , and Precipitation y j , were the temperature and precipitation in a given year (y) at a given meteorological station ( j). b 1 was the coefficient of year, and Year y was the year of y. b 0 was the intercept, and a j was the random intercept. Both e of equation (4) and e j of equation (5) were the residuals, and residuals assumed homogeneity in equation (4) but heterogeneity in equation (5) between the stations. The model selection and assumption examination followed the procedures in Zuur et al (2009). The model selection was based on Akaike information criterion of models with different variance and covariance (or random) structures, but fixed-effect structure was the same since only the year effect was used in this case and thus restricted maximum likelihood was used to compare the models and estimate model parameters. Meanwhile, maximum likelihood was used to compare models with different fixed-effect structure (e.g. with and without the interaction of elevation and survey year in the 2.5.2). The violation of homogeneity and independence were examined by checking residual plots along the fitted values and years. The independence assumptions (i.e. temporal independence along the years) were further examined by checking correlogram plots of standardized residuals of the above models. The above statistical analyses were made by R program (Ball 1980) and the 'nlme' package (Pinheiro et al 2019).
2.5.2. Increment rates of recruitment and basal area Linear mixed-effects models were used to analyze temporal (i.e. three survey periods used for recruitment, and four used for basal area) and elevational patterns of increment rates of recruitment (equation (6)) and basal area (equation (7)). Note that the last period was dropped for recruitment analysis since >95% of the observations were zero.
where R e j s , , and BA e j s , , were the increment rates of recruitment and basal area at a given elevation (e) from a given study site ( j) during a given survey period (s), respectively. R e j s , , was log transformed and 0.1 was added to avoid zero. b 1 was the coefficient of elevation, and Elevation e j , was the elevation a given study site ( j). b s 2, was the coefficient of one survey period, and Survey s was one of the survey periods, including 1989-1994, 1994-1999, 1999-2005, and 2005-2009. b 0 was the intercept, and a e j , was the random intercept. In equation (6), e e j s , , was the residual, and assumed heterogeneity between study sites and survey periods, heterogeneity within each study site with elevation and the heterogeneity strength with elevation was different for each study site. Beside, e e j s , , included a temporal correlation along the three survey periods. In equation (7), e e j , was the residual, and assumed heterogeneity per study site, heterogeneity within study site with elevation and the heterogeneity strength with elevation was the same between each study site. Beside, e e j , included a temporal correlation along the four survey periods. The model selection and assumption examination followed the same procedure in 5.1 (see the correlograms in figures A3 and A4).
On the basis of above models (with two fixed factors of survey period and elevation, and one random factor of study site), partial residual plots (or componentplus-residual plots) of the two increment rates (i.e. recruitment and basal area) were used to visualize temporal and elevational patterns in the presence of all the other factors, respectively (Chatterjee and Hadi 1988).

Species mean elevation (SME)
Linear mixed-effects models were used to analyze the temporal (i.e. five survey years) patterns of SME (equation (8) where SME i j y , , was the mean elevation of a given species (i) at a given study site ( j) in a given year (y). b y 1, was the coefficient of one survey year, and Year y was one of the survey years, including 1989, 1994, 1999, 2005, 2009. b 0 was the intercept, and a i j , was the random intercept. e i j , was the residual and assumed heterogeneity between species and study sites. Beside, e i j , included a temporal correlation along the survey years. The model selection and assumption examination followed the same procedure in 2.5.1, and result visualization followed the same procedure in 2.5.2.

Species elevational shift (SES)
During the 20 year survey period, the elevational patterns of recruitment rates were analyzed for each study species in this study by linear mixed-effects models (equation (9), Note: the differences of this recruitment analysis from the previous one in 2.5.2 is that (1) this analysis was based on the species level instead of grouping species together at each elevation; (2) this analysis was based on the 20 year survey period instead of 5 year). From the equation (9), the elevation coefficients (b i 1, ) and their p-values were calculated for each species. The relationship between elevation coefficients (b i 1, of the equation (9)) and SES was quantified by a linear regression model (equation (10)) weighted by the variances of SES between the five study sites, i.e. a species with lower variance had larger weights in the regression model. The weighting focused this analysis more on the species with more consistent shifting patterns between the five sites. where R i and SES i were the increment rate of recruitment and SES of a given species (i), respectively. In equation (9), b i 1, was the coefficient of elevation, Elevation e j i , , was a given elevation (e) of a given species (i) at a given study site ( j), b i 0, was the intercept, a j i , was the random intercept, and e e j i , , was the residual (Note: the random factor and variance and covariance structures of model residuals may vary by species, and were determined by the same procedure in the~2.5.1). In equation (10) (coefficient of elevation in equation (9)), and b 0 was the intercept. e i was the residual and its variance would vary with var , i which is the SES variance of a given species (i) between the J sites (J=5). w i was the weights of linear regression. SES i was the mean SES of a given species (i) across the J sites.

Temporal and elevational patterns of increment rates of recruitment and basal area
From 1985 to 2013, temperatures of the study region showed a significant warming pattern (p<0.0001, figure 2(A)). In contrast, precipitation did not show any significant trends over time ( figure 2(B)). The increment rates of both recruitment and basal area showed significant increases during the survey periods (p<0.0001, figure 3). However, there were few new individuals observed in the last survey period (2004)(2005)(2006)(2007)(2008)(2009), and the increment rate of basal area also leveled off. Moreover, both of the two increment rates decreased significantly with elevation (p<0.0001, figure 4), and the slopes of elevational changes (b 1 in equations (6) and (7)) were −0.0013 and −0.0006 for recruitment and basal area, respectively. These decreasing elevational patterns suggests faster forest recovery at lower elevations.

SME and SES
Averaging the 29 study species from five study sites (equation (8)), SMEs decreased over time (p<0.0001, figure 5(A)), but leveled off during the last two survey years (2004 and 2009). Moreover, the SESs were significantly related to the elevation coefficients (b i 1, ) in equation (9) (p<0.0001, figure 5(B)), suggesting that the negative (or downslope) SESs might be driven by negative values of the elevation coefficients, i.e. higher recruitment rates at low elevations. Moreover, eight species, labeled in figure 5(B), showed significant elevational patterns of recruitment rates (p<0.05 based on equation (9)). In these eight species, six had negative coefficients of elevation (or high recruitment rates at low elevations), and five of these six species (Quercus variabilis, Quercus glauca, Quercus chenii, Celtis sinensis and Litsea coreana) showed downslope shifts.

Temporal and elevational patterns of forest recovery
Due to the implementation of new conservation efforts, deforestation in our study region ceased in 1970. The subsequent recovery of forests is evident in the temporal and elevational patterns of forest recruitment and growth. Specifically, the increment rate of basal area increased until the third survey period and then leveled off after 2004, suggesting that the forests had reached or were approaching stand maturation (Zhai et al 2015). Stand maturation was also indicated by the fact that very few new individual trees recruited into the plots during the last survey period. As stands mature, the closing forest canopies can lead to a light limitation in the understory and decrease the occurrence of new seedlings.
Forest recovery was faster at low elevations as indicated by the decreasing increment rates of recruitment and basal area with elevation (i.e. higher rates at low elevations). The faster recovery may result from lower competition after the cessation of disturbances (Suttle et al 2007) and reduced abiotic stress (e.g. more benign temperatures) at low elevations (Vittoz et al 2009, Küchler et al 2015. Additionally, the faster natural recovery of forests may be attributable to the seeds of low-elevation species generally recruiting faster than high-elevation ones (Kueppers et al 2017). Notably, with stand maturation and canopy closure, the linear decrease in increment rates with elevation (figure 4) may switch to other forms in the future as plants in lowland forests suffer from amounting stresses, such

SME and SES
On average, the SMEs decreased over time for the study species, and the SESs were highly related to the elevational patterns of recruitment such that the downslope shifts might be driven by higher recruitment rates at lower elevations. The fact that some species showed significant downslope shifts despite consistently rising temperatures is likely attributable to the forest recovery at low elevations following the cessation of deforestation in 1970. This supports our hypothesis that land-use changes can dampen or reverse climate-driven species shifts towards higher elevations. In addition to the downslope shifts, it is likewise possible that some upslope shifts are driven by forest recovery following agricultural land abandonments and cessations of forest harvesting at high elevations ( To understand and predict these contingencies, we will need to be able to disentangle the impacts from different processes, particularly when the plant recoveries are non-random in relation to species ranges (i.e. faster plant recovery at low elevations may contribute to the observed lag in the responses of plant communities to climate warming (Lenoir et al 2010b, Bertrand et al 2011). It is noteworthy that although the downslope shift was a general pattern, there were large differences between individual species in their elevational shifts. For example, in the eight species with significant elevation coefficients (labeled species in figure 5(B)), one species shifted upslope, five species shifted downslope, and two species did not show clear or consistent shifts. We offer three possible explanations for these speciesspecific responses: (1) some species are more responsive to warming climate, and thus are more likely to shift upslope; (2) some species are more responsive to the cessation of disturbances, and thus are more likely to shift downslope as low elevations are allowed to recover; (3) other factors that were not accounted for in this study play an important role in determining species shifts. In this study region, the upslope shifting species, Platycarya strobilacea occurs at relatively high elevations (400-1400 m) (Lu et al 2002), and may be more sensitive to the warming climate that alleviate cold stress in high elevations. In contrast, the downslope-shifting species are dominated by oaks (Q. chenii, Q. variabilis and Q. glauca). Q. chenii occurs at lower elevations (below 600 m) (Huang et al 2002), and thus may be more sensitive to the disturbance cessation at low elevations. The other two oak species occur throughout the elevational ranges of our study (Huang et al 2002), and thus their downslope shifts may be affected by other processes. For example, the oak species are relatively adapted to dry sites (Battaglia et al 2000, Johnson et al 2009), which tend to be more common at low elevations (Sundqvist et al 2013). Moreover, downslope shifts of plant species could be related to their strategies of seed dispersal (Cannone and Pignatti 2014), since oak species have relatively long dispersal distances (e.g. more than 150 m (Johnson et al 2009)) which may facilitate their ability to capitalize on new habitat available following the cessation of disturbance.
In addition to the differences in shifting direction, the magnitude of downslope shifts was generally higher than the upslope ones (i.e. the magnitude of downslope shift reached ∼150 m, but upslope shift was ∼50 m in figure 5(B)). Faster downslope shifts might provide additional support for the fact that downslope shifts were related to the forest recovery following disturbance cessation, since the high magnitudes of downslope shift are beyond what would be expected from climate change over the 20 years. These species-specific responses should be considered in future studies to enhance our understanding of how species will shift under the joint impacts of climate and land-use changes

Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request. The data are not publicly available for legal and/or ethical reasons.  (8)). (B) The relationship between species elevational shifts (SESs) and coefficients of elevation (b i 1, in the equation (9)). The labeled specie had the significant coefficients of elevation (p<0.05). The different sizes of the scatter points indicated the variances of SESs between the five study sites for each species, i.e. larger sizes indicated smaller variance (See details in the Method 2.5.4). The regression equation (its line and standard errors), p value and r 2 value were based on the equation (10).  (6) in the Method 2.5.2. Figure A4. Correlogram of equation (7) in the Method 2.5.2.