Trait-based responses to forestry and reindeer husbandry modify long-term changes in forest understories Running title: Trait-based responses to forestry and reindeer husbandry

Aim:​ Land use is the foremost cause of global biodiversity decline, but species do not respond equally to land-use practices. Instead, it is suggested that responses vary with 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 plant diversity are related to changes in functional composition. Location:​ Finland Time period:​ 1968–2019 Major taxa studied:​ Vascular plants Methods:​ We resurveyed 245 vegetation plots in boreal herb-rich forest understories, and used hierarchical Bayesian linear models to relate changes in diversity, species composition, average plant size, and leaf economic traits to reindeer abundance, forest management intensity, and changes in climate, canopy cover and composition. Results:​ Forestry decreased species richness and increased turnover, but did not affect functional composition. Increased reindeer densities corresponded with decreased height and specific leaf area, and increased leaf dry matter content, evenness and diversity. Successional changes in the canopy were associated with increased specific leaf area and decreased leaf dry matter content and height over the study period. Effects of reindeer abundance and canopy density on diversity were partially mediated by vegetation height, which had a negative relationship with evenness. Climate change had no discernible effect on any variable. Main conclusions:​ Functional traits are useful in connecting vegetation changes to the mechanisms that drive them, and provide unique information compared to turnover and diversity metrics. These trait-dependent selection effects could inform whether species benefit or suffer from land use changes and explain observed vegetation shifts 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 (i.e. species richness) do not necessarily capture its full extent (Hillebrand et al., 2018) . For example, although species extinction rates are estimated to be high above background levels globally (Barnosky et al., 2011) , several studies have shown that regional and local-scale diversity can actually remain stable or even increase with time (Blowes et al., 2019;Steinbauer et al., 2018;Vellend et al., 2013) . Compositional turnover has often been studied using temporal dissimilarity metrics, which are informative about the rate change (Blowes et al., 2019) . However, as dissimilarity metrics are directionless, they only indirectly inform us about selection pressures that land-use change exerts on communities. To better connect observed biodiversity changes to the processes that drive them, more than diversity and dissimilarity metrics are needed.
Functional traits have been suggested as one set of essential biodiversity variables needed for monitoring biodiversity change (Pereira et al., 2013) . However, monitoring is only effective if we know the potential causes of changes in the monitored attributes, otherwise observed changes cannot inform effective interventions. One way to move beyond describing patterns of diversity change towards understanding the mechanisms behind them is to link changes in the functional composition of communities to potential selective forces, such as different types of land-use change. Changes in average functional traits can be directly informative of changing selection pressures on vegetation, and resulting ecosystem consequences (Lavorel & Garnier, 2002;Violle et al., 2007) . However, traits have so far been applied in only a few studies to demonstrate the importance of land-use legacies on vegetation changes in forests (Depauw et al., 2020;e.g. Hedwall et al., 2019;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 can provide valuable clues to identifying important variables affecting changes in biodiversity. 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 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) . Altered average vegetation height may in turn alter diversity and evenness due to the well-known unimodal relationship between diversity and productivity (Grime, 1973) .
In this study we apply a trait-based approach to study vegetation changes in the biodiversity hot-spots of boreal forests in Finland. Boreal forests constitute almost a third of global forest cover and most of these forests are managed for industrial wood production (Gauthier, Bernier, Kuuluvainen, Shvidenko, & Schepaschenko, 2015) . In Fennoscandian boreal forests, forest management has intensified manyfold from low intensity wood harvesting to industrial-scale forestry after the 1950s 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 those average community traits that control susceptibility to grazing, such as plant size and leaf palatability (Díaz et al., 2007) .
One of the main large herbivore species in boreal forests is the reindeer/caribou ( Rangifer tarandus ), which in many parts of Eurasia is semi-domesticated free-ranging livestock. Several local scale studies have shown reindeer herbivory to cause declines in biomass and in favoured forage species in the boreal coniferous heath forests and on tundra heaths (Bråthen & Oksanen, 2001;Olofsson, Moen, & Östlund, 2010;Sundqvist et al., 2019) . Emerging evidence suggests that selection by reindeer is related to plant functional traits (Kaarlejärvi, Eskelinen, & Olofsson, 2017) , which suggests that large-scale variations in reindeer densities could be an important driver of long-term vegetation dynamics.
Here we study the effects land-use on long-term vegetation changes in boreal forest understories by employing diversity-, dissimilarity-, and functional trait -based approaches. Specifically, we ask four questions: 1) have plot-scale species diversity, species composition, and functional composition changed in boreal herb-rich forests over 40-50 years, 2) are there spatial differences in these changes 3) are the changes related to forest management intensity or semi-domesticated reindeer abundance, and 4) are changes in functional composition 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 hierarchical Bayesian regression modelling.

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 (Maliniemi, Happonen, & Virtanen, 2019) , are an important habitat for many threatened species (Kouki et al., 2018) , and are thus important to monitor. They occur typically as small patches and cover only a fraction of the forested area in northern Scandinavia reflecting the scattered distribution of calcareous bedrock and soil. The field layer 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 typical for 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, during both the original survey and the resurvey.
Using detailed notes on the locations of the original sites, 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, which lacked detailed enough 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 (Kapfer et al., 2016) , 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. 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 and diversity were not used as measures of management intensity. Signs of forestry were present in ca. 70% of the resurveyed plots, whereas the rest ca.
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 abundances in the different reindeer herding districts were provided by Natural Resource Institute Finland. Reindeer densities in the herding area had increased ca. 40% during the study period, from 1 to 1.4 km⁻² ( Fig. S2.1). We 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

