Trait-based responses to forestry and animal husbandry modify long-term changes in forest understories

Aim Land use is the foremost cause of global biodiversity decline, but species do not respond equally to land-use practices. Instead, responses are mediated by species traits, but long-term data on the trait-mediated effects of land-use on communities is scarce. Here we study how forest understory communities have been affected by common land-use practices during 4–5 decades, and whether changes in diversity are related to changes in functional composition. Location Finland Time period 1968–2019 Major taxa studied Vascular plants Methods We resurveyed 245 semi-permanent vegetation plots in boreal herb-rich forest understories, and used path analysis to relate changes in diversity, species composition, average plant size, and leaf economic traits to reindeer husbandry, forest management intensity, and changes in canopy cover and canopy traits. Results Forestry affected understories indirectly by increasing canopy shading, which increased understory SLA and decreased LDMC over the study period. Intensive management also decreased species richness and increased turnover. In areas with reindeer husbandry, reindeer density had increased along with understory evenness and diversity. Areas with reindeer husbandry also had lower temporal community turnover. Plant height had increased in areas without reindeer, but this trend was suppressed or even reversed within the reindeer herding area. Main conclusions Functional traits are useful in connecting vegetation changes to the mechanisms that drive them. Forest management causes directional selection on light-interception traits by altering shade. Reindeer husbandry seems to buffer forest understory communities against compositional changes by altering selection on whole-plant traits such as size. These trait-dependent selection effects could inform which species benefit and which suffer from different types of land use, and point to the potential usefulness of large herbivores as tools for managing vegetation changes under global change.


Introduction
Decreasing habitat quantity and quality are the lead drivers of biodiversity change in terrestrial ecosystems, caused primarily by human land use and land-use change (Díaz et al., 2019) . However, the most common metrics used to monitor biodiversity change do not necessarily capture the processes leading to biodiversity change. Although species extinction rates are estimated to be high above background levels globally (Barnosky et al., 2011) , several studies have shown that regional-to-local-scale diversity can actually remain stable or even increase with time (Blowes et al., 2019;Sax & Gaines, 2003;Steinbauer et al., 2018;Vellend et al., 2013) . Compositional change has often been studied using temporal dissimilarity metrics, which are informative about the rate of compositional turnover (Blowes et al., 2019) . However, dissimilarity metrics are directionless, and are at best indirectly informative about changing selection pressures on species affected by anthropogenic stressors such as land-use change.
One way to capture this information more directly is to monitor changes in the functional composition of communities. Unlike dissimilarity metrics, changes in average functional traits are informative of changing selection pressures on vegetation (Violle et al., 2007) , and resulting ecosystem consequences (Lavorel & Garnier, 2002) . Indeed, functional traits have been suggested as one set of essential biodiversity variables needed for monitoring biodiversity change (Pereira et al., 2013) , yet so far they have been applied in only a few studies to demonstrate the importance of land-use legacies on vegetation changes in forests (e.g. Hedwall et al., 2019;Maes et al., 2020;Perring et al., 2018;Verheyen, Honnay, Motzkin, Hermy, & Foster, 2003) .
The above-ground traits of vascular plant species and communities vary on two independent axes that are important for different kinds of environmental constraints. Plant size determines competitive hierarchies in relation to light, while the leaf economics spectrum describes trait covariation related to e.g. survival-growth trade-offs (Bruelheide et al., 2018;Díaz et al., 2016) . For example, increased light and nutrient availability due to recent forest cuttings should lead to understory communities composed of higher plants, since only taller species can reach the light in conditions of high resource availability (Blondeel et al., 2020) .
In this study we apply trait-based approaches to study vegetation changes in the biodiversity hot-spots of Finnish boreal forests. Boreal forests constitute almost a third of global forest cover, (Brandt, Flannigan, Maynard, Thompson, & Volney, 2013) , and most of these forests are managed for industrial wood production (Gauthier, Bernier, Kuuluvainen, Shvidenko, & Schepaschenko, 2015) .
There have been relatively few studies about the long-term effects of forest management on the structure and diversity of boreal forest understories (see e.g. Hedwall et al., 2019) . Boreal forest understory species have strong species-specific responses to tree cuttings that are related to shade tolerance (Tonteri et al., 2016) . Consequently, one would expect the average leaf economic traits in communities to change in response to altered light levels in managed forests. In Fennoscandian boreal forests, forest management has intensified manyfold from low intensity wood harvesting to industrial-scale forestry after the 1950s (Löfman & Kouki, 2001) , with detrimental effects on forest biodiversity (Hyvärinen, Juslén, Kemppainen, Uddström, & Liukko, 2019) . More information is needed on how different forest management practices affect understories in the long-term (Su, Wang, Huang, Fu, & Chen, 2019) , especially since management effects on community composition can sometimes lag decades behind changes in management practices (Muurinen, Oksanen, Vanha-Majamaa, & Virtanen, 2019) .
In addition to forestry, humans also affect understories by controlling herbivore densities, for example by practicing animal husbandry. Since herbivores have strong effects on understory vegetation (Bråthen & Oksanen, 2001;Oldén, Raatikainen, Tervonen, & Halme, 2016) , human-mediated changes in herbivore density can be regarded as land-use change. According to theory, the densities of large herbivores should be reflected in the average community traits that control susceptibility to grazing, such as plant size and leaf palatability (Díaz et al., 2007) , but interpreting grazing effects in the light of quantitative functional traits is still rare outside grasslands.
One of the main large herbivore species in boreal forests is the reindeer ( Rangifer tarandus ), which in many parts of Eurasia is semi-domesticated free-ranging livestock. Several local scale studies have shown reindeer herbivory to cause decreases in biomass and declines in favoured forage species such as herbs, deciduous shrubs, and lichens in the boreal coniferous heath forests and on tundra heaths (Bråthen & Oksanen, 2001;Olofsson, Moen, & Östlund, 2010;Sundqvist et al., 2019;Väre, Ohtonen, & Oksanen, 1995) . Connecting these effects to plant functional traits could shed light on whether these effects scale to larger geographical extents, and how changes in reindeer densities are reflected in long-term vegetation dynamics.
Here we study the effects of land-use on long-term vegetation changes in boreal forest understories by employing diversity-, dissimilarity-, and functional trait -based approaches. Specifically, we ask three questions: 1) how have plot-scale species diversity, species composition, and functional composition changed in boreal herb-rich forests over 40-50 years, 2) how are the changes correlated with changes in forest management intensity and semi-domesticated reindeer grazing pressure, and 3) are changes in functional traits correlated with changes in diversity across both time and space, hinting at the existence of general processes that link traits to diversity?
We answer these questions with a vegetation resurvey study along gradients of management intensity and grazing pressure intensification in boreal herb-rich forests, using Bayesian regression modelling and path analysis.

Vegetation data and resampling
We conducted vegetation resurvey in boreal herb-rich forests in Northern Finland (Fig. 1). The resurveyed herb-rich forest sites were originally surveyed in 1968-75 by Eero Kaakinen (Kaakinen, 1971(Kaakinen, , 1974Unpublished) . The purpose of the original surveys was to identify phytosociological species associations in mature herb-rich forests. For this reason, the surveys were conducted in relatively undisturbed habitats. These forests harbour higher species richness than the surrounding less fertile and homogenous boreal forests. They occur typically as small patches and cover only a fraction of the forested area in northern Scandinavia, resulting primarily from the scattered distribution of calcareous bedrock and soil. The field layer of these forests is dominated by herbs such as Geranium sylvaticum and Filipendula ulmaria , ferns such as Gymnocarpium dryopteris and Athyrium filix-femina , and graminoids such as Milium effusum and Elymus caninus . Occasionally they also host species typically found in less fertile boreal forests, such as the dwarf-shrubs Vaccinium myrtillus and Vaccinium vitis-idaea . The three most abundant tree species Picea abies , Alnus incana , and Betula pubescens formed more than 80% of tree cover in our study plots, both during the original and resampling time periods.
Using detailed notes on the locations of the original sites, which in many cases included additional information on topographical position (elevation, aspect), we relocated and resurveyed the majority of the original 336 sites in 2013-14 and in 2019. Most surveys were done in teams of ≥2 investigators.
We omitted sites with vague descriptions of location, and some sites that were located very far from other sites. The final number of plots, after applying exclusion criteria described below, was 245.
Relocation error is inherent in vegetation resurveys, but was minimized as much as possible by the relocation help of the original surveyor, the patchy distribution of the studied vegetation type, and the fact that herb-rich forests are easily told apart from the boreal forest matrix (Kapfer, Hédl, & Jurasinski, 2017) . This study thus applies all of the measures listed by Verheyen et al. (2018) for minimizing observer and relocation error. During both the original survey and the resurvey, species composition was estimated from 5 m x 5 m quadrats for the field layer and from 10 m x 10 m quadrats for shrub and canopy layers. The smaller quadrats were nested in the larger quadrats. The abundances of field and shrub-layer species were estimated as absolute percentage covers, and the covers of tree-layer species were estimated as relative covers. In addition, total canopy cover was estimated visually on a three-point ordinal scale (0-30%, 31-70%, 71-100%). In this study, we focus only on the vascular plants of the field and tree layers (for species list and trait values, see Table S1.1 in Appendix S1 in Supporting Information).

Environmental changes
We estimated forest management intensity during the resurvey based on the protocol used for the assessment of threatened habitat types in Finland (Kouki et al., 2018, p. 180) . The original protocol assigns sites on an ordinal scale of habitat quality ranging from zero to four. Here, we inverted the scale to describe management intensity. We also omitted the few sites that were destroyed by deforestation and land conversion to arrive at a management impact variable ranging from one to four.
The management intensity assessment criteria include the structure of the canopy and vegetation layers, disturbance of the soil, and the amount of deadwood, and are described in Table S2.2 (See Appendix S2). These criteria are applied in a one-out-all-out manner, meaning that if any of the criteria for a more intensive management class were fulfilled, the site was assigned in that class even if it would have been a less intensively managed site based on all other criteria. The criteria are independent of the other variables used in this study, meaning that species composition was not used as a measure of management intensity. Signs of forestry were present in ca. 70% of the resurveyed plots, whereas the rest ca. 30% were deemed unmanaged ( Fig. S2.1).
Reindeer are the major large grazer in our study area. In Finland, semi-domesticated reindeer occur only inside a designated reindeer herding area. Data on reindeer density trends in the different reindeer herding districts were provided by Jouko Kumpula from the Natural Resource Institute of Finland. Reindeer densities in the herding area had increased ca. 40% during the study period ( Fig.   S2.1). The population sizes of other large herbivores in the study area, as inferred from hunting statistics, are negligible (Natural Resources Institute Finland, 2020, Fig. S2.1) , with the exception of the moose ( Alces alces ). Although the moose population in Finland has grown due to changes in forestry practices and hunting, the growth has been rather uniform across the country (Nygrén, 2009) .
Moreover, herbaceous field layer vegetation constitutes less than 5% of moose summer diets (Hjeljord, Hovik, & Pedersen, 1990) . Thus, while moose density changes might have affected forest regeneration and shrub layer composition, the direct effects of increased moose densities should be rather small and uniform across the study area.
Because of our focus on land-use change, we intentionally left climate out of our models. On average, increases in minimum and maximum summer temperatures in our study plots have been 0.6℃ and 0.8℃, respectively ( on field-layer vegetation, these will not be ignored but will manifest as unexplained trends in our models.

Community properties
All data processing and analyses were done in R (version 3.6.2 , R Core Team, 2019) .
As species diversity measures, we calculated species richness, Shannon diversity, and species evenness. By Shannon diversity we mean the so-called true diversity of order 1, which is a unit conversion of the Shannon diversity index from entropy to species (Jost, 2006) . As species evenness, we used Pielou's J, or the ratio of Shannon entropy to log-transformed species richness, which is a measure of relative evenness ranging from zero to one (Jost, 2010) .
We calculated community weighted means (CWM) of log-transformed trait values for one size-structural trait and two leaf economic traits: vegetative height, specific leaf area (SLA), and leaf dry-matter content (LDMC). Trait values were log-transformed because generally, traits follow a log-normal distribution (Bjorkman, Myers-Smith, Elmendorf, Normand, Thomas, et al., 2018) .
Species relative covers were used as weights. Trait values were retrieved from the TRY database (Kattge et al., 2020(Kattge et al., , version 5, 2011 and the LEDA database (Kleyer et al., 2008) , and supplemented with our own measurements for common species that were not found from the databases. In four plots, some of the traits were available for less than 80% of total cover. These plots were excluded from analyses. In the remaining plots, trait values were available for >99.5% of total cover for all traits.
We also calculated a CWM for the SLA of the canopy layer. SLA was not available for the non-native We subtracted the diversity measures and trait CWMs of the original sampling period from those of the resamples to get estimates of community change between the time-periods.
We used the R package vegan (version 2.5-6, Oksanen et al., 2019) to calculate temporal turnover in community composition using the version of Jaccard distance that takes into account species abundances. Plot-scale turnover was logit-transformed before analyses. We also calculated temporal turnover for the entire dataset.

Modelling
Unless otherwise specified, models were built using the R package for Bayesian modelling brms We used those sites for which we had the full set of response and predictor variables: changes in species richness, species diversity, species evenness, height, LDMC, SLA, canopy cover, and canopy SLA, location within the reindeer herding area, and management intensity. The final dataset thus consisted of 245 resampled sites.
We used two multivariate linear regression models to analyze plot-level changes in species diversity and composition. The first model tested average changes across the study area. It was an intercept-only model, with observed plot-level differences in species richness, species diversity, species evenness and functional composition, and turnover as the responses. We also modelled residual correlations among the variables. The second model tested the effects of canopy changes and increased herbivore density on the same seven community change metrics. In this second model, we added as predictors the observed changes in canopy cover and canopy SLA, and whether the plot was in the reindeer herding area. Since forest management can consist of varying practices e.g. ditching and other practices that have effects on soils, the residual bias of this second model was further studied in the different management intensity classes to see if management intensity had additional effects on community characteristics besides those mediated by changes in the canopy layer.
We built two additional models to study the effects of forest management on canopy characteristics.
First, we modelled tree-layer canopy cover in the plots as a function of plot identity, sampling time, and management intensity. Since canopy cover was measured on a three-level ordinal scale, we used an ordinal regression model with a cumulative logit link (Bürkner & Vuorre, 2019) . Plot identity was added as a group-level (random effect), and management intensity, sampling time, and their interaction were added as population-level (fixed) effects. Second, we modelled canopy SLA with the same predictors, but using ordinary linear regression with normally distributed errors.
We used the canopy models, i.e. the established relationship between management intensity and canopy SLA and height, to predict the indirect effects of management intensity on community diversity and composition via changes in canopy characteristics. We made 4000 draws from the posterior distributions of the model parameters, and used them to acquire fitted values of canopy cover and canopy SLA in each management intensity class during original sampling and resampling. For the ordinal model, this was the probability of each of the canopy cover classes. We then used these fitted values to calculate the posterior probability of canopy cover transitions (canopy cover decreased, canopy cover stayed the same, canopy cover increased), and the posterior probability distribution of canopy SLA changes. We then used these posterior distributions of canopy change as predictors in the model of community change described above to infer the indirect effects of forest management intensity on community properties via effects on canopy characteristics.
To assess the generality of the covariation between traits and diversity metrics, we inspected whether temporal correlations resembled the spatial covariance between community properties. If the temporal covariation between community properties is driven by mechanistic linkages, the same correlation structure should be evident spatially as well. To quantify spatial correlations in community properties, we built a model explaining community diversity and functional composition with plot as a group-level (random) effect, time as a population-level (fixed) effect, and modelled residual correlations in the response variables. We further investigated the relationship between vegetation height and species richness, species diversity, and species evenness at different times using single predictor generalized additive models with the R package mgcv (version 1.8-31 Wood, 2011) .
In all models, priors were left as the brms defaults. All population-level effects were given a noninformative flat prior over the reals. The standard deviations of all group-level effects were given a half-t prior with 3 degrees of freedom, each prior was further scaled with a constant that equaled 10 standard deviations of the response. Correlations among responses were given a noninformative LKJ prior. There were no divergent transitions in any of the MCMC chains in any of our models, which would have indicated unreliable estimates of the parameter posterior distributions (Carpenter et al., 2017) . In addition, all parameters converged to a stationary posterior distribution (R-hat < 1.01, Gelman & Rubin, 1992) .

Data availability
The data and scripts used to produce these results have been deposited to Zenodo (Happonen et al., 2020) .

Results
Across the study area, we detected no changes in species richness (Fig. 2 a), but observed increases in plot-level Shannon diversity (Fig. 2 b) and evenness (Fig. 2 c) during the sampling interval. There was no detectable change in average species height or SLA (Fig. 2 d,f), but average LDMC had increased, indicating slower leaf economics (Fig. 2 e). Average turnover (Fig. 2 g) between the sampling periods was about 80%. For the full dataset (i.e. the relative abundances of all species across the study area), temporal Jaccard dissimilarity was 37%, which, when divided over the average sampling interval (44 years), corresponds to a yearly turnover of 2.2%.
Management intensity modulated changes in the canopy. In the three lowest management intensity classes, canopy SLA decreased, while in the highest management intensity class canopy SLA remained unchanged (Fig. 3 a). SLA of the tree layer was highly negatively correlated with the relative cover of Picea abies (Fig. S2.4), the main evergreen tree species in these forests. In all the management intensity classes, increasing canopy cover was always more probable than decreasing canopy cover (Fig. 3 b). In the highest and lowest management intensity classes, lack of canopy cover class transitions was the most probable outcome, while in the two middle classes increases in canopy cover were much more probable than other outcomes (Fig. 3 b).
Neither canopy changes nor location within the reindeer herding area had high probability of affecting species richness changes (Fig. 4 a). However, increases in Shannon diversity and species evenness were both positively associated with the presence of reindeer (Fig. 4 b,c). Plots within the reindeer herding area were associated with a more negative trend in plant height (Fig. 4 d). Increased canopy cover was associated with increased SLA and decreased LDMC, while canopy SLA changes were positively correlated with LDMC changes and negatively correlated with changes in SLA (Fig. 4 d,e).
Community turnover was higher in plots outside the reindeer herding area, and quite probably also in plots where canopy SLA had decreased (Fig. 4 g).
In addition to these explainable trends, we detected increasing trends in height and LDMC that were unassociated with our selected predictors (base change in Fig. 4 d & e).
Examining the bias in the residuals of the model explaining changes in community properties with canopy changes and the reindeer reindeer herding area revealed that management intensity seemed to have additional direct effects on some variables in addition to effects modulated by canopy changes (Fig. 5). Specifically, higher management intensity was associated with decreasing species richness ( Fig. 5 a), increasing height (Fig. 5 d), decreasing SLA (in the highest management intensity class, Fig. 5 f), and increased turnover (Fig. 5 g).
Simulating the indirect effects of management intensity on field layer diversity and composition via changes in canopy characteristics revealed that the indirect effects were the strongest in the two middle classes of management intensity (Fig. S2.5). This was due to canopy cover changes being the most pronounced in these classes.
Spatial correlations and temporal covariance of community properties were highly similar (Fig. 6).
Diversity and evenness measures correlated positively with each other, but negatively with height.
LDMC had a positive correlation with height, but a negative correlation with SLA.
Vegetation height was unimodally related to all diversity measures across both sampling times ( Fig.   S2.6).
The direct and indirect effects of land use and land-use change identified in this study are summarized in Fig. 7.

Discussion
In our study, most observed changes in functional composition were directly linked to land-use in ways that are consistent with theory, highlighting the suitability of these metrics for biodiversity monitoring in the anthropocene (Pereira et al., 2013) and for conservation planning (Keddy, 1992) . As examples of trait-modulated effects, forest management increased canopy shading in the long-term, followed by shifts in leaf traits related to resource-use strategy and light interception, while increased grazer density decreased average plant size and, consequently, increased species evenness and species diversity.
Our analysis revealed that many community characteristics were stable over the 4-5 decades study interval across the study area, but displayed trends when conditioned on land use. Apparent stability can thus be the sum of opposing trends across different land-use types, a type of ecological Simpson's paradox. Reading too much into these raw changes is thus discouraged. However, for the sake of comparison, we report that plot-scale species richness was very stable as in a previous meta-analysis (Vellend et al., 2013) , even though plot-scale species turnover was very high ( ca. 80%). Our results are thus consistent with the hypothesis that local species richness is highly constrained by regional processes such as compensatory colonization-extinction dynamics, meaning that lost species are replaced by species better suited for new conditions (Supp & Ernest, 2014) . Furthermore, despite high plot-scale turnover attributable to high species richness of these biodiversity hotspots and inherent relocation error in resurvey studies, temporal turnover for the full dataset was of the same magnitude as average turnover in terrestrial ecosystems in a recent meta-analysis (2% year −1 , Blowes et al.,

2019) .
Our results strongly suggest that increasing reindeer density in areas with reindeer husbandry is affecting the change-trajectory of field-layer vegetation in these diverse northern forests by selectively retarding the relative abundance of large plants with high SLA. Strong effects of the reindeer herding area on observed field-layer changes are consistent with reported experimental effects of grazing on plant communities in Scandinavian mountain forests (Sundqvist et al., 2019) . In the absence of canopy changes, plots within the reindeer herding area had decreasing or stable plant size, decreased SLA, and increased Shannon diversity compared to increasing plant size, stable SLA and stable diversity outside the area. SLA has been shown to correlate positively with nutrient concentrations (i.e. forage quality) in leaves (Díaz et al., 2016) , while large plants have been shown to be more susceptible to mammalian herbivory (Díaz et al., 2007) . Our data thus support the suggestion that large size and high SLA could make species especially vulnerable to herbivory.
Communities in the reindeer herding area also had lower turnover than communities outside it. Since the reindeer herding area is generally in higher latitudes than areas outside it, it is possible that the observed effect is due to forests in colder climates having slower growth, smaller species pool, and thus slower community dynamics. Such differences in turnover rates between fertile and less diverse infertile sites have previously been reported for boreal forests (Maliniemi, Happonen, & Virtanen, 2019) . A further explanation could be that changes induced by increased reindeer densities are counteracting changes unrelated to them. We found this to be the case with plant size -vegetation height increased outside the reindeer herding area, but was stable or decreased within it. These types of effects have previously been reported in the tundra, where reindeer and other large herbivores have been found to inhibit shrubification, tree-line creep, and associated ecosystem changes (Olofsson & Post, 2018) . Our study thus suggests that the buffering effect of large herbivores against vegetation changes extends from the tundra to forest understories as well, because they select against large species favored by climate change.
Forest management had direct and indirect, canopy-mediated effects on understory properties.
Directly, turnover was greater in more intensively managed forests, agreeing with previous studies from other biomes that have reported accelerated biodiversity change following human disturbance (Barlow et al., 2016;Lake et al., 2000) . Management intensity also had clear negative effects on plot-scale species richness, which could be the result of concentrating resources to the tree-layer by practices such as ditching, which increase the growth of the tree-layer but decrease resource availability in the field layer. In addition, SLA had decreased in the highest management intensity class which included recent clearcuts, highlighting the fact that SLA is a key trait in light-availability -based species sorting (Dahlgren, Eriksson, Bolmgren, Strindell, & Ehrlén, 2006) . Indirectly, forest management affected light-interception traits in the understory. These effects were most strongly mediated by increased rate of canopy closure in the two middle classes of management intensity, which include mature human-impacted stands but not recent clearcuts. Increased canopy cover had led to increased SLA and decreased LDMC in the understory. SLA is a succession trait, which increases in response to shading (Dahlgren et al., 2006) , whereas LDMC has probably decreased because of its negative covariance with SLA as part of the community-level leaf economics spectrum (Bruelheide et al., 2018) . Increased light availability following harvesting, along with its associated vegetation changes (Tonteri et al., 2016) , can thus be transient and reversed in the following decades.
In addition, plots had increased field-layer LDMC in the absence of canopy changes, indicating slower leaf economics. This unexplained change could be an effect of some global change driver we did not include in the models. Since warmer climates are generally associated with faster leaf economics (Bjorkman, Myers-Smith, Elmendorf, Normand, Rüger, et al., 2018;Wieczynski et al., 2019) , which would correspond with lower LDMC, climate change is not a likely driver for this trend.
Further research is thus needed on the causes of this trend.
Temporal and spatial covariation of community properties was extremely similar, hinting that the diversity and functional composition of these forest understory communities are not independent of each other, but are controlled by some mechanistic constraints. Unsurprisingly, leaf traits seemed to covary according to the leaf economics spectrum. However, many of the strongest correlations involved average height; plots with higher vegetation had lower Shannon diversity and species evenness, and higher LDMC. The fact that richness, diversity and evenness were all unimodally related to vegetation height across both sampling times confirmed that the diversity-biomass relationship is at least partially behind this pattern. Communities composed of larger species often have fewer plant individuals and exhibit decreased evenness and richness due to the dominance of large plants (Oksanen, 1996) . Reindeer grazing effects are most likely a result of this coordination: tall species are being eaten (Kaarlejärvi, Eskelinen, & Olofsson, 2017;Sundqvist et al., 2019) and dominance thus reduced, which increases diversity (Grime, 1973) . In this study we did not take into account intraspecific variation in plant height, and since grazing obviously affects vegetation via both compositional and intraspecific effects, the true effects of grazing on vegetation height and diversity are probably even stronger.

Conclusions
Our study demonstrates that functional trait-based approaches are highly useful in connecting vegetation changes to land use in a globally applicable way. In the long-term, forest management densens canopies and decreases light availability in the understories of boreal herb-rich forests. This leads to field-layer communities with leaf traits that maximize intercepted light per invested biomass, such as high SLA. In addition, forest management decreases plot-scale diversity and increases compositional turnover. In the absence of herbivory or canopy changes, we found that field-layer vegetation had grown in height. However, this trend was reverted or even reversed within the reindeer herding area, where the human-controlled density of grazing ungulates had increased during the sampling interval. Reindeer husbandry thus seems to reduce plot-scale turnover and dominance, and increase diversity, inviting further study into the use of large herbivores to manage vegetation changes in northern ecosystems.  , species evenness (E, c), functional composition (d-f), and species composition (g), across the study area. The points are medians, the narrow and thick lines are 95% and 66% credible intervals, respectively, and the gray shading is a kernel density plot of the parameter's posterior distribution. LDMC = leaf dry matter content, SLA = specific leaf area.  The effects of changes in canopy cover and canopy SLA, and location within the reindeer herding area on changes in species richness (D0, a), Shannon diversity (D1, b), species evenness (E, c), functional composition (d-f), and species composition (g). The points are medians, the narrow and thick lines are 95% and 66% credible intervals, respectively, and the gray shading is a kernel density plot of the parameter's posterior distribution. 'Base change' corresponds to the intercept of the model, and depicts changes in the absence of canopy changes and reindeer. 'Canopy SLA increase' depicts changes resulting from increasing log canopy SLA by 1.

Fig. 5:
The posterior distribution of residual bias in different management intensity classes from the model explaining changes in community properties with canopy changes and reindeer. The narrow and thick lines are 95% and 66% credible intervals, respectively, and the gray shading is a kernel density plot of the residual bias' posterior distribution. D0: species richness, D1: Shannon diversity, E: species evenness, LDMC: leaf dry matter content, SLA: specific leaf area.  Temporal changes in the diversity and composition of boreal herb-rich forest understory communities, and the long-term direct and indirect contributions of forest management and reindeer husbandry to those changes. This figure is a simplified summary of our key findings. Ambient trends are trends that are not attributable to the explanatory variables used in this study.