Contrasting plant responses to multivariate environmental variations among species with divergent elevation shifts.

The general predictions of climate impacts on species shifts (e.g., upward shift) cannot directly inform local species conservation, because local-scale studies find divergent patterns instead of a general one. For example, our previous study found three shift patterns with elevation (strong down-, moderate down-, and up-slope shifts) in temperate mountain forests. The divergent shifts are hypothesized to arise from both multivariate environmental variations with elevation and corresponding species-specific responses. To test this hypothesis, we sampled soils and leaves to measure elevation variations in soil conditions and determined plant responses using discriminations against heavier isotopes, carbon (13 C) and nitrogen (15 N). Functional traits of the species studied were also extracted from a public trait dataset. We found that: (1) With low soil water contents at low elevations, only the leaves of up-shifters had lower 13 C discriminations at low vs. high elevations; (2) With low soil P contents at high elevations, only the leaves of moderate down-shifters had higher 15 N discriminations at high vs. low elevations; (3) The leaves of strong down-shifters did not show significant elevation patterns of the discriminations; (4) The contrasting responses among the three types of shifters agree with their functional dissimilarity, suggested by their separate locations in a multitrait space. Taken together, the divergent shifts are associated with the elevation variations in environmental conditions and contrasting plant responses. The contrasting responses could result from the functional dissimilarity among species. Therefore, a detailed understanding of both local environmental variations and species-specific responses can facilitate accurate predictions of species shifts to inform local species conservation.

Introduction

Some biogeographical studies showed general patterns of species range shift (upward and poleward shift) associated with warming temperatures (Parmesan and Yohe 2003, Lenoir and Svenning 2015), but other studies have highlighted that temperature is only related to range limits of no more than one‐third of species studied (Normand et al. 2009) and even weakly related to observed range limit shifts in Rumpf et al. (2019). Moreover, local studies observed shifts of different directions and rates, e.g., multiple shift patterns with elevation (Guo and Lenoir 2018). Our previous study also showed divergent shift patterns across woody species along the same elevation gradient (Zhang et al. 2019). In that study, we found some oak species moved down‐slope as forest recovered from disturbance, despite warming temperatures. Therefore, the divergent shift patterns at local scale could be driven by other factors in addition to temperature, and understanding of these potential factors is a crucial knowledge gap in the study of species shifts under changing climate.

The divergent shift patterns could result from multivariate nature of a changing environment (Kling et al. 2020), including nutrient and water. Nutrient limitation is a fundamental limitation of primary producers (Elser et al. 2007, Vitousek et al. 2010), therefore fertilization is used to mitigate effects of nutrient limitation on forest growth, e.g., pine plantation (Zhai et al. 2015). However, nutrient limitation has not received sufficient attention regarding its impacts on species shifts (Fajardo and Piper 2017). Nutrient limitation is potentially important, because plant nutrient status could determine both plant responses to changing climate, e.g., drought (Gessler and Schaub 2017), and predictions of species distributions (Buri et al. 2017). The variation in nutrient distribution can be found along elevation gradients in mountains, e.g., soil nutrient contents could be lower at higher elevations (Wang 2018). Therefore, decreasing foliar nutrient contents (foliar %N and %P) with elevation have been found in some studies, e.g., Asner (2017), but this elevation pattern cannot be applied globally (Körner 2007). The low nutrient availability at high elevations further constrains plant growth in Alpine ecosystems where there are early‐successional soils (Hagedorn and Gavazov 2019). Moreover, drought‐driven forest shifts were widely observed (Fei et al. 2017). In mountain systems, moisture availability can mediate warming‐caused up‐slope shifts of treeline (Sigdel et al. 2018). Water availability could be reduced at lower elevations (Normand et al. 2009), and variation in water availability account for plant shifts in the mountains of California, USA (Rapacciuolo 2014). In temperate mountain forests of our study areas, the water availability at low elevations might be reduced under a warming climate and consistent rainfall pattern (Zhang et al. 2019). Therefore, the reduced water availability was hypothesized to lead to up‐slope shifts of some woody species (Fig. 1). Down‐slope shifts of the other woody species are hypothesized to arise from lowland forest recovery with more growth space available following cessation of historical harvests, particularly when nutrients are limiting at high elevations (Fig. 1). Although these hypothesized impacts of nutrient, water, and growth space might agree with the field conditions in our study area, there is a lack of rigorous hypothesis tests, leading to a limited understanding of the impacts of multivariate environment on plant elevation shifts.