Community properties
All data processing and analyses were done in R (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 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. References to the original trait datasets are listed in Appendix S3.
We also calculated a CWM for the SLA of the canopy layer. SLA was not available for the non-native Larix sibirica (present on 1 site), for which we used the SLA of Larix decidua instead.
Summary statistics of the studied properties during original sampling ca. 1970 are presented in Tables S2.3 & S2.4, which help in interpreting the observed rates of change.
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 to satisfy the assumption of homoskedasticity. We also calculated temporal turnover for the entire dataset.

Modelling
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, reindeer densities, management intensity, and climate. The final dataset thus consisted of 245 resampled sites. These sites were then assigned to five biogeographical districts (BGD, Fig. 1) based on the Finnish biogeographical province division to account for spatially structured confounders in modelling. The original provinces were slightly modified for this study by merging districts with very few plots to neighboring districts, and by moving plots between districts so as to make each district either completely within or outside the reindeer herding area.
We used hierarchical bayesian regression models to analyze plot-level changes in species diversity and composition. Two kinds of models were deployed based on the nature of the response variable.
For species richness, evenness, diversity, and the three functional composition variables, the response variables could be modelled as a function of time, whereas temporal turnover could only be modelled with spatially varying effects, because turnover already represents the difference between two timepoints. We built four models with increasing complexity to answer our stated questions. there were no average changes in species richness, and because species diversity is the product of richness and evenness, trends in diversity could not be strongly driven by changes in richness. Model 4 was thus only fitted for evenness. Furthermore, we could not think of a mechanism by which leaf traits would affect diversity directly. If a correlation between these variables existed, it would probably be related to the evolutionary history of the species pool, and thus be unique to each area. In contrast, average species size is directly related to the productivity and physical structure of the community, with believably generalizable effects on diversity. Model 4 thus includes average species height as the only predictor describing functional composition. All responses were standardized by subtracting the mean and dividing by standard deviation before analyses. The structure for the expectation function of different models is presented in Table 1.
All models assumed a gaussian error distribution. In the models in Table 1 decrease, no change, increase. The effect of no canopy cover change was fixed to zero.
Bayesian hierarchical models were built using the R package for Bayesian modelling rethinking (version 2.11, McElreath, 2020) , an interface to the probabilistic programming language Stan (Carpenter et al., 2017) . In all models, coefficients for continuous variables were given a Normal(0,0.5) prior distribution, coefficients for categorical variables were given a Normal(0,1) prior distribution, all standard deviations were given an Exponential(1) prior distribution, and temporal correlations in the BGD-level random effect were given an LKJ (2)  To test whether vegetation height is related to evenness in time (modelled above) as well as in space, we modelled evenness as a function of height using a generalized additive model (Wood, 2011) . The effect of height was entered as a thin plate spline with basis dimension 20. We also tested whether the height-evenness relationship was different between sampling times (whether there was a time-height interaction), but a comparison of AICs indicated this was not the case.

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

Overall vegetation changes
Across the study area, we detected no changes in species richness, but observed increases in plot-level Shannon diversity and evenness during the sampling interval. There was no detectable change in average species height or SLA, but average LDMC had increased. Average plot-scale turnover between the sampling periods was about 71% (Fig. 2). 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) using the equation for compound interest, corresponds to an annual turnover of 2.2%.

Geographical variation in vegetation changes
Average plot-level vegetation changes were not representative of changes in all the BGDs (Fig. 2).
Increases in species diversity and evenness had a clear geographic trend, with no changes in the south and clear increases in the north. Height changes also had a weak latitudinal gradient, with uncertain increases in height associated with more southern areas, and vice versa for the north. Patterns in species richness and SLA change and turnover were more idiosyncratic with respect to geography.
LDMC was the only community property whose changes did not differ geographically; it increased uniformly in all the BGDs.

