Site fertility drives temporal turnover of vegetation at high latitudes

Abstract Experimental evidence shows that site fertility is a key modulator underlying plant community changes under climate change. Communities on fertile sites, with species having fast dynamics, have been found to react more strongly to climate change than communities on infertile sites with slow dynamics. However, it is still unclear whether this generally applies to high‐latitude plant communities in natural environments at broad spatial scales. We tested a hypothesis that vegetation of fertile sites experiences greater changes over several decades and thus would be more responsive under contemporary climate change compared to infertile sites that are expected to show more resistance. We resurveyed understorey communities (vascular plants, bryophytes, and lichens) of four infertile and four fertile forest sites along a latitudinal bioclimatic gradient. Sites had remained outside direct human disturbance. We analyzed the magnitude of temporal community turnover, changes in the abundances of plant morphological groups and strategy classes, and changes in species diversity. In agreement with our hypothesis, temporal turnover of communities was consistently greater on fertile sites compared to infertile sites. However, our results suggest that the larger turnover of fertile communities is not primarily related to the direct effects of climatic warming. Furthermore, community changes in both fertile and infertile sites showed remarkable variation in terms of shares of plant functional groups and strategy classes and measures of species diversity. This further emphasizes the essential role of baseline environmental conditions and nonclimatic drivers underlying vegetation changes. Our results show that site fertility is a key determinant of the overall rate of high‐latitude vegetation changes but the composition of plant communities in different ecological contexts is variously impacted by nonclimatic drivers over time.


| INTRODUC TI ON
A pioneering experimental study on British grasslands showed that community-level changes under simulated climate change are dependent on site fertility and associated plant community composition. Communities on fertile sites changed rapidly in response to warming and water manipulations, while minor change was observed in communities on infertile sites (Grime et al., 2000). The explanation provided by the authors was relatively simple: Species exhibiting fast dynamics characterize fertile sites, whereas species with slow dynamics dominate in infertile sites and resisted shortterm experimental treatments. Thus, site fertility and associated plant functional types could determine the rate of plant community changes under climatic change. In the longer term, the communities on infertile sites showed more pronounced changes in composition but it was suggested that these communities would still respond relatively weakly to climatic changes (Fridley, Lynn, Grime, & Askew, 2016;Grime et al., 2008). These findings from British grasslands led to a development of a general view upon the role of site fertility as a key property influencing climate sensitivity of plant communities and might represent one of the few emerging generalizations on how plant communities respond to climate change (Harrison, Damschen, Fernandez-Going, Eskelinen, & Copeland, 2015). However, it remains unknown whether these findings can be extrapolated to a wider range of natural ecosystems along broad geographical gradients. Such knowledge would be imperative for identifying links to experimental manipulations, ecosystem properties, and long-term effects of ongoing climatic warming (Parmesan & Hanley, 2015). The idea behind differential responsiveness of communities on fertile versus infertile sites emerge from plant life histories under different resource levels and evolved plant traits. These underlying foundations of communities have been strongly investigated from 1970s and pioneered by Grime (1974Grime ( , 2001. The broad generalizations of evolved plant strategies are widely applicable in global temperate-boreal ecosystems (Grime, 2001;Grime, Hodgson, & Hunt, 2007). Yet, details in underlying processes and their mechanistic interpretations have been the subject of debate (Craine, 2005). However, it is not known whether plant strategy-based approaches provide any useful predictions for plant community changes under climate warming across broader spatial scales. Despite advances in theory, modeling, and experimentation (Buri et al., 2017;Díaz et al., 2016;Pellissier et al., 2018;Reich, 2014), this remains an unresolved issue, recently been stressed also by Parmesan and Hanley (2015).
It is widely acknowledged that predicting the effects of climate change on plant communities is notoriously difficult due to complex relations between environmental conditions and multispecies assemblages (Fridley et al., 2016). Despite these difficulties, resourcebased plant strategy theories may allow some general predictions of climate change effects. In areas with short growing season, climatic warming can be expected to favor species that can take advantage of the prolonged growing season, which can be accompanied by improved availability of soil resources. Such conditions likely favor species with an ability to gain higher abundance under warmer conditions (Tilman, 1988), for instance, taller species that are superior competitors for light. This could further cause decreases of low-statured species. In the context of Grime's (1974Grime's ( , 2001 plant strategy theory, this could mean an increase of species with the competitive strategy and a decrease of stress-tolerant species.
However, the responses of communities along fertility gradient are likely variable because of differing initial community compositions that result from distinct resource acquisition and conservation traits.
Communities on fertile sites consist to a large extent of plants with fast dynamics (i.e., rapid lifecycles), while communities on infertile sites contain mainly species with slow dynamics (Grime, 1974(Grime, , 2001Reich, 2014). Long-term observational studies at high latitudes have provided tentative evidence that the magnitude of plant community changes under recent climate warming is linked to a gradient from more productive (fertile) to less productive communities (Virtanen et al., 2010), while some broader scale comparative studies failed in finding any clear link between long-term community changes and productivity (Kapfer et al., 2013). It is thus unsettled how fertility modulates the responses of local communities to climate change in different contexts.
In this study, we used forest vegetation resurvey data from several sites distributed over a latitudinal gradient to test whether site fertility is a general ecosystem property that influences the magnitude of plant community responses under long-term climate warming. Specifically, we tested the hypothesis that long-term compositional changes are consistently greater in communities on fertile than on infertile sites. We analyzed temporal turnover of communities and changes in plant strategy types, morphological groups, individual species, and species diversity, and examined how these changes are linked to trends in climatic and nonclimatic drivers.
We searched for historical data on the composition and structure of forest vegetation from 3-5 decades ago at fertile and infertile sites across the boreal zone in Finland. We found several forest vegetation studies from northern Finland that were sampled during 1957-1981 (Table 1) and encompass a south-north gradient from middle boreal to northern boreal zones (see Appendix S1a). More specifically, the surveyed sites represent the following four subregions (with abbreviations): the middle boreal zone (MB), the southern part of the northern boreal zone (NBs), the middle part of the northern boreal zone (NBm), and the northern part of the northern boreal zone (NBn; northern arboreal as defined by Ahti et al. (1968) and Hämet-Ahti (1963)). All the original surveys contained detailed location information (e.g., specific area, description of the surroundings, elevation, aspect, specific vegetation type) allowing sufficiently accurate relocation of vegetation plots (Kopecký & Macek, 2015), which are thus regarded quasi-permanent (sensu Kapfer et al., 2017).
Relocation error and consequent effects of spatial heterogeneity on observed vegetation changes (see Verheyen et al., 2017) were further minimized as the original surveyor helped in relocating fertile sites in subregions MB, NBs, and NBm. For the infertile sites of same subregions, precisely drawn maps on plot locations existed.
In subregion NBn, large plot size increased the possibility to relocate close to original plot location (Table 1). In 2013-2014, we relocated and resurveyed infertile and fertile forest sites, considering the original phenological date. In each original survey, the complete species composition (vascular plants, bryophytes, and lichens) was recorded using a sample plot method and species covers were estimated using percentage cover scales, which were applied in a similar manner during resurveys. The magnitude of human disturbance was estimated visually on site at plot scale and from land use raster maps classified by human activity for the surroundings using a radius of 100 m (further details in Walz & Stein, 2014) in ArcMap software (v.

10.2; ESRI).
Prior to statistical analyses on compositional changes, we explored the spatial distribution of study plots and disturbance levels.
For this study, we included sites that met the following criteria: Each subregion had comparable number of vegetation plots per fertile and infertile sites (Table 1) and plots within each subregion covered about equal spatial extent. Thus, spatially strongly isolated plots (typically more than 100 km apart from the other plots) were left out. The individual plots within each site were situated at distances of 26-2,200 m (the average distance to nearest neighboring plot is 367 m). Moreover, we included only sites without signs of human disturbance (e.g., cuttings, forest management, constructions).
Based on these preset selection criteria, a total of 78 plots were left for analyses (Table 1). Of these, 38 plots represent fertile herb-rich forest communities and 40 plots infertile heath forest communities, which were distributed equally over the latitudinal MB-NBn gradient.
Before analyses, original and resurveyed species data were checked for consistency and some similar species were combined Species data were pooled into morphological plant groups (shrubs, dwarf shrubs, forbs, graminoids, pteridophytes, bryophytes, and lichens), and vascular plants were further categorized into Grime's (1974Grime's ( , 2001 54-57 (1957-1960) Hämet- Ahti (1963) Note: Subregion (MB = middle boreal, NBs = northern boreal/southern part, NBm = northern boreal/middle part, NBn = northern boreal/northern part), forest type (herb-rich = fertile, heath = infertile), soil type (specified in original survey and resurvey), number of resurveyed vegetation plots, plot sizes used in cover estimations for shrub layer species (shrub lr.) and for field and ground layer species (field lr.), average soil pH (based on our own analyses of soil samples collected during resurvey), years from the original sampling and the original data. a Partly paludified, peat-like soil in the resurvey. b Averaged from five 1 × 1 m vegetation grids. c Averaged from two to four 1 × 1 m grids. d Averaged from four 1 × 1 m grids.
stress-tolerant competitors, and generalists) after Pierce et al.
(2017) whose globally calibrated CSR strategies allow comparisons between and within biomes.

| Trends in climatic and nonclimatic drivers in the study area
Daily temperature and precipitation values over the study period were obtained from E-OBS raster dataset (Haylock et al., 2008) for each site. Annual thermal sums were summed from daily temperatures of the growing season that exceeded the +5°C threshold. Growing season was considered to start and to end from the day following at least ten subsequent days with a daily average temperature ≥+5°C to mark the start and <+5°C to mark the end.
Annual precipitation was summed from the daily data. Linear regressions were formed for annual thermal sums and precipita- Reindeer (Rangifer tarandus) is the main herbivore in the study area excluding study sites in the subregion MB that is located outside reindeer herding districts. Data on annual reindeer numbers were obtained from Natural Resources Institute Finland. In the subregion MB, moose densities have been relatively high, with small populations of forest reindeer and roe deer (Siira, Keränen, & Heikkinen, 2009). The fertile sites in subregion MB have likely been grazed by free-ranging cattle long ago, but this activity ended before the first vegetation survey (Happonen, 2016;Kaakinen, 1974). At a European scale, northernmost Finland receives the lowest N deposition (Harmens et al., 2011), but it cannot be ruled out that the changes in atmospheric chemistry and associated nutrient addition or changes in sulfur deposition would have some impacts on vegetation in the longer term (Fowler et al., 2007). As our sampling sites represent forested sites with minimal direct human disturbances, changes in understorey composition and structure can be attributed to successional development toward more closed canopy. Alternatively, canopies may have opened up due to reindeer browsing preventing tree regeneration or insect outbreaks causing tree mortality. To document possible changes in overstoreys between the surveys, data from overstorey structure and cover were obtained from the original studies, and repeated samplings for the same parameters were conducted.

| Statistical analyses
The magnitude of compositional turnover in plant communities over time was quantified by calculating Morisita-Horn similarities between each pair of original and resurveyed plot using the pack- Team, 2018). The Morisita-Horn similarity index was chosen because it is not strongly sensitive to variation in species richness (Chao, Chazdon, Colwell, & Shen, 2006), that is affected to some degree by the varying plot sizes between the study sites (see preliminary analyses in Appendix S2), and to minimize the errors due to observer bias (Kapfer et al., 2017). A Bayesian linear regression model with normally distributed errors was fitted to test the effects of site fertility on the magnitude of compositional turnover using package brms (Bürkner, 2019), an interface to the Bayesian modeling framework Stan (Carpenter et al., 2017). This package was also used to fit the following models. In additional models, linear, nonlinear, and interactive effects of sampling time (e.g., the temporal difference between the original survey and resurvey) were further included as covariates to test if they affect the magnitude of turnover. Nonlinear effects were included as thin-plate splines and their interaction with fertility as a factor smooth interaction. Plot size was not independent of sampling time and thus was left out from the models. Different models were ranked using pareto-smoothed importance sampling leave-  to account for repeated measurements. The plot-specific random effects and residuals were allowed to covary between strategy classes, as a covariance matrix was fitted for both. The proportional covers of strategy classes were logit-transformed before modeling to avoid bias due to boundaries in proportional data and to ensure that model residuals were approximately normal. The model was fit only for main strategy class C and S because of identifiability-the relative cover of R must always be 100% − (C + S).
The results of the model were visualized by calculating average proportional cover of each main strategy class for each time period from the posterior distribution of the multivariate model.
Error variance and plot-specific random effects were omitted from the visualizations.
Changes in species diversity per plot over time were quantified by calculating effective number of species for Simpson diversity measures (Jost, 2006(Jost, , 2007 using the package vegetarian. The Simpson diversity measure gives more weight to more abundant species and thus has the same benefits as the Morisita-Horn index, used for analyzing temporal turnover. The effective number of species is the number of equally abundant species needed to achieve the diversity measure and thus corresponds to a unit conversion of Simpson's diversity index from probability to species (Jost, 2006(Jost, , 2007. Effective species numbers were calculated for vascular plants, Lastly, the direction of beta diversity changes was analyzed, because temporal turnover within each site can be due to either rising compositional similarity (homogenization) or dissimilarity (heterogenization). First, a multiple-sample Morisita-Horn index was calculated for each site in the original survey and resurvey, using the package vegetarian. Resulting value ranges from zero (no turnover between samples) to one (complete turnover between all samples).
The sign of beta diversity change was then calculated as a difference between each pair of original and resurveyed sites (total n = 8). The multiple-sample Morisita-Horn index was used in the calculations to make analyses comparable to those of temporal turnover and alpha diversity change.

| RE SULTS
Thermal sums showed relatively similar increasing long-term trends throughout the study sites (Appendix S4a,b), rising between +182 and +253°Cd over the study period. An exception was the northernmost subregion NBn, where a minor increase (+71°Cd) was recorded.
Annual precipitation increased between +30 and +153 mm but decreased slightly (−10 mm) at infertile site in subregion NBm (Appendix S1a). Reindeer grazing pressure increased over the study period but showed considerable fluctuations especially in the northernmost subregion NBn (Appendix S4c). Estimates of overstorey coverage changed over time in both fertile and infertile sites (Appendix S5).
Substantial increases (MB and NBm) or decreases (NBs and NBn) in canopy cover occurred at fertile sites, whereas more modest decreases took place at infertile sites (canopy cover estimate was not available for the original survey in subregion NBm).
In agreement with our predictions, the magnitude of compositional turnover between original and resurveyed vegetation plots was greater on fertile than on infertile sites across all subregions Because of the diverging trends between plant group diversities, only a weak increase was seen in total diversity (Appendix S10b). Within-site beta diversity was higher in fertile sites during both surveys across the study area (Table 2)
However, our results suggest that the larger turnover of fertile com-  Grime et al. (2000), which is likely due to rather late successional stages of studied communities. Moreover, Grime et al. (2008) found that in the longer term, plant communities F I G U R E 1 Simulated Morisita-Horn turnover from a Bayesian regression model for fertile and infertile sites with 95% credible intervals. The observed means of each subregion during the original survey and resurvey are indicated as gray-scale lines. Turnover is Morisita-Horn dissimilarity index (1-Morisita-Horn index) that was calculated for each pair of original and resurveyed vegetation plot F I G U R E 2 Nonmetric multidimensional scaling (NMDS) ordination grouped by original and resurveyed communities on fertile and infertile sites in subregions middle boreal zone (MB), northern boreal zone/southern part (NBs), northern boreal zone/ middle part (NBm), and northern boreal zone/northern part (NBn). Labels within the confidence ellipses show the locations of group centroids also on infertile sites showed more turnover. Therefore, it is possible that site fertility better predicts the rate of long-term community changes rather than specific long-term compositional outcomes.
Given that all subregions experienced a warming trend and often increase in precipitation, we expected that species with the competitive strategy (e.g., tall herbs, shrubs) increased at the expense of competitively inferior species (e.g., low-growing species) through compensatory competition (Tilman, 1999). Contrary to this expectation, we did not observe a consistent climate-driven increase of com-  (2016) found no temperature-related vegetation changes. A growing amount of evidence is indicating that compositional shifts in forests may lag behind climatic warming due to the microclimatic buffering effect from the canopy layer (De Frenne et al., 2013, 2019 or unchanged snow cover conditions (Kreyling, Haei, & Laudon, 2012).
Thus, stronger compositional shifts and more direct responses to warming may emerge as climate change proceeds.
The lack of general trends in changes of plant morphological and specified strategy groups, individual species and diversity measures strongly suggests that compositional changes in forest understorey communities are linked to nonclimatic drivers. We emphasize that changes in overstorey structure, insect outbreaks, or reindeer grazing likely underlie the observed local changes in communities and interplay with recent climatic changes. Essentially, as has been observed elsewhere (Bernhardt-Römermann et al., 2015;Hédl, Kopecký, & Komárek, 2010;Prach & Kopecký, 2018), these nonclimatic drivers may explain some of the most evident "anomalies" in vegetation changes. In addition, compositional shifts within vegetation types can be highly variable in different geographic contexts in general (see, for example, Maliniemi, Kapfer, Saccone, Skog, & Virtanen, 2018). communities. Thus, the larger turnover on fertile sites may be related to secondary succession due to continuous disturbance (Pierce et al., 2017). An exception to this is the subregion MB that locates outside the reindeer grazing districts. Here, slight increase of competitive species and the release from long-term grazing by cattle may have enhanced the turnover (e.g., Czortek et al., 2018).
It is possible that the impact of nonclimatic drivers discussed above is more pronounced on fertile than on infertile sites, and this could have contributed to the observed differences in the overall magnitude of community changes.
Even though infertile communities showed in general less temporal turnover, rather substantial subregion-specific changes in their species and species-group abundances had occurred. Among infertile communities, the largest turnover had taken place in the sub- Our resurvey analyses revealed that among all studied sites local species diversity had either increased or remained stable over the past few decades. However, we did not observe any pattern in the directional beta diversity changes (homogenization-heterogenization) neither on fertile nor on infertile sites.
The decrease of compositional variation in fertile communities in subregion NBs, which is also associated with a change from formerly herb-rich (Geranium, Filipendula) to recent graminoid-rich (Elymus, Melica) composition, may be indicative of grazing-driven increases of species diversity and compositional homogenization.
In this subregion, grazing increased between 1980 and 2010, and studies elsewhere have shown that increased grazing pressure can contribute to preservation of high local-scale species diversity (Chollet, Baltzinger, Le Saout, & Martin, 2014;Kaarlejärvi, Eskelinen, & Olofsson, 2017). Even so, relatively intense grazing pressure can also lead to more homogeneous community composition (Rooney, 2009). A prominent increase in species diversity was also observed in fertile sites of subregion NBn, likely resulting from grazing and increased light availability as mentioned above, but contrary to fertile sites in subregion NBs, no signs of homogenization were seen. This can be due to high initial species richness and relatively large species pool of the area, buffering the community against homogenization. The increased species richness in reindeer-grazed subareas NBs, NBm, and NBn is related to increasing abundance of local graminoid and bryophyte species. suggests that site fertility can be regarded as a general ecosystem property that predicts the magnitude of compositional turnover in high-latitude plant communities over time.

ACK N OWLED G M ENTS
This research was financially supported by the Academy of Finland (project # 259072). We acknowledge the E-OBS dataset from the EU-FP6 project ENSEMBLES (http://ensem bles-eu.metof fice.com) and the data providers in the ECA&D project (http://www.ecad.eu).
We thank Jouko Kumpula from Natural Resources Institute Finland for providing data on reindeer numbers and Eero Kaakinen, Anu Skog and Karoliina Huusko for assisting in resurveys. We thank for insightful comments given by Kari-Anne Bråthen, Robert Björk and Radim Hédl.
TA B L E 2 Beta diversity, as a community turnover rate with standard error inside brackets, within studied communities in the original survey (β old ) and resurvey (β new ) on fertile and infertile sites in each subregion Note: Rates range from zero (no turnover between sample plots of a site) to one (complete turnover between samples). Difference between old and new turnover value (Δβ) indicates whether community compositions show signs toward heterogenization or homogenization.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
TM and RV initiated the study. KH and TM analyzed the data. TM led the writing of the manuscript. All authors contributed to data collection, designing the analyses and reading and revising the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data used in this study will be deposited in BioTIME (Global database of biodiversity time series, biotime.st-andrews.ac.uk) and is available upon request from the corresponding author (tuija.malin-iemi@oulu.fi).