Fig. 1

Conceptual figure of the elevation patterns (from 60 to 1,150 m) of environmental conditions (soil water and nutrient contents and growth space) and the associated species shift patterns. The elevation change in soil water is hypothesized to drive up‐slope shifts, and the elevation changes in soil nutrient and lowland growth space are hypothesized to drive down‐slope shifts. Note: Seed‐related processes (seed dispersal and germination) and parent tree effect were also considered in this study to understand the shift patterns.

Plant responses to the elevation changes in water and nutrient conditions can be determined by discrimination against heavier carbon (13C) and nitrogen (15N) isotopes, respectively, from their sources to leaf materials in conjunction with other measurements (e.g., foliar nutrient contents and soil isotope signatures). Atmospheric CO2 diffuses into air‐filled spaces within leaves and is fixed by photosynthesis, and both the diffusion and biochemical processes discriminate against 13C in favor of 12C. The discrimination against 13C within leaves declines with increasing water stress, because it reduces atmospheric CO2 diffusion into leaves by closing stomata. Then, a greater proportion of 13CO2 within leaves is fixed during photosynthesis (Farquhar and Sharkey 1982). Therefore, lower 13C discrimination is observed as a plant experiences more limited water availability, e.g., at low elevations (Lajtha and Getz 1993). Similarly, plant nitrogen uptake discriminates against 15N in favor of 14N, and the discrimination increases with lower plant N demand or root growth/activity (Fry et al. 2000, Robinson et al. 2000). The low N demand could arise from plant P limitation based on optimum N:P ratios for plant growth (Ågren and Wetterstedt 2012). Therefore, higher 15N discrimination is observed as plant response to a more limited P availability (Inglett et al. 2007, Zhai et al. 2018). Together, our study measured the 13C and 15N discriminations to indicate plant responses to the elevation changes in soil water and nutrient contents, respectively.

The divergent shift patterns could also arise from different responses to environmental changes among species (Kling et al. 2020), and the different responses could be caused by species‐specific functional traits, affecting plant growth with changing environment (Violle et al. 2007). For example, plant responses to water and nutrient limitations can be defined by leaf thickness (LT) and specific leaf area (SLA) (Griffin‐Nolan 2018, He et al. 2018), and leaf dry matter content (LDMC) and leaf size (Hodgson 2011), respectively. Based on these functional traits, species showed different probabilities of colonizing different elevations (Matteodo et al. 2013). Therefore, species with divergent shift patterns may arise from their functional dissimilarities, determined by a multivariate analysis of functional traits (Stahl et al. 2013). The functional dissimilarities can be visualized by separate locations of species in a multidimensional space, defined by plant functional traits. To test the functional dissimilarities, our study extracted functional trait data of the species studied from open‐access plant trait databases. These public data have been widely used similarly by other studies quantifying functional attributes of local species, e.g., quantifying functional diversity (Ouyang 2019).

To test if the divergent shift patterns were attributed to the multivariate environmental variations and species‐specific plant responses, we measured the elevation variations in soil conditions (soil water and nutrient contents), and we also determined the plant isotope responses with leaf samples from eight woody species characterized with divergent shift patterns (strong down‐, moderate down‐ and up‐shifters). Furthermore, to understand the species‐specific responses, we evaluated the functional dissimilarities among the species studies using their functional traits from a trait database. In detail, we tested the following hypotheses: (1) Soil water content increases with elevation, but soil nutrient contents decrease with elevation; (2) The 13C discrimination of up‐shifters and 15N discrimination of down‐shifters increases with elevation; (3) Strong down‐shifters have higher growth rates than moderate ones; (4) There are functional dissimilarities among the three types of shifters.

Materials and Methods

Study area

The study area was located in mountain forests in Liu'An of Anhui province, eastern China (115°22′–116°11′E longitude and 31°06′–31°48′N latitude). Based on monthly climate records between 1996 and 2014, this area receives a minimum 28.6 ± 1.1 mm, maximum 253.3 ± 88.5 mm, and average 90.8 ± 77.9 mm precipitation, and a minimum 2.7 ± 1.3°C, maximum 30.8 ± 7.0°C, and average 16.4 ± 8.8°C. From 1996 to 2014, the mean temperature for the study area increased gradually from 16.0°C to17.7°C, but no significant change in precipitation has been recorded in this area from 1996 to 2014 (Zhang et al. 2019).

There were 89 permanent forest plots (25.82 × 25.82 m) established at regularly spaced sites along an elevation gradient from 60 to 1,150 m. Before 1970, this area was harvested for timber and fuel woods. The lowland plots suffered more from the harvests (i.e., a higher number of stems harvested) compared with ones located at higher elevations, but the harvest and natural disturbances (e.g., fire) have been strictly limited by local government agencies since 1970, and this conservation policy contributed to the natural forest recovery in our study site (Zhang et al. 2019). Beginning in 1989 and every five years through 2009 (1989, 1994, 1999, 2004, and 2009), in total 3,316 living woody stems ≥5 cm in diameter at breast height (DBH) of 18 woody species were tagged, measured (DBH), recorded according to growth status (live or dead) and identified to species level.

Sample collection and measurements

Woody species

In our study site (i.e., Liu'An), the species shift directions and rates were quantified by changes in species mean elevations, that is the mean elevation weighted by the number of individuals per elevation in each census. The calculation details can be found in Zhang et al. (2019). We sampled eight woody species, showing three shift patterns: (1) Strong down‐shifters (shift rate < −5 m/yr): Dalbergia hupeana, Cerasus pseudocerasus, and Quercus glauca; (2) Moderate down‐shifters (−5 m/yr ≤ shift rate ≤ 0 m/yr): Pinus massoniana, Lindera glauca, and Cunninghamia lanceolata; (3) Up‐shifters (shift rate > 0 m/yr): Pinus taiwanensis and Liquidambar formosana. The eight species are dominant ones, collectively representing >80% of the individuals in our study site. In addition to this study site (Liu'An), our previous study reported the different shift patterns among species with data from Liu'An and the other four study sites in the same region (Zhang et al. 2019). This study aimed to investigate mechanisms of the different shift patterns among species along the same elevation gradient.

Plots

In September 2016, we collected samples in 16 plots from the 89 plots established in 1989. These plots were uniformly distributed along the same slope from 310 to 1,150 m, leading to one plot per 52 m with elevation (Appendix S1: Fig. S1). Moreover, these plots were set up with similar aspects to control the effects of aspect. In each sampling plot, we chose three random subplots (5 × 5 m), which lay along a diagonal line crossing the sampling plot. In each subplot, two soil samples were collected from the surface soil (0–20 cm) and subsurface soil (20–40 cm) (after forest floor litter was removed). Soil samples were sealed by preservative film and cellulose tape and individually stored in plastic zip bags to avoid evaporation. In addition, in each sampling plot, ≥10 fully expanded and sun‐exposed leaves were sampled from the upper canopy of an individual stem, and three stems were chosen for each of the eight selected species. The leaf samples were individually stored in plastic bags wrapped with wet gauze. These samples were collected by a single team of researchers using standardized collection protocols.

Measurements of soil and plant leaves

To measure gravimetric soil moisture, we removed stones and roots in the soil core sample by sieving through a coarse sieve. Then, we calculated soil moisture by the difference between fresh and dry mass (dried at 105°C until a constant mass), and the difference was standardized with dry mass. Remaining soil samples were air‐dried (at ˜15°C), sieved (2 mm mesh), and stored in air‐tight plastic bags to measure: (1) soil total nitrogen (TN, g/kg) by an elemental analyzer (Vario EL III, ELementar, Germany); (2) total phosphorus (TP, g/kg) using a UV–visible light spectrophotometer.

Foliar and soil samples were oven dried (60°C for at least 48 h) and ground. We loaded both soil and foliar samples (≥5 mg) into individual tin cups to measure N content (%N) and stable isotope ratios. The stable isotope measurements were made by a continuous flow isotope ratio mass spectrometer (Costech Analytical Technologies, USA) at the isotope analysis laboratory, Nanjing Forestry University. The stable isotope ratios of carbon and nitrogen were reported using δ13C and δ15N, respectively, applying Eq. 1:

The isotope discriminations against 13C and 15N relative to their respective sources (atmospheric CO2 and soil N) were calculated using Eq. 2: where y is the atomic mass of the heavy isotope (C: 13 and N: 15), X is the atom of interest (C and N), and are the measured ratios of the sample and a reference standard. δ13C values are expressed relative to the Vienna Belemnite of the Pee Dee formation and δ15N values are expressed relative to standard air nitrogen. is the isotope ratio of sources, and is the isotope ratio of leaves. The δ13C of atmospheric CO2 of our study site (−8.31‰) was extracted from nearby stations affiliated to the World Data Centre for Greenhouse Gases, World Meteorological Organization (WMO) (https://gaw.kishou.go.jp/). The elevation pattern of δ13C of atmospheric CO2 was assumed to be constant in our study and others, e.g., Wu et al. (2017), because it had a marginal effect on foliar δ13C (Körner and Farquhar 1988). There was a steady decrease in the atmospheric CO2 partial pressure at a rate of 9.5% per 1,000 m elevation (McElwain 2004), but the decreasing partial pressures only caused minor changes in the 13C discrimination. For example, the effect of CO2 partial pressure was greatest in the plots at the highest elevation where the discrimination was only reduced by less than 10% of the original value (Seibt et al. 2008). Moreover, at the same elevation, different species are exposed to the same CO2 partial pressure. Together, we assumed that CO2 partial pressure does not affect our comparison among the species studied. For soil samples, we measured the δ15N of N content.

Functional traits

From the 18 woody species in our study area, functional trait data were extracted based on 11 species: (1) strong down‐shifters: Quercus chenii, Q. glauca, D. hupeana, and Pistacia chinensis; (2) moderate down‐shifters: P. massoniana, Cunninghamia lanceolate, L. glauca, and Diospyros lotus; and (3) up‐shifters: L. formosana, Litsea cubeba, and L. glauca. The species selection was based on data availability from the TRY plant trait database (Kattge 2020). Our functional trait analysis was based on the following functional traits: SLA, leaf area (LA), leaf phenology (LP), LT, LDMC, seed dry mass (SDM), and maximum height (HT). For each of these traits, the mean value of one species was calculated by all observations found from the trait database. The maximum heights were compiled from Wu and Raven (1999). Additionally, the growth rate of basal area (BA.Growth) was added to the functional trait analysis. BA.Growth was calculated by mean growth rate across the measured individuals from 1989 to 2009. Moreover, the BA.Growth was standardized with initial BA, as a way to account for the effects of tree size associated with differences in elevation‐related productivity and tree age. The calculation details can be found in Zhang et al. (2019). To summarize the information above, we made a table to describe details of the data used for the functional trait analysis, including species name, shift group, functional traits, etc. (Appendix S1: Table S1).

Statistical analysis

Elevation patterns of soil and foliar attributes

Both of the elevation patterns of soil (water content and N:P) and foliar (13C and 15N discriminations) attributes were analyzed by generalized least squares (GLS) and linear mixed‐effects (LME) models, respectively (Appendix S1: Section S1). The residual structures of these models may vary by different dependent variables to allow for the spatial correlation with elevation and the heterogeneity in model residuals (Appendix S1: Section S1). Akaike information criterion was used to select models with different structures of the fixed effects and model residuals. Maximum likelihood was used to compare models with different fixed‐effects structures (e.g., with or without the interactions of soil depth and elevation in the models of soil attributes). Restricted maximum likelihood was used to compare models with different residual structures and estimate model parameters. We removed influential points and outliers (<5% total observations) identified by Cook's distance and standardized residuals, respectively. The violations of linearity, equal variance, normality, and independence (i.e. spatial independence with elevation) were examined by the residual plots along the fixed‐effects factors and fitted values, the normal quantile plot, and the variogram of standardized residuals, respectively. The above statistical analyses were implemented in the R program (R Core Team 2020) and the nlme package (Pinheiro et al. 2019).

Relationship between foliar 15N discriminations and foliar nitrogen content

This relationship was calculate using Model II regression, because both foliar 15N discriminations and foliar nitrogen content were measured values with random error (Sokal and Rohlf 1995). As the distributions of these two variables were not bivariate normal and their relationship was linear, Model II regression was calculated using the ordinary least squares method (Legendre and Legendre 1998). Foliar nitrogen content on the base‐10 logarithmic scale was used in the analysis. The above statistical analysis was implemented using R and the lmodel2 package (Legendre and Oksanen 2018).

Functional trait analysis

Functional traits were analyzed with a multivariate analysis using non‐metric multidimensional scaling (NMDS). NMDS was performed on a matrix of Euclidean distances among the nine species. Two dimensions were used in the NMDS, as adding more dimensions did not significantly reduce the stress value. The stress value of the NMDS was 0.06. The above statistical analysis was implemented in R and the vegan package (Oksanen 2020).

Results

Soil water content and 13C discrimination

There were lower soil water contents at low vs. high elevations (Fig. 2a, P < 0.01). With the low soil water contents at low elevations, there were no significant elevation patterns of 13C discrimination for the strong (Fig. 2b) and moderate down‐shifters (Fig. 2c). However, the up‐shifters had lower 13C discrimination at low vs. high elevations (Fig. 2d), suggesting their significant response to the lowland water limitation.

Fig. 2

The relationships between elevation (m) and (a) soil water content; (b) 13C discrimination from air to leaf of the strong down‐shifters; (c) 13C discrimination from air to leaf of the moderate down‐shifters; and (d) 13C discrimination from air to leaf of the up‐shifters. The relationships of elevation with soil water content and 13C discrimination are denoted by the fitted lines based on the population means of the Eqs. 1 and 2 from Appendix S1: Section S1, respectively. The solid lines denote significant elevation effects at a significance level of 0.05. The significances of the elevation effect are denoted by the P‐values based on the above equations. One point in the figures is one measurement of an individual sample, and its color is coded by genera.

Soil N:P, 15N discrimination, and foliar nitrogen

Soil N:P ratios increased with elevation (Fig. 3a), suggesting high P limitations at high elevations. With high P limitations at high elevations, neither the strong down‐shifters (Fig. 3b) nor up‐shifters (Fig. 3d) showed a significant elevation pattern for the 15N discrimination. However, the moderate down‐shifters had greater 15N discrimination at high vs. low elevations (Fig. 3c), indicating their significant nutrient response. Furthermore, the significant nutrient response of moderate down‐shifters was suggested by their greater 15N discriminations with lower foliar %N, i.e., the negative correlation in Fig. 4a.

Fig. 3

The relationships between elevation (m) and (a) natural logarithm of soil nitrogen and phosphorus ratio; (b) 15N discrimination from soil to leaf of the strong down‐shifters; (c) 15N discrimination from soil to leaf of the moderate down‐shifters; and (d) 15N discrimination from soil to leaf of the up‐shifters. The relationships of elevation with soil N:P and 15N discrimination are denoted by the fitted lines based on the population means of Eqs. 3 and 4, respectively, reported in Appendix S1: Section S1. The solid lines denote significant elevation effects at a significance level of 0.05. The significances of the elevation effect are denoted by the P‐values based on the above equations. One point in the figures is one measurement of an individual sample, and its color is coded by genera.

Fig. 4

The relationships of the moderate down‐shifters between the foliar 15N discrimination and foliar nitrogen content at the base‐10 logarithmic scale. The relationship is denoted by the fitted line that was based on the linear II regression analysis. The significance of the relationship is denoted by the P‐value. One point in the figures is one measurement of an individual sample.

Plant functional traits

The two‐dimensional space was downscaled from the multidimensional space built by the functional traits (Fig. 5). Axis 1 was positively related to SLA, but negatively related to LT, LP (from deciduous to evergreen species), and BA.Growth. The second axis was positively related to LDMC and SDM, but negatively related to LA and HT. In this trait space, the species showing divergent shift patterns were grouped separately. Among the three shifter groups, the up‐shifters had the highest LA and SLA (axis 1), the strong down‐shifters had the highest LDMC and SDM (axis 2), and the moderate down‐shifters had the highest BA.Growth (axis 1), greatest HT (axis 2), and were most dominated by evergreen species (axis 1) (Fig. 5).

Fig. 5

Two‐dimensional NMDS ordination diagram with the three groups of shifters (including strong down‐shifters, moderate down‐shifters, and up‐shifters) and three ellipses showing the 0.95 confidence regions for the locations of group centroids. BA.Growth, growth rate of the basal area; HT, maximum tree height; LA, leaf area; LDMC, leaf dry matter content; LP, leaf phenology; LT, leaf thickness; SDM, seed dry mass; SLA, specific leaf area.

Discussion

Up‐slope shifts are associated with low soil water contents at the low elevations

Only the up‐shifters showed significantly lower 13C discriminations at the low vs. high elevations. The isotope response aligns with the low soil water content at the low elevations in our study area. Other factors may confound the ability of 13C discrimination to estimate the plant response to soil water variation with elevation, such as temperature‐related changes in leaf traits (Vitousek and Field 1990, Maxwell and Silva 2018), but these changes contributed to higher 13C discrimination at low vs. high elevations, which is opposite to our results. Therefore, these confounding influences tend to be offset by lowland water limitation (Van de Water and Leavitt 2002), leading to declining 13C discriminations at low elevations, e.g., Li et al. (2009), which were also observed from the up‐shifters in our study. The other confounding factor is the anisohydric degree of these species. Although direct measurements of anisohydric degree are not possible, we compiled published studies in Appendix S1: Table S2, suggesting that most of the species studied reduce stomatal conductance or 13C discrimination under limited water availability. Therefore, there are minimal anisohydric effects on the species studied, and water availability limitation provided an important constraint on the lower ends of plant distribution (Colwell et al. 2008, Rigling 2013), but only for the up‐shifters in our study.

Moderate down‐slope shifts are associated with low soil nutrient contents at the high elevations

The significantly increasing soil N:P ratios with elevation in our study area suggested P limitation of the high elevations (Elser et al. 2010). P limitation could be related to faster N deposition (Lovett and Kinsman 1990) or different mineralization rates between N and P (Wang et al. 2018) at high elevations. Importantly, P limitation could reduce plant N demand contributing to lower N uptake, further increasing 15N discriminations. Additionally, the lower N uptake and resultant discrimination may arise from reduced root activity, due to the decreasing water viscosity and cell membrane permeability in a cold climate characterizing high elevations (Mayor 2017). The higher discrimination with lower N uptake is further suggested, as foliar %N declined with higher discrimination in our study and others, e.g., Brookshire et al. (2020). Potential confounding factors on 15N discrimination is a switch of plant N source from inorganic to organic forms of N by declining temperature with elevation (Schimel and Bennett 2004), but the switch may lead to lower discriminations (Averill and Finzi 2011), instead of the higher discrimination found in our study. Therefore, the confounding effects are marginal on the species comparisons in our study.

As P limitation can constrain plant growth by decreasing investments in P‐rich RNA required for organismal growth (Elser et al. 2010), it constrained the up‐slope shifts of responsive species, i.e., the moderate down‐slope shifting species in our study. Notably, highland nutrient stress varies among distantly geographical regions, e.g., highland P limitation in northeastern China (Zhao et al. 2016), but N limitation in the Peruvian Andes (Fisher et al. 2013). These highland nutrient limitations are further supported by field experiments in which fertilization stimulated highland forest growth (Möhl 2019). Therefore, highland nutrient limitation could determine upper limits to plant distributions (Rehm and Feeley 2015). In our study, with the limited potential of up‐slope shifts by the highland nutrient stress, species may leverage the lowland growth space available from the cessation of historical harvests, contributing to their down‐slope shifts.

Three types of shifters were grouped separately in the trait space

Compared with the up‐shifters and moderate down‐shifters, the strong down‐shifters did not show the significant isotope responses with elevation. This finding might suggest their tolerance to the variations in soil water and nutrient with elevation. Notably, compared with the other shifters, the tolerances of strong down‐shifters to the lowland water and high nutrient limitation were also supported by their relatively low SLA and high LDMC, respectively (Peìrez‐Harguindeguy 2013). With the tolerances, the strong down‐shifters could take advantage of the lowland growth space following the cessation of historical harvests in our study site (Zhang et al. 2019), and the harvest cessation‐caused down‐slope shifts of forests are also found in the mountains of the northeastern USA (Wason and Dovciak 2017). Therefore, post‐harvest regeneration or land‐use history should be considered to predict plant‐distribution shifts (Miller and McGill 2018, Wang et al. 2019). Compared with strong down‐shifters, the moderate down‐shifters had faster growth rates. Therefore, instead of growth rate, other mechanisms could explain their differences in the shift rate. One mechanism might be related to seed‐related processes, including seed dispersal and germination. A long dispersal can facilitate elevation shifts and space colonization (Cannone and Pignatti 2014), characterizing the dominant species of strong down‐shifters, e.g., oak species (Q. chenii, Q. glauca) (Johnson and Shifley 2009). Dispersal distance tends to be negatively related to plant height (LaRue and Holland 2018), therefore the shorter heights of strong down‐shifters could contribute to their long dispersal compared with the moderate down‐shifters. Moreover, seed germination initiates plant establishment, potentially affecting plant shifts. Seed germination is positively affected by SDM (e.g., Stanton 1984) which was higher on the strong down‐shifters compared with moderate ones.

The three types of shifters were grouped separately in the trait space. Axis 1 of the space is related to traits more responsive to water stress, such as LT and SLA. A similar dimension is also found in another study analyzing intraspecific trait variations as they change with elevation (Umaña and Swenson 2019). The other axis is related to traits more responsive to nutrient stress, such as LDMC (Hodgson et al. 2011). In the trait space, the separate locations of these species groups suggested that their functional dissimilarities may account for the divergent shift patterns among species along the same elevation gradient. For example, based on the dissimilarity in leaf traits suggesting tolerance to water stress (e.g., SLA, LT, LA), the down‐shifters showed higher tolerance than the up‐shifters, and the tolerance differences aligned with the lower soil water content in low elevations compared with high ones. Therefore, the general shift patterns (e.g., upward shift) with warming climate could be mediated by species‐specific plant responses to environmental variations. These species‐specific responses have also been found in other elevation‐based studies. For example, water stress has stronger effects on highland vs. lowland species from the Bavarian Alps in Germany (Rosbakh et al. 2017). In addition to the functional differences, the different shift patterns might be related to remaining trees after the harvest was ceased, particularly for the species of strong down‐shifters in the low elevations, because these trees could be parent trees contributing to the downshift. However, the species of strong down‐shifters made up a small percentage of trees in the low elevations in the beginning, based on the survey data in 1989 (Appendix S1: Table S3). Therefore, the parent tree effect cannot be supported by the forest survey data.

Conclusion and management implications

Our study suggests that the divergent shifts among co‐occurring species along the same elevation gradients could be attributed to both multivariate environmental variations and species‐specific plant responses. Multivariate environmental variations are characterized with elevation variations in soil water, soil nutrient, and growth space in our study. The species‐specific responses arose from the functional dissimilarities. Therefore, integrating information on environmental variation and plant functional traits could inform predictions of species shifts (e.g., Zhai et al. (2016)). The trait‐based predictions of species distribution have already been applied on a global scale. For example, van Bodegom and Douma (2014) derived global trait maps based on the association between traits and multiple environmental drivers. This study also characterized nine globally representative vegetation types based on their trait combinations. Importantly, valid distribution predictions of these vegetation types were made by the trait maps. The potential of trait‐based prediction for local conservation is well supported by the findings from our study. Notably, prediction at a local scale can use plant traits to directly characterize species of interest, instead of the artificial vegetation types used on the global scale.

Appendix S1

Click here for additional data file.

Data S1

Click here for additional data file.