Long-term monitoring reveals forest tree community change driven by atmospheric sulphate pollution and contemporary climate change

AIM: Montane environments are sentinels of global change, providing unique opportunities to assess impacts on species diversity. Multiple anthropogenic stressors such as climate change and atmospheric pollution may act concurrently or synergistically in restructuring communities. Thus, a major challenge for conservation is untangling the relative importance of different stressors. Here, we combine long‐term monitoring with multivariate community modelling to estimate the anthropogenic drivers shaping forest tree diversity along an elevational gradient. LOCATION: Camels Hump Mountain, Vermont, USA. METHODS: We used Generalized Dissimilarity Modelling (GDM) to model spatial and temporal turnover in beta diversity along an elevational gradient over a 50‐year period and tested for spatiotemporal shifts in density and elevational distribution of individual species. GDMs were used to predict community turnover as nonlinear functions of changes in elevation, climate and atmospheric pollution. RESULTS: We observed significant shifts in elevational range and density of individual species, which contributed to an overall reduction in the elevational gradient in beta diversity through time. GDMs showed the combined effects of sulphate deposition and temperature as drivers of this temporal reduction in beta diversity. Spatiotemporal changes differed among species, with shifts observed both up‐ and downslope. For example, in a reversal of a previous upslope range contraction, red spruce (Picea rubens Sarg.) increased in density and shifted downslope since the 1990s, occupying warmer, drier climates. MAIN CONCLUSION: Our results demonstrate that global change is impacting the stratification of forest tree diversity along elevational gradients, but the responses of individual species are complex and variable in direction. We suggest abiotic drivers may directly impact individual species while also indirectly altering species interactions along elevational gradients. Our approach modelling the drivers of compositional turnover quantifies the rate and amount of change in beta diversity along environmental gradients and serves as a powerful complement to studying species‐specific responses.


| INTRODUC TI ON
Understanding the impacts of contemporary global change on the diversity of forest tree communities represents a critical challenge for conservation biologists, as trees are foundational components of terrestrial ecosystems that provide extensive ecological and economic value (Ellison et al., 2005;Myneni et al., 2001). Both forest tree community composition and the distribution of individual species are strongly tied to climatic conditions through species-specific physiological tolerances to the abiotic environment, including temperature and water availability (Grubb, 1977;Woodward & Williams, 1987). As such, anthropogenically induced climate changes can exert strong influences on forest tree distribution and the biodiversity of forest communities. Climate-driven range shifts have been documented through poleward and upslope migrations in response to warming temperatures (Parmesan, 2006;Parmesan & Yohe, 2003).
However, poleward range shifts depend upon dispersal over large geographic distances. In contrast, montane forest communities are more sensitive to global change because steep environmental gradients compress variation in climate and other stressors that occur over relatively short spatial scales (e.g., Beckage et al., 2008;Kelly & Goulden, 2008;Walther, Beißner, & Burga, 2005).
While plant taxa have demonstrated the ability to track suitable climatic conditions (e.g., Martínez-Meyer & Peterson, 2006), species responses are often a manifestation of interactions with other agents including atmospheric pollution, land use history and other biota (Bytnerowicz, Omasa, & Paoletti, 2007;Peltzer, Allen, Lovett, Whitehead, & Wardle, 2010;Tylianakis, Didham, Bascompte, & Wardle, 2008;Walther, 2010). For example, in the north-eastern United States atmospheric pollution has become a critical environmental stress in forested ecosystems as a result of anthropogenic activities dramatically increasing nitrogen and sulphate deposition compared to pre-industrial levels (Clark et al., 2018;Driscoll et al., 2001). Elevated acidic atmospheric deposition increases the leaching of base plant nutrient cations and mobilizes aluminium (Al) in the soil of forested ecosystems, which results in long-lasting consequences of soil acidification and a cation imbalance in plant photosynthetic tissue (Driscoll et al., 2001). Consequently, atmospheric pollution can compound the effects of climate change to alter tree physiology and the competitive regime within the forest tree community, causing shifts in species composition that may not strictly correspond to prevailing climate change predictions.
In the Appalachian Mountains of the north-eastern United States, empirical studies documenting elevational range shifts have reported conflicting findings. For example, the temperateboreal ecotone (TBE)-the transition zone from lower elevation hardwood forests to higher elevation spruce-fir coniferous forests-occurs over a relatively small spatial scale (e.g. ~100 m) in which community composition changes rapidly with elevation relative to higher or lower elevations. The location of this transition is thought to be largely driven by climatic factors such as mean summer temperature (Allen & Breshears, 1998;Cogbill & White, 1991), yet the location of the TBE has been shown to have moved both up-and downslope across the region, suggesting drivers other than climate alone are influencing the location of the ecotone Foster & D'Amato, 2015).
Beginning in the mid-1960s, P. rubens showed a regional decline in the northeast (Johnson, Cook, & Siccama, 1988;Siccama, Bliss, & Vogelmann, 1982;Vogelmann, Perkins, Badger, & Klein, 1988) attributed to winter foliar injury linked to acid deposition-induced calcium (Ca) depletion (DeHayes, Schaberg, Hawley, & Strimbeck, 1999;Schaberg & DeHayes, 2000). While recent reports suggest P. rubens is making a recovery (Battles, Fahey, Siccama, & Johnson, 2003;Kosiba, Schaberg, Rayback, & Hawley, 2018;, it is unknown how climate and/or atmospheric deposition have impacted forest tree community composition as a whole across the elevational gradient encompassing the ecotone. Investigating patterns of biodiversity along environmental gradients often focuses on modelling either individual species or higher-level entities in aggregate, such as ecosystems, vegetation functional types or communities Ferrier, 2002;). An alternative approach involves modelling emergent properties of these communities as they change, or "turnover," along geographical and environmental gradients (Ferrier, Manion, Elith, & Richardson, 2007). Turnover is related to the concept of beta diversity, which quantifies differences in community composition between samples. Turnover can be thought of as the rate of change in beta diversity through either space or time (e.g., Blois, Williams, Fitzpatrick, Jackson, & Ferrier, 2013) and can be modelled as pairwise changes in species composition (e.g. dissimilarity) among sampling points as a environmental gradients and serves as a powerful complement to studying speciesspecific responses.

K E Y W O R D S
beta diversity, elevational gradients, generalized dissimilarity modelling, global change, montane environments, species distributions function of their differences along one or more environmental gradients. As such, high turnover suggests community diversity is structured through space or time by one or more environmental drivers, while low turnover is indicative of a more homogenous or stable community. A benefit of modelling beta diversity is that it captures changes in the diversity of the system as a whole, including all available data across all species, regardless of the number of records per species, and including the contribution from rare species that may be difficult to model individually. Moreover, this approach does not assume species will respond in a cohesive way to environmental and climatic change. Thus, modelling turnover along spatially and temporally changing environmental gradients can quantify emergent properties of biological responses to global change and serves as a powerful complement to studying species-specific responses.
Forest tree community responses to global change are complex and extend over long temporal periods as a result of slow migration rates and the long life span of trees. While anthropogenic stressors are known to impact forest tree species and communities, estimating the relative importance of different stressors driving community change is critical for the conservation of beta diversity. In this study, we couple multivariate modelling of biodiversity turnover with data from long-term monitoring of forest communities along an elevational gradient in the north-eastern United States. Specifically, we seek to (a) characterize how the elevational gradient in forest composition has shifted over a 50-year period, both in terms of beta diversity and the distribution and density of individual species, and (b) determine the environmental drivers contributing to temporal changes in forest tree community composition and how these drivers may be differentially impacting low versus high-elevation forest communities. Collectively, the use of a multivariate approach with long-term inventory data allowed us to unravel the drivers of beta diversity turnover along spatial and temporal environmental gradients in north-eastern forests.

| Study area and monitoring study
In north-eastern North America, montane forest communities are structured along an elevational gradient in three relatively distinct zones (Siccama, 1974;Figure 1a In 1964, a long-term forest monitoring study was established on Camels Hump, a 1,248 m peak in the Green Mountains of Vermont (44°19′N, 72°53′W). Prior to the mid-1950s, there was selective logging in the low-elevation hardwood forests, but cutting was absent in the study area (Siccama, 1974;Vogelmann, Badger, Bliss, & Klein, 1985;Whitney, 1988). Inventory stands were established at intervals of approximately 60 m along an elevational gradient from 549 to 1,159 m on the mountain's western slope (Siccama, 1974). Stands were systematically established to have an unbiased sampling of the plant distribution along the elevation transect. At each stand (N = 11), five to ten permanently marked survey plots, each 3.0 × 30.5 m, were placed parallel to each other and roughly perpendicular to the contour line (separated by approximately 6 m), so that each plot has a unique elevation value. There are ten plots per stand below 884 m and five plots per stand above this elevation to maintain a westerly aspect. Within each plot (N = 85), species identity and diameter at breast height (dbh; 1.37 m above ground) of all woody stems > 2 cm dbh were recorded during each of nine census years (1965, 1979, 1983, 1986, 1990, 1995, 2000, 2004 and 2015).

| Environmental data
Environmental predictors of biodiversity turnover were chosen based on their known or predicted influence on forest tree distributions in north-eastern montane ecosystems. High-resolution climate data are not available at the fine spatial scale of our sampling (a linear distance of ~2.5 km). Therefore, we extrapolated plot-level climate data by calculating lapse rates based on daily summaries of temperature and total precipitation (rain and snow meltwater equivalent) from NOAA Climate Data for Burlington International Airport (100.6 m elevation, 44.46°N, 73.14°W) and the summit of Mount Mansfield (1,204 m elevation, 44.52°N, 72.81°W), a nearby peak north of Camels Hump in the Green Mountains (data available from https ://www.ncdc.noaa.gov; Figure 1b). Daily mean temperatures were calculated as the mean of the daily minimum and maximum temperature and were subsequently used to calculate mean annual temperature at each station. Likewise, daily precipitation values were summed to estimate annual precipitation (see Appendix S1: Figure S1.1). For each census year, we applied linear extrapolation (i.e. a lapse rate) between the low-and high-elevation climate stations to assign an estimated mean annual temperature and annual precipitation value to each inventory plot, based on its elevation (see Beckage et al., 2008;Siccama, 1974). On average, temperature decreased by 0.53°C and precipitation increased by 9.4 cm for every 100 m increase in elevation. It should be noted that the precipitation lapse rate does not account for precipitation input from cloud deposition, which can be significant at high elevations (Siccama, 1974). Temperature lapse rates estimated by extrapolation closely matched rates based on direct temperature measurements along the transect during 2017 (J. Butnor, personal communication). Across the nine censuses, the mean annual temperature and annual precipitation averaged across the entire elevational gradient on Camels Hump ranged from 2.38 to 4.31°C and from 117 to 203 cm, respectively. We obtained regional atmospheric pollutant nitrogen (combined see Appendix S1: Figure S1.1).

| Generalized dissimilarity modelling
In order to test for spatial and temporal changes in the forest tree community, we fit Generalized Dissimilarity Models (GDMs) using the "gdm" package in R (Manion et al., 2018;R Development Core Team, 2017). GDM is a multivariate form of nonlinear matrix regression that uses monotonic I-splines and generalized linear modelling to model turnover in biological composition (or biological dissimilarity) between sites as a function of environmental and geographic distances (Ferrier et al., 2007). This technique models complex ecological changes in community structure by using flexible splines to individually transform each predictor variable to provide the maximum likelihood fit between pairwise environmental and compositional dissimilarities (Ferrier et al., 2007). Unlike other community modelling approaches, GDM does not assume species respond similarly to environmental and climatic changes.
Response matrices were constructed using pairwise dissimilarity in species abundance (i.e., absolute number of individuals/sample) using the Bray-Curtis index of pairwise dissimilarity (d ij ), defined as the difference in absolute abundance of each species occurring in sample (i) compared to sample (j) divided by the total species abundance between both samples (Bray & Curtis, 1957). We performed two types of GDMs: (a) spatial GDMs that analysed dissimilarity between survey plots (N = 85) across the elevational transect within each census year, and (b) temporal GDMs that analysed dissimilarity of transects across the census years (N = 9). For the spatial GDMs, we used pairwise elevational differences between survey plots to create a single predictor matrix (per year) to analyse shifts in forest tree community zonation along the spatial elevational gradient. For the temporal GDMs, we used pairwise temporal differences across census years in climate (mean annual temperature and total annual precipitation) and atmospheric pollution (pollutants N and S) to create predictor matrices that analyse shifts in community compo- mean annual temperature and annual precipitation per year for each elevation group, and the annual pollutants N and S estimates were used for both elevation groups within a year.
Fitted I-spline turnover functions from GDMs provide two types of information: (a) the shape of the splines indicates the rate of forest tree community change explained by the position along the environmental gradient(s), and (b) the maximum height of the spline indicates the overall magnitude of turnover (i.e., total beta diversity) explained by the gradient (Ferrier et al., 2007). The amount of dissimilarity in species composition between any two points on the environmental gradient is reflected by their corresponding difference along the turnover function, with dissimilarity (i.e., beta diversity) accumulating monotonically. Thus, cumulative beta diversity represents the total amount of dissimilarity explained across the environmental gradient.
Furthermore, turnover functions can have different heights even if the two most distant sampling points are comprised of completely different species. This is attributable to differences in the importance of a given environmental gradient in explaining variation in dissimilarity between pairs of samples.
For each predictor variable, the asymptotic maximum value of the I-spline turnover function was used to quantify the relative contribution of that predictor in explaining forest tree community turnover. We also assessed model and predictor variable significance using matrix permutation (N = 1,000 replicates) and estimated the importance of each predictor by quantifying the decrease in percent deviance explained by a particular variable when it was randomly permuted among the samples (Ferrier et al., 2007). Minimal adequate models were obtained by permuting one variable at a time and reassessing variable importance and significance. At each step, the least important variable was dropped to eliminate non-significant predictors (Manion et al., 2018).

| Species-level responses
To better understand the response of individual species over the sampling period, we compared the density (number of woody stems > 2 cm dbh/m 2 ) of four dominant canopy species, including two broadleaf species common in the low-elevation northern hardwood forest (Ac. saccharum and F. grandifolia) and two coniferous species common in the high-elevation boreal forest (P. rubens and Ab. balsamea) across the elevational gradient for three census periods (1965, 1990 and 2015). Likewise, significance of rank shifts in median elevation and density between these years were tested using Wilcoxon signed-rank tests adjusted within species using Bonferroni corrections to account for the multiple temporal comparisons performed (α = 0.05/3 = 0.02).

| Forest tree community and species turnover along an elevational gradient
Spatial GDMs showed that the single spatial predictor of elevation was sufficient to explain between 53.53% and 63.01% of F I G U R E 2 (a) GDM-fitted I-splines displaying the rate and variation in tree community turnover across the elevational gradient for each census year.
We assessed the uncertainty around the fitted I-spline for the 2015 census with bootstrapping (N = 1,000 iterations).
There was significantly less turnover in the most recent census compared to previous years, as indicated by the error band (± one standard deviation). (b) Inset reflects the asymptotic (maximum) value of each spline to better differentiate the total amount of turnover explained by elevation in each census year the overall deviance in forest tree community composition for each census period (see Appendix S1: Figure S1. Elevational shifts in the range and density of individual species were also evident across the census period (see Appendix S1: Table   S1.1). Among low-elevation species, increases in density and ups-

| Anthropogenic influences on temporal forest tree community change
We used temporal GDMs to elucidate potential anthropogenic drivers of forest tree community change across the 50-year study period. In all temporal GDMs, variation in annual precipitation and pollutant N were not identified as significant predictors during model selection, and were subsequently removed. However, variation in pollutant S and mean annual temperature were both significant predictors in explaining forest tree community turnover through time when considering the whole-mountain (i.e., all plots pooled; Table 1). When analysing temporal GDMs separately for low-and high-elevation plots, pollutant S and mean annual temperature were again significant predictors of composition turnover through time for low-elevation forests ( Figure 5). In contrast, while climate and atmospheric deposition contributed to explaining compositional turnover in high-elevation forests, the overall model and individual predictors were not significant (Table 1). In both the whole-mountain and low-elevation temporal models, pollutant S was found to be the most important driver of turnover (Table 1; Figure 5). Turnover was strongest as temporal differences in pollutant S increased and gradually saturated at the highest pollutant S concentrations (4 mg/L; Figure 5a,e). In contrast, turnover was moderate and showed little sign of saturation in response to temporal differences in temperature over the last 50 years.

| Long-term studies of forest communities along elevational gradients reveal drivers of global change
Species range shifts to higher elevations have been correlated with rising temperatures, both in the past from natural climatic warming (Davis & Shaw, 2001) and more recently with human-induced warming (Parmesan, 2006;Parmesan & Yohe, 2003). The use of elevational gradients is a powerful approach to investigate climate change impacts on biodiversity distribution because climate, and therefore climatically determined species distributions, change dramatically across relatively fine spatial scales (e.g., Brusca et al., 2013;Kelly & Goulden, 2008;Walther et al., 2005). Similar studies have documented elevational shifts in range along fine-scale gradients, but typically do so comparing a pair of discrete census periods (e.g., Damschen, Harrison, & Grace, 2010;Savage & Vellend, 2015;Scherrer, Massy, Meier, Vittoz, & Guisan, 2017).
Our study is unique in providing replicated censuses conducted on the same long-term monitoring plots arrayed across a steep elevational gradient over the course of 50 years. We observed an overall reduction in beta diversity across the elevation gradient in the latest census (2015), reflecting a more homogeneous forest tree community than previously reported. Homogenization of the spatial components of species diversity would be an anticipated consequence of synchronous upslope movement of forest species, as might be predicted from expectations of climate warming alone, but we found little evidence for coordinated upslope shifts in species ranges (Figure 4). Rather, shifts in elevational range and density appeared less cohesive and unique to individual species (e.g., Pucko, Beckage, Perkins, & Keeton, 2011;, with higher densities of common low-elevation species such as F. grandifolia occupying the temperate-boreal ecotone (TBE) while high-elevation species such as P. rubens have expanded downslope, both contributing to increased homogenization and decreased separation of community composition along the elevational gradient.
Our multivariate temporal GDM analysis strongly points to joint effects of sulphate deposition and temperature as drivers of this temporal change in beta diversity, with the greatest changes over time occurring paradoxically at lower elevations ( Figure 5).
Notably, sulphate pollution was found to be the most important driver of temporal forest tree community turnover, which highlights the slow recovery of forests from the accumulated legacy effects of atmospheric pollution (Driscoll et al., 2001).
Previous studies conducted on Camels Hump and throughout the north-eastern region have investigated the effects of climate change on forest communities (e.g., Beckage et al., 2008;Tang & Beckage, 2010), and the contribution of climate and atmospheric pollution on range shifts and growth of individual species (e.g., Kosiba, Schaberg, Rayback, & Hawley, 2017Wason, Dovciak, Beier, & Battles, 2017). Our results corroborate these other studies in finding temporal changes in species densities and distributions, but we are the first to parse out the joint effects of different anthropogenic influences on forest tree community turnover.
The temporal parsing of anthropogenic drivers of changes in beta diversity has led to some novel insights. For example, warming temperatures were previously implicated in models of temporal changes in the elevational position of the ecotone on Camels Hump  and our results further confirm the importance of temperature. Additionally, we identified atmospheric pollution loading as a driver of community change, which is known to have detrimental impacts on below-and above-ground processes in plants. Because Ca plays a critical role in plant carbon regulation and stress response systems, environmental Ca depletion induced by atmospheric pollution can lead to plant Ca deficiencies that predispose trees to decline-often following exposure to a secondary stressor (Schaberg, DeHayes, & Hawley, 2001;Schaberg, Miller, & Eagar, 2010). Decline symptoms associated with acid deposition-associated Ca deficiency have been noted in multiple species including P. rubens (Hawley, Schaberg, Eagar, & Borer, 2006), Ac. Saccharum (Huggett, Schaberg, Hawley, & Eagar, 2007), Halman, Schaberg, Hawley, Hansen, & Fahey, 2015). Likewise, a recent study found Ab. balsamea growth is sensitive to precipitation pH, but the mechanism underlying this sensitivity is unknown . Considering the dominance of these species in montane ecosystems across the region, the change in health and productivity attributed to atmospheric pollution likely played a major role in shifting the ecological community over the past 50 years. Atmospheric pollution has decreased since its regional peak in the mid-1990s and the passing of the Clean Air Act and subsequent amendments (Driscoll et al., 2001;Likens, Driscoll, & Buso, 1996). However, it should be noted that the decline in pollution has been more dramatic and consistent for S deposition (4 to 0.4 mg/L) in comparison with N deposition (1.95 to 1.15 mg/L), which could explain the importance of this atmospheric pollutant in explaining species turnover.
The lack of high-elevation response to changes in atmospheric pollution and climate was unexpected given the increased rates of climate warming, atmospheric deposition, and winter injury that occur at high elevations (Cogbill & Likens, 1974;Lazarus, Schaberg, Hawley, & DeHayes, 2006) and the expected upslope shift in the ecotone . This may reflect the limitation of forest species establishment and survival at high elevations because of high abiotic stress. Likewise, species richness has been shown to peak at mid-elevations (e.g., Kessler, Kluge, Hemp, & Ohlemüller, 2011;Rahbek, 2005;Siccama, 1974) and is generally decreasing towards high elevations. As a result, high-elevation forests have lower forest tree species diversity that could result in inherently less observed turnover through time compared to low-elevation forests. This appears to be a likely explanation for the finding of significant turnover in temporal GDMs for low-but not high-elevation communities, as low-elevation plots shifted in species dominance through time while also receiving new establishment from higher elevation species.
Had our analysis been confined to a single pair of census periods, or if monitoring had ceased after the early 2000s, we would have missed the important effects of both atmospheric pollution and temperature at lower elevations, as well as the downslope shifts in P. rubens and the homogenization of the TBE. Continued monitoring is critical to document further shifts in the forest community as the environment and climate change in a complex, multivariate manner.
However, a caveat of our statistical approach is that we are unable to explicitly account for time-lagged responses to a particular environmental driver. Therefore, the value of our temporal GDMs is on parsing the relative importance of different drivers of compositional turnover as observed across the entire study period rather than between any pair of census years. An additional complication is the potential for spatial and temporal autocorrelation to influence our results. GDM can accommodate this by including geographic (temporal) distance as a covariable, but in our dataset, this created collinearity with environmental predictors. While the shift in community composition could result from natural succession, manipulative experiments have shown the direct effects of atmospheric pollution on forest tree physiology and health , which provides additional evidence that sulphate deposition is a major driver of forest community change.
TA B L E 1 Results from temporal Generalized Dissimilarity Models (GDMs) assessing forest tree community turnover across the entire elevational gradient and separately within low-and high-elevation forests

| Potential drivers of the downslope movement of a high-elevation species
The spatiotemporal changes of P. rubens are not unique to Camels Hump, but are reflective of regional change occurring throughout the north-eastern United States in recent decades (Battles et al., 2003;Foster & D'Amato, 2015;. The decline of P. rubens in the north-eastern United States began in the mid-1960s and was attributed to winter injury induced by atmospheric pollution (Schaberg & DeHayes, 2000). Recent studies of P. rubens, however, document a recovery of the species throughout the region, evidenced by increased growth rates as a result of warming temperatures and reduced atmospheric deposition (Kosiba et al., 2018;. Warming temperatures and the gradual overcoming of physiological impairments following reductions in pollution loading may suggest that P. rubens has been released from its period of decline. Recovery from harvesting during the mid-19th and early 20th centuries (Oosting & Billings, 1951;Whitney, 1988) may also partly explain the downslope movement of P. rubens. While extensive logging did not occur in our study area, harvesting in the surrounding forest may have reduced recruitment through dispersal at lower elevations. Thus, the absence of P. rubens from lower elevation Picea rubens is a late-successional, shade-tolerant species, and as a result, the reproduction and spread of the migration front into low elevations that we document in this study was likely a slow process.
Importantly, this spread into lower elevations is not an expansion of P. rubens into previously unoccupied climates, but rather likely reflects the re-colonization of its historic niche. Records of mature Picea witness trees from pre-settlement surveys (1620-1825) indicate that P. rubens historically occupied lower elevations (<850 m) across the north-eastern United States (Thompson, Carpenter, Cogbill, & Foster, 2013 While our analyses implicated climate and atmospheric pollution as explaining temporal shifts, we cannot rule out many correlated effects, including biotic interactions. Species occurring along elevational gradients often span the lower and upper margins of their physiological tolerances, wherein the lower margin is defined by antagonistic species interactions (e.g., competition and herbivory) and the upper margin is largely regulated by abiotic stress (Bertness & Callaway, 1994;Brooker & Callaghan, 1998). A reduction in lower margin competition due to climatic changes adversely affecting one of the two interacting species could allow for an expansion of the realized niche into unoccupied regions of the fundamental niche for a species with a distribution biotically limited at its lower margin (Lenoir et al., 2010).
Our results demonstrate that P. rubens is able to establish and survive in the warmer and drier conditions found at lower elevations on Camels Hump, and there is ongoing potential for expansion of P. rubens into low elevations, possibly because of a reduction in the abundance of Ac. Saccharum or another low-elevation species experiencing decline (Bishop et al., 2015;. However, while long-term monitoring studies can document important correlative trends, identifying the importance of biotic interactions or other drivers of elevational shifts requires experimental validation (Alexander, Diez, Hart, & Levine, 2016).

| CON CLUS I ON S AND FUTURE DIREC TIONS
A temporal perspective is critical to understanding forest community responses to global environmental change and is particularly important for tree species due to slow migration rates and long life spans. Using long-term monitoring coupled with a multivariate approach to analyse drivers of forest tree community turnover, our results suggest multiple anthropogenic stressors are driving shifts in species density and distribution, which have collectively translated into a more uniform community across the elevational gradient despite the variability of individual species responses. Our work highlights the complexity of species-and community-level responses to ongoing global change and also argues for a better understanding of how anthropogenic influences collectively structure species distributions and how biotic interactions may indirectly and directly be influenced by anthropogenic stressors. The use of single and multispecies common gardens could be used to experimentally assess the impact of altered biotic interactions on plant fitness under different climates, providing a potential explanation for the complex and asynchronous responses of species to environmental changes occurring along steep elevational gradients.

ACK N OWLED G EM ENTS
We appreciate the numerous field crews who helped collect the forest inventory data. Paul Schaberg, John Butnor and Ali Kosiba provided helpful advice and discussion of red spruce physiology and ecology. We

DATA AVA I L A B I L I T Y S TAT E M E N T
All data analysed and presented in the paper are publicly available on Stephen R. Keller https://orcid.org/0000-0001-8887-9213