Causes of vegetation changes
There was no discernible effect of climate change on any of the studied variables (Fig. 3a).
Increasing reindeer densities were associated with increasing species diversity, evenness and LDMC, and decreased SLA and height (Fig. 3b). In addition, plots within the reindeer herding area had more positive trends in diversity and evenness and more negative trends in height compared to those outside the area.
Forests with signs of intensive management had more negative species richness trends and increased temporal turnover compared to unmanaged forests (Fig. 3c). Unmanaged forests had an uncertain positive richness trend, with a mean change of +0.3 standardized units and an 80% probability of the trend being positive.
Increases in canopy SLA decreased field-layer SLA, and decreased temporal turnover (Fig. 3d).
Canopy SLA changes were almost completely explained by changes in the dominance of the evergreen tree Picea abies (Fig. S2.5), with increasing Picea corresponding to decreasing canopy SLA.
Canopy cover changes also affected these communities, with cover increases and decreases changing communities in different directions (Fig. 3e). The differences in responses to increasing and decreasing canopy cover suggest that denser canopies cause increased diversity, evenness, SLA, and turnover, and decreased LDMC and height.

Effect of height on diversity change
Height had a nonlinear but mostly negative spatial relationship with evenness (Fig. 4a). The temporal relationship between height and evenness was also negative (Fig. 4b). Including height as a predictor somewhat decreased the inferred effects of reindeer and canopy cover changes on evenness ( Fig.   S2.7).

Discussion
As an answer to our first three study questions, we found that temporal trends in plant community properties have regional and fine-scale variation that is partially explained by differences in land-use.
Most evidently, increased reindeer densities decreased average SLA and decreased plant height, and resulted in increased species evenness and diversity. The effects of forestry, in contrast, manifested as decreased species richness and increased turnover, demonstrating that land-use effects on vegetation may be missed if composition, turnover, and diversity are not all taken into account. In addition to land-use change, functional composition also responded to canopy development, because light-availability in the understory affected both vegetation height and leaf traits. Our results show that diversity changes were partially modulated by changes in plant size, illustrating a direct effect of changing functional composition on diversity, and answering our fourth study question. Thus, most observed changes in functional composition were directly linked to land-use and canopy changes in ways that are consistent with theory, and provided complementary information compared to changes in diversity, which highlights the suitability of functional metrics for biodiversity monitoring in the anthropocene (Pereira et al., 2013) and for planning conservation interventions (Keddy, 1992) .

Vegetation trends vary in space
Our analysis revealed that many community characteristics were stable over the 4-5 decades study interval across the study area, but displayed trends when vegetation changes were conditioned on geographic location. Apparent stability at larger scales can thus be the sum of opposing trends in different regions. However, we report that plot-scale species richness was stable as in a previous meta-analysis (Vellend et al., 2013) , even though plot-scale species turnover was very high ( ca. 70%).
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 readily replaced by new species from the metacommunity (Supp & Ernest, 2014) . Furthermore, despite high plot-scale turnover, most likely attributable to high species richness (Table S2.3), a long sampling interval, and inherent relocation error, 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) .

Land-use and canopy dynamics change the composition of field-layer vegetation
Reindeer husbandry was connected to more community change metrics than forestry. Our results strongly suggest that increased reindeer density selectively diminished the relative abundance of large plants with high SLA and thereby changed the functional composition of the communities. These strong effects of reindeer grazing are consistent with experiments in Scandinavian mountain forests (Sundqvist et al., 2019) . 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) . Supporting this general pattern, our results provide clear evidence that reindeer control forest community composition by trait-based selective herbivory. Such evidence has so far been very limited, especially on longer timescales. Furthermore, reindeer presence seemed to affect vegetation changes even when animal densities remained the same, since changes in vegetation height were more negative in BGDs that were inside the reindeer herding area. However, this effect could also be related to other confounding factors that change along the latitudinal gradient.
Forest management effects on composition were mostly direct i.e. not mediated by management effects on canopy cover and composition, as evidenced by the low difference between management effects in models that included vs. excluded canopy variables. Managed forests had higher turnover, which agrees with previous studies from other biomes that have reported accelerated biodiversity change following human disturbance (Barlow et al., 2016;Lake et al., 2000) . However, changes in functional composition between unmanaged and managed forests did not differ markedly, despite management affecting species richness and turnover. Thus, long-term species losses and gains caused by forest management do not seem to be strongly controlled by plant height or leaf traits, but some other properties of plants. However, the most intense management intensity class that included recent clearcuts had more negative SLA and more positive LDMC trends compared to other managed forests, which supports previous findings of harvesting having a signal on the light-interception strategy of the understory community (Tonteri et al., 2016) , but suggests that this effect is transient and undetectable in later forest development phases.
Canopy development also affected functional composition. Canopy cover had on average increased (see Fig.S2.6), leading to increased shading and decreased abundance of tall, light-demanding species with high SLA and low LDMC. SLA is a succession trait, which increases in response to shade (Dahlgren, Eriksson, Bolmgren, Strindell, & Ehrlén, 2006) , whereas LDMC has probably decreased because of its negative covariance with SLA (Bruelheide et al., 2018) . Canopy composition, here reflecting the dominance of Picea (see Fig.S2.5), also influenced the functional composition of the understory. However, the effects on understories can also be partly due to some other quality of Picea

Changes in diversity are partially mediated by vegetation height
Changes in species diversity were mostly caused by changes in evenness, and those in turn were directly related to changes in vegetation height and its environmental drivers, namely reindeer herbivory and changes in canopy cover. Furthermore, the relationship between evenness and plant height was consistent both in space and time, providing further evidence that changes in functional composition can be a direct mechanism behind diversity changes. Communities composed of larger species often have fewer plant individuals and exhibit decreased evenness due to the dominance of large plants (Oksanen, 1996) . Reindeer grazing and canopy cover effects are most likely a result of this coordination: tall species are being eaten or suppressed by lack of light. This in turn reduces dominance, which increases diversity (Grime, 1973) . In these forests, both canopy cover and reindeer densities have increased on average. According to our results, both these changes favor shorter plants and result in increased evenness and diversity, thus partially explaining the observed large-scale changes in diversity. These results agree with Depauw et al. (2020) , who also found that increased shading is responsible for increased evenness in the herb-layers of temperate European forests.
The effects of reindeer and canopy variables on evenness were not fully mediated by vegetation size (Fig. S2.7). One explanation may be that our study does not account for intraspecific changes in height, which could happen in response to altered light availability or grazing pressure (Jessen, Kaarlejärvi, Olofsson, & Eskelinen, 2020) . Thus, the effects of reindeer herbivory and canopy cover on evenness via vegetation height are probably even stronger than detected here.
Findings discussed above are summarized in Fig. 5.

Conclusions
Our study demonstrates that large-scale vegetation changes are not necessarily caused by any single driver, and that observed change or lack thereof can be a function of many, sometimes opposing trends. However, we show that functional traits can be highly useful in connecting vegetation changes to local and regional drivers such as land-use and canopy changes, as well as explaining detected diversity changes in a globally comparable way. In conclusion, diversity, functional composition, and species turnover provide complementary information on how communities respond to environmental changes, which underlines the importance of using multiple metrics to monitor the effects of global change on biodiversity.  Table 1: Structure of the expectation function of models used in this paper. All models had a gaussian error distribution. Community property refers to standardized species richness, species diversity, species evenness, and the abundance-weighted averages of log-transformed SLA, LDMC and height. Turnover refers to standardized logit-transformed temporal Jaccard distance. See methods for descriptions of the different parameters.

Fig. 2:
Temporal changes in boreal herb-rich forest understory species richness (RIC), Shannon diversity (DIV), species evenness (EVE), lead dry matter content (LDMC), specific leaf area (SLA) and height (HEI), and average turnover (DIS), across the biogeographical districts (BGD). The points are medians, the narrow and thick lines are 2/3 and 20/21 credible intervals, respectively, corresponding to 2/1 and 20/1 odds for the true parameter value to lie within the given interval, or 66.67% and 95.24% probability. The gray violin plot in the background is a kernel density plot of the posterior distribution of average change across all plots as inferred from model 1. All response variables except DIS are shown on a standardized scale. BGDs are arranged according to their latitude, from south (1) to north (5).

Fig. 3:
Standardized direct effects of climate change (a), reindeer (b), forest management intensity (c), canopy SLA increase (d), and canopy cover change (e) on community-level properties and turnover in boreal herb-rich forest understories. Reindeer effects are displayed as the direct effect of reindeer density change and the combined effects of density change and being situated inside the reindeer herding area (average difference in trend compared to a plot outside the area). Management intensity effects are contrasts against unmanaged forests (intensity class 1). The points are medians, the narrow and thick lines are 2/3 and 20/21 credible intervals, respectively, corresponding to 2/1 and 20/1 odds for the true parameter value to lie within the given interval. Response abbreviations as in Fig. 2. Traits were log-transformed and turnover logit-transformed before analyses. Fig. 4: Spatial (a) and temporal (b) relationships between community evenness and average plant height. The spatial relationship is an unadjusted GAM of evenness against height. The temporal effect is a marginal standardized effect from hierarchical linear model 4 (see methods). In a), the line is the mean effect and the shading is the 20/21 confidence interval. In b), the point is the median effect, and the narrow and thick lines are 2/3 and 20/21 credible intervals, respectively. Refer to methods for more information.