Nitrogen deposition causes eutrophication in bryophyte communities in central and northern European forests

Abstract

Key message

Our results indicate that nitrogen deposition is likely to adversely affect forest bryophyte communities, having negative impacts in terms of increased dominance of nitrophilic species at the expense of N-sensitive species and a decrease in evenness.

Context

Elevated atmospheric deposition of nitrogen (N) has long been recognised as a threat to biodiversity and, despite declines in European emission levels, will remain a threat in the future.

Aims

It has proven difficult to show clear large-scale impacts of N deposition on vascular forest understorey species, and few studies have looked at impacts on forest bryophytes. Here, we assess the impact of nitrogen deposition on forest bryophyte communities.

Methods

We used data from 187 plots included in European monitoring schemes to analyse the relationship between levels of throughfall nitrogen deposition and bryophyte taxonomic and functional diversity and community nitrogen preference.

Results

We found that nitrogen deposition is significantly associated with increased bryophyte community nitrogen preference and decreases in species evenness.

Conclusion

Our results indicate that nitrogen deposition is likely to adversely affect forest bryophyte communities, having negative impacts in terms of increased dominance of nitrophilic species at the expense of N-sensitive species and a decrease in species evenness.

1 Introduction

Elevated atmospheric deposition of nitrogen (N) has long been recognised as a threat to biodiversity and, despite declines in European emission levels, will remain a threat in the future (Dirnböck et al. 2018). Numerous experimental studies have clearly shown that nitrogen addition can result in changes in forest understory vegetation (Bobbink et al. 2003; Nordin et al. 2005; Gilliam et al. 2016). Many observational studies have also demonstrated atmospheric N deposition mediated changes in forest floor vegetation at stand level, an early example being Thimonier et al. (1992) who showed an increase of nitrophilous species. Further studies found N effects, for example increased homogeneity of understory vegetation was found using 20 years of data from a monitoring site in Austria (Hülber et al. 2008), and an analysis of Swedish data over a similar time period demonstrated a eutrophication effect seen in the decline in dwarf shrubs and increase in mesotrophic forbes and ferns found in the temperate zone (where deposition levels were around 12.5 kg N ha⁻¹ yr⁻¹) but not in the boreal zone (deposition levels around 2 kg ha⁻¹ yr⁻¹) (Hedwall and Brunet 2016). Similar effects have been found at landscape level, such as the decline in several indices of plant diversity at high N deposition levels found at sites across Switzerland (Roth et al. 2015). Nitrogen has also been found to have strong effects on tree growth in terms of stand volume (Etzold et al. 2020), and to negatively impact ectomycorrhiza (Lilleskov et al. 2001, 2002; van der Linde et al. 2018) and epiphytic lichens (Giordani et al. 2014), amongst other impacts.

Despite the above results showing eutrophication effects at stand and landscape scales, it has proven difficult to detect broad-scale (e.g. European level) N deposition effects on forest floor vegetation comparable to the negative effects of N deposition on grassland species richness that have been demonstrated at continental scale (Stevens et al. 2010). The last decade has seen several studies that have aimed to demonstrate such an effect in forest understorey, generally focusing on vascular plants. Studies using Europe-wide data from long-term monitoring schemes have found some evidence of eutrophication effects on vascular plant community composition (Seidling and Fischer 2008; Seidling et al. 2008) and a significant but small shift towards nitrophilous species attributable to anthropogenic N deposition (van Dobben and de Vries 2016). Another study of seminatural European temperate woodland understory using long-time scales found significant shifts in species composition but no directional change that could be attributed to N deposition (Verheyen et al. 2012). While Dirnböck et al. (2014) also found no Europe-wide trend across all sites, by focusing on sites where the critical load for N deposition is exceeded significant relationships emerged between increased N deposition above critical load levels and decreasing cover of oligotrophic plants, and a higher proportion of nitrophilous plants were found among newly occurring species in surveys. Recently, Staude et al. (2020) found a broad-scale N-driven exchange of species with small range size by species with large range sizes.

It has been hypothesised (Jonard et al. 2015; Gilliam 2019) that N deposition can result in increased canopy growth which in turn supresses responses to increased N availability in forest floor vegetation due to increased light limitation (Binkley and Högberg 2016). Bryophytes however are generally less sensitive to shading than vascular plants (although there is variation between species (Marschall and Proctor 2004)). Light-response curves for a range of bryophytes show saturation of photosynthesis at around 20% of full sunlight (Rastorfer 1970; Kallio and Heinonen 1975; Alpert and Oechel 1987). Consequently, a eutrophication signal might be easier to find among forest floor bryophytes than among vascular plants in the same habitat if increased shading is indeed inhibiting the response of vascular plants.

That most studies focus on vascular species is understandable given that bryophyte data is sparse in monitoring schemes. Non-vascular plants are, however, an important part of the understory plant community in many European forests (Turetsky 2003). In addition, the basic biology of vascular plants and bryophytes suggests that they should show different sensitivity to changes in the environment originating in the atmospheric deposition of pollutants, with clear reasons to expect a stronger response in bryophytes. Since bryophytes absorb water through their above-ground surfaces rather than via roots in soil, they are more sensitive than vascular plants to the chemical composition of rainfall (Bobbink et al. 2010). Thus, changes in bryophyte composition can be particularly strong with additional N, with some species being highly sensitive (Nordin et al. 2005; Strengbom et al. 2001), and these changes can be very long lasting, even with no further N inputs (Bobbink et al. 2003). Strong responses of bryophytes have also been shown where vascular plants at the same locations have failed to respond to N addition (Mäkipää 1998), while a resurvey study after 33 years found greater community N affinity in bryophytes but not vascular plants (Becker Scarpitta et al. 2017). Given the combination of bryophyte sensitivity to N deposition and their general shade tolerance, evidence of eutrophication effects at a large geographic scale may be more detectable than among vascular plants.

While changes in both community N preference and taxonomic diversity can be expected, a complementary investigation of functional traits can also be informative, particularly where large geographic areas are involved. Such an analysis of functional diversity can be helpful when considering responses across areas which incorporate different species pools (Diaz et al. 2004). In N-limited conditions, there is considerable spatial heterogeneity in N availability, contributing to plant diversity. Long-term N deposition however implies a more homogenous availability of nitrogen and a corresponding homogenisation of the plant community (Gilliam et al. 2006; Hülber et al. 2008). We therefore hypothesised that diversity measures would be negatively associated with N deposition.

Long-term monitoring schemes that regularly assess forest floor vegetation and levels of N deposition are a rare and valuable resource for investigating these questions. Beginning in the late 1970s, concerns about the potential impact of air pollution on European forest ecosystems led to signing of the United Nations Economic Commission for Europe (UNECE) Convention on Long-range Transboundary Air Pollution (CLRTAP or Air Convention). During the 1980s, the International Co-operative Programme on Assessment and Monitoring of Air Pollution Effects, including ICP Forests (Ferretti 2013) and ICP Integrated Monitoring (ICP-IM 1998), was installed. By using resources from the ICP Forests and ICP IM monitoring schemes such as data on bryophyte abundances and throughfall deposition, an investigation of the potential impacts of N deposition on bryophyte taxonomic and functional diversity, and on community N preference, is possible, including an investigation of whether a eutrophication signal can be seen in community-weighted mean Ellenberg values. We hypothesised that N deposition would be associated with a shift in bryophyte community composition towards more nitrophilous species and a decrease in diversity.

2 Material and methods

The evaluation was based on data that was collected by the UNECE ICP Forests Network and the ICP IM Network. N deposition data is based on measurements of throughfall deposition as this is the most realistic measurement of deposited N available to understorey vegetation. Ground vegetation surveys are undertaken at least once every 5 years. Bryophyte abundance is not a core focus of the monitoring, and at some locations, bryophytes are a very minor part of the ground cover. Consequently, bryophyte abundance data has not always been recorded in the monitoring schemes, and some sites which participate in the schemes cannot be included in the present study for this reason. Also, there is some concern about the quality of bryophyte data (e.g. Seidling et al. 2020; Vittoz et al. 2010). In particular, large-scale data on bryophytes can be affected by several sources of error, i.e. the inherent difficulty in detection (specimens of some moss species can be tiny and often with small cover — Seidling et al. 2020), uncertainty with respect to taxonomy and proper identification (Molina and Marcot 2007) and observer skills and bias. Some of these are inherent in long-term vegetation monitoring whether of bryophytes or vascular plants, where teams of observers change between sites and at the same sites over time. Unfortunately, there is limited information on the actual quality of bryophyte data at ICP Forests and IM sites (see Seidling et al. 2020), but it is likely that data used in this paper can be affected by the above-mentioned error sources (for reference, the data included here have a mean species richness of 9.2 and a coefficient of variation of 68%). However, it is also likely that errors have a random distribution, and this may reduce potential major bias in the model outputs. Nonetheless, potential issues in data quality must be considered when evaluating the results.

The monitoring schemes also cover a very wide geographical area, and we deliberately excluded sites located in the Mediterranean zone as a biome with distinct climate, forest structure and species composition, using the European Environment Agency biogeographical regions dataset containing the official delineations used in the Habitats Directive (92/43/EEC). Full details of the methodology used in all aspects of the monitoring programmes are available (ICP-IM 1998; ICP-Forests 2016).

Despite these limitations, data on N throughfall deposition levels (NO₃ and NH₄), bryophyte abundance (ground cover, expressed as a proportion) and the other variables included here were available from 187 plots spread across mostly central and northern Europe (Fig. 1), covering the period 1994–2016 (for a total of 537 plot surveys). Although there are currently 545 plots in ICP Forests recording atmospheric deposition, many were unable to be included due to lack of other data required here. Note that the participation of countries and the maintenance of plots have not been consistent over this time period, with some dropping out and others joining later, meaning that data for many plots is available only for a subset of the years covered. The number of observations per year is also not constant, and many sites keep to the same 5-year cycle for repeat samples resulting in peaks in the number of observations in 1995, 2000, 2005, etc. The sites used cover a gradient in N deposition, from high levels (e.g. over 30 kg/ha/year) in some parts of Central Europe to very low levels (e.g. under 1 kg/ha/year) in northern Scandinavia. A few observations showed exceptionally high levels of N deposition (presumably due to a local point source) and were excluded as outliers by removing the 100th percentile of the data based on N throughfall deposition levels (i.e. all data < = 0.99 is included). One plot with exceptionally high combinations of both NH₄ and NO₃ values was further removed as an outlier due to a disproportionately increased bias in the subsequent additive quantile regression models and to avoid overfitting in areas with very low sample density. This resulted in removing 3 observations (of 540, leaving 537). Where multiple measurements per year were available for a variable, these are aggregated to an annual mean value.

Fig. 1

Location of plots used

The effect of N deposition on bryophyte communities was investigated by calculating the species richness and evenness (Pielou evenness (Pielou 1966), as recommended in, e.g., Jost 2010) for each plot/year combination. We also calculated the community-weighted mean Ellenberg N value for each plot/year combination using the vegdata R package (Jansen and Dengler 2010) with values obtained from the BryForTrait database (Bernhardt-Römermann et al. 2018). Ellenberg N values are a classification of species on a 9-point scale according to the location of their ecological optimum along a nutrient (mainly N availability) gradient (Ellenberg et al. 2001). As a fourth response variable, we used functional diversity. This was calculated as a community-weighted mean value for Rao’s quadratic entropy (which is a measure of the dispersion of species, weighted by relative abundance, in a multidimensional functional trait space (Botta-Dukát 2005), based on three broad morphological traits (growth form, life form, and life strategy). Trait values were acquired from the BryForTrait database (Bernhardt-Römermann et al. 2018), with Rao’s quadratic entropy calculated using the “FD” R package (Laliberté et al. 2014).

To assess relationships between deposition and the bryophyte community response variables, we used smooth additive quantile regression models (R package qGAM) (Fasiolo et al. 2020) that fit the desired response by smooth functions of predictor variables. As qGAMs are not yet commonly used in ecology, we include a relatively detailed description here (see Appendix), but readers are referred to Fasiolo et al. (2020) for full details. One main benefit of the approach is that (as with generalised additive models, GAMs) the effect of a predictor variable on the response is not limited to a linear relationship. In comparison with classical regression models like GAMs however, these nonparametric models allow us to fit any quantile of interest, and do not require a pre-set error distribution and perform well even with heteroscedasticity. Furthermore, quantile regression is more robust to outliers compared to classical regression models (Koenker 2005). For this study, we fit the median response. The syntax of such models is equal to classical generalised additive models (Wood 2012) implemented by the well-established R-package mgcv.

This offers high flexibility in model complexity and makes it possible to include terms that account for spatiotemporal autocorrelation or complex non-linear interactions. Details of the calculation of p-values based on Wald tests in these models can be found in Wood (2012).

In this study, we model the spatiotemporal autocorrelation of the observations (i.e. the same locations being resampled over time and the similarity of locations nearby one another) by including a tensor term of survey year, latitude and longitude of samples (Wood 2006) in the models. The models include the throughfall nitrogen deposition of the two forms of nitrogen NO₃ and NH₄ (Table 1). These are modelled as an interaction smooth, as their interaction is likely to be ecologically relevant (i.e. we expect bryophyte communities will not necessarily react in the same way to NH₄ and NO₃). Potential colinearity among variables was checked using the concurvity function of the mgcv package. Concurvity is essentially a generalisation of colinearity which is applicable to GAMs, causing similar problems of interpretation as colinearity in a multivariate linear regression (Buja et al. 1989). Explanatory variables showing worst-case concurvity over 0.8 were further investigated by refitting models with those terms excluded. For all smooth terms, we used thin plate regression splines (Wood 2003).

The variables considered beyond N deposition were chosen based on factors which are well known to affect vegetation composition at large spatial scales: annual mean temperature and precipitation, light availability (i.e. canopy cover, which can exceed 100% due to overlapping levels), forest age, soil pH and forest type (i.e. broadleaved or coniferous). Forest type is modelled as a fixed effect as it is not continuous and has only two levels, making it unsuitable for inclusion as a random effect.

Model selection was done through comparison of AIC values to determine the inclusion of the fixed effects followed by the use of an additional penalty option for the smooths, which allows for variables which do not improve the model to be effectively penalised to zero (the model intercept), which can be interpreted as zero or at least negligible effect of that variable (Marra and Wood 2011). As this variable selection is an integral part of the model fitting process, nonsignificant terms (based on Wald tests) are present in the final models but with effects penalised towards zero. In all cases, model fitting was followed by checks of goodness of model fit via diagnostic plots.

The final models are as follows:

$${q}_{0.5}\left(Ellenberg\ {N}_{i}\right)\sim {f}_{1}\left({NH_4}_{i},{NO_3}_{i}\right)+{f}_{2}\left({Longitude}_{i},{Latitude}_{i},{Year}_{i}\right)+{f}_{3}\left({precipitation}_{i}\right)+{f}_{4}\left({canopy}_{i}\ \right)+{f}_{5}\left({age}_{i}\right)+{f}_{6}\left({pH}_{i}\right)+{\psi foresttype}_{k_i}$$

$${q}_{0.5}\left({Evenness}_i\right)\sim {f}_1\left({NH_4}_i,{NO_3}_i\right)+{f}_2\left({Longitude}_i,{Latitude}_i,{Year}_i\right)+{f}_3\left({precipitation}_i\right)+{f}_4\left({canopy}_i\ \right)+{f}_5\left({age}_i\right)+{f}_6\left({pH}_i\right)+{\psi foresttype}_{k_i}$$

$${q}_{0.5}\left({Richness}_i\right)\sim {f}_1\left({NH_4}_i,{NO_3}_i\right)+{f}_2\left({Longitude}_i,{Latitude}_i,{Year}_i\right)+{f}_3\left({precipitation}_i\right)+{f}_4\left({canopy}_i\ \right)+{f}_5\left({age}_i\right)+{f}_6\left({pH}_i\right)+{\psi foresttype}_{k_i}$$

$${q}_{0.5}\left({RaoQ}_i\right)\sim {f}_1\left({NH_4}_i,{NO_3}_i\right)+{f}_2\left({Longitude}_i,{Latitude}_i,{Year}_i\right)+{f}_3\left({precipitation}_i\right)+{f}_4\left({canopy}_i\ \right)+{f}_5\left({age}_i\right)+{f}_6\left({pH}_i\right)+{\psi foresttype}_{k_i}$$

where q_0.5 is the median, Ellenberg N is community-weighted mean Ellenberg N value, RaoQ is functional diversity, Evenness is Pielou evenness, Richness is species richness, NH_4i is throughfall NH₄, NO_3i is throughfall NO₃, Longitude_i is longitude, Latitude_i is latitude, Year_i is year, precipitation_i is precipitation, canopy_i is canopy cover, age_i is stand age, pH_i is soil pH, and ψforesttype is forest type. f₁, etc. indicate smooth functions.

3 Results

In the following sections, we present the results of the qGAMs, analysing how Ellenberg N, species evenness and richness and functional diversity are affected by our explanatory variables. Grey areas of the heat maps are N deposition levels not present in the data. Note that the fitting procedure effectively penalises the effect of a term that does not improve the model to zero. In prediction of the partial effect of this term on the response scale, this leads to a flat line at the model intercept.

3.1 Community mean N preference

The mean N preference of the bryophyte community changes significantly with deposition levels of both NH₄ and NO₃ (p < 0.001), although the gradient is not even, indicating a variation in community responses to the two forms of nitrogen. However, the strongest increase in Ellenberg N values is seen with substantial deposition levels of both N forms combined (upper red parts of Fig. 2a). The relationship between N preference and shading is not significant and is effectively selected out of the model by the penalty term (Fig. 2b). There is a negative relationship between N preference and precipitation (p < 0.001, Fig. 2c). Forest age is significantly related to N preference (p < 0.001, Fig. 2d), with higher values in younger forests (< 75 years), and coniferous forests have significantly (p < 0.001) lower N values than broadleaf-dominated forests (Table 1).

Fig. 2

QGAM smooths for mean Ellenberg N (y-axes), with 95% confidence intervals indicated. In all plots, grey areas of the heat map are N deposition levels/combinations not present in the data. a Heat map of the smooth term showing the relationship between Ellenberg N and N deposition levels (kg N/ha/yr), b smooth for canopy cover (%), c smooth for precipitation (mm), d smooth for forest age (years), and e smooth for soil pH

Table 1 Model results for community mean N preference (Ellenberg N). Approximate significance of the intercept, any parametric terms and the smooth terms for single variables and tensor product smooth terms for interactions

3.2 Species evenness

Pielou evenness shows a significant decline with increasing levels of NH₄ and NO₃ (p < 0.001) (Fig. 3a). Evenness shows a significant positive relationship with canopy cover (p < 0.001, Fig. 3b). There is a positive relationship between diversity and precipitation (p = 0.004, Fig. 3c), although the confidence intervals are wide and only just clear the median, indicating a weak relationship. Forest age is not significantly related to evenness (Fig. 3d) and is penalised to a flat response. The apparent negative relationship with pH is also not significant (Fig. 3e). Finally, we found that coniferous forests have significantly (p < 0.001) lower evenness than broadleaf-dominated forests (Table 2).

Fig. 3

QGAM smooths for Pielou eveness. All other details as in Fig. 2

Table 2 Model results for species evenness. Approximate significance of the intercept, any parametric terms and the smooth terms for single variables and tensor product smooth terms for interactions

3.3 Species richness

Species richness does not change significantly with N deposition (Fig. 4a). There is a significant negative relationship with canopy cover (p < 0.001, Fig. 4b), although there is an increase at the highest levels of cover (albeit with very wide confidence intervals in this area). There is a positive association with precipitation levels (p < 0.001, Fig. 4c) and negative relationships with forest age (p < 0.001, Fig. 4d) and pH (p = 0.008, Fig. 4e). The confidence intervals for pH are wide and only just clear the median however, indicating a weak relationship. Forest type shows a significant difference in species richness (p < 0.001), with coniferous forests having higher richness than broadleaf-dominated forests (Table 3).

Fig. 4

QGAM smooths for species richness. All other details as in Fig 2

Table 3 Model results for species richness. Approximate significance of the intercept, any parametric terms and the smooth terms for single variables and tensor product smooth terms for interactions

3.4 Functional diversity

Rao’s quadratic entropy as a measure of functional diversity apparently decreases with increasing N deposition, although this effect is narrowly nonsignificant (p = 0.07, Fig. 5a, Table 4). As the p-values especially for penalised terms are only approximate, this effect should not be entirely neglected, however. Canopy cover shows an overall positive relationship (p < 0.001) with functional diversity (Fig. 5b). Precipitation, forest age, and pH (Fig. 5c and d) are not significant, and age and pH are modelled with a flat response and effectively removed from the model by the extra penalty. Finally, forest type is significant (p < 0.001), with coniferous forests having lower functional diversity than broadleaf-dominated forests (Table 4).

Fig. 5

QGAM smooths for Rao’s quadratic entropy as an index of functional diversity. All other details as in Fig. 2

Table 4 Model results for functional diversity (Rao’s Q). Approximate significance of the intercept, any parametric terms and the smooth terms for single variables and tensor product smooth terms for interactions

4 Discussion

We modelled data from bryophyte communities across central and northern Europe taking into account the effects of differences in canopy cover, forest age and type, temperature, pH and precipitation. We found that N deposition is significantly associated with increased bryophyte community mean Ellenberg N values, and decreased species evenness, on a European scale. In terms of effects, we find a decline of at most ca.15% in evenness attributable to N deposition (as seen in the heat map, Fig. 3a). It may be useful to compare this to the coefficient of variation (the ratio of the standard deviation to the mean, which indicates the extent of variability around the mean in the data) of 41%. The impact of N deposition on mean Ellenberg N preference is at most a ca. 20% increase (Fig. 2a), compared to a coefficient of variation of 27%.

Since different bryophyte species show varying responses to increased N, a shift in community-weighted mean Ellenberg N value is an expected result of higher levels of N deposition. The varying responses underlying community shifts can, for example be seen in Salemaa et al. (2008) where experimental addition in laboratory conditions of NH₄ and NO₃ to three common forest bryophyte species found that all had a unimodal response but with different optima and width of tolerated range. Experimental additions in field studies have also found changes in community composition, with an increasing dominance of nitrophilous species at the expense of others (Gordon et al. 2002). High turnover of species within bryophyte communities (Zechmeister et al. 2007) suggests that the time frame of the current study is more than adequate to capture changes in taxonomic and functional diversity that we expected to see due to the rapid life cycle of many bryophyte species (Turetsky 2003).

Previous studies focussed on vascular forest understory plants have found only limited evidence of eutrophication effects or effects in subsets of sites or species (Seidling and Fischer 2008; Seidling et al. 2008; Dirnböck et al. 2014; van Dobben and de Vries 2016; Roth et al. 2020). Several reasons have been proposed to explain these limited effects. Much of Europe has experienced decades of elevated N deposition, and high levels of background deposition can make it difficult to provoke changes with N additions. Gilliam et al. (2006) for example explain the lack of response in a temperate hardwood forest by referring to the 10 kg N ha⁻¹ yr⁻¹ wet N deposition, a suggestion supported by Hedwall et al. (2013) who found a reduction in the effect of N addition in experiments performed in higher N deposition areas. While strongly N-limited sites (and acidic soils (Simkin et al. 2016)) are particularly vulnerable to N deposition-induced vegetation changes, N-rich sites are relatively unaffected by the small size of N depositions compared to the high rate of N mineralization and/or the large existing N pool (Diwold et al. 2010). Given the prolonged deposition of nitrogen over much of Europe, it is to be expected that many sites are affected by this phenomenon, hindering efforts to find continental scale eutrophication effects.

A decrease in light availability due to fertilisation-induced increase in canopy cover has also been suggested as an explanation, with N uptake by canopy species resulting in increased shading which in turn suppresses the response of understorey species (Jonard et al. 2015; Gilliam 2019). Consequently, attempts to find an N-driven response in understorey plants should take into account changes in the overstorey (Binkley and Högberg 2016), as under closed canopies light is often the most important limiting factor (Gilliam and Roberts 2003). We found that a eutrophication effect can be seen even when accounting for changes in light regime, but light limitation is still a plausible explanation for the relatively small size of some of the observed effects. Light may also play a role in the unexpected negative relationship between stand age and richness if the lower richness in broadleaf forests (Table 3) is driven by light limitation in old dense stands. A negative relationship emerged between N preference and stand age, which may be partly explained by a higher N availability in young stands where these are created by extensive cutting or other disturbance (Hannerz and Hånnel 1997; Prescott 1997).

Not only the amount but also the form of N deposition can be important. Boreal plant species for example are adapted to low N soils with slow N mineralisation rates where NO₃− is rare and the N available for uptake is generally organic forms of N. Consequently, many low N-adapted species are only able to make use of deposited NO₃ to a limited extent, while nitrophilous competitors are well adapted to make use of it, giving them a competitive advantage where nitrate is relatively abundant (Nordin et al. 2001; Bobbink et al. 2010). While NH₄ shows a slightly weaker effect than NO₃ in our models of Ellenberg N (Fig. 2a), the heat map of the interaction of NH₄ and NO₃ clearly shows that both forms should be considered.

While no change in species richness was found that could be attributed to N deposition, a clear decline in evenness with increasing N was found, indicating an increasing dominance of locally common species. This also suggests that the increase in community-weighted mean Ellenberg N value is generally a result of the increased dominance of existing nitrophilous species. Given the issues noted in the methods section regarding data quality, it is likely that there are species missing which are difficult to find or identify, a problem which could particularly affect the species richness analysis, where we acknowledge that a loss or a gain of such species may have occurred undetected. However, while these may be of particular conservation concern, they are unlikely to be present in abundances which could be expected to greatly shift the community-weighted mean N preference. These results should generally be interpreted as changes in those species that are more abundant and therefore more easily (and more likely correctly) identified. We hypothesised that N deposition would be associated with an increase in community N preference, which was confirmed. Our hypothesis of decreased taxonomic and functional diversity was only partly confirmed, with an apparent (but narrowly nonsignificant) decrease in functional diversity, and no change in species richness, but a decrease in species evenness.

5 Conclusion

While N deposition levels are declining in most of Europe, they continue to threaten biodiversity and are likely to do so into the future. Our results show that N deposition is likely to adversely affect forest bryophyte communities, having negative impacts in terms of increased dominance of nitrophylic species at the expense of N sensitive species. Although no evidence of an impact on species richness was found, species evenness also decreased with N deposition.

Appendix

1.1 qGAMs

Generalised additive models (GAMs) are an advanced regression technique where the response variable is fitted by a sum of smooth functions of the covariates (Wood 2017). In contrast to linear regression, the smooth functions could have any possible but continuous shape and level of complexity and thus are not limited to linear relationships only. A trade-off between the level of complexity of each smooth term is controlled by the addition of a penalty usually based on the second derivative multiplied by a respective smoothing parameter to avoid overfitting (Wood 2017). A very large smoothing parameter would result in a penalization towards a straight line, while a very low smoothing parameter would lead to a very wiggly curve. Objective and automated smoothing parameter estimations are available using prediction error criteria or likelihood-based methods like in popular “R” packages such as “mgcv” to automatically pick the best smoothing parameter. The easiest way to generate a smooth function of a covariate is to use piecewise polynomials placed and aligned at knots across the covariate range. More advanced techniques such as thin plate regression splines (Wood 2003) avoid the selection of the position of knots and are the default in packages such as “mgcv”. Sometimes, the interaction between multiple covariates has an additional effect on the response variable. Smooth functions can be extended to more dimensions so that the potentially non-linear main effects and interactions are incorporated into the model. If the units of the covariates are the same, this is modelled as extension of, e.g. a thin plate spline into multiple dimensions and called an isotropic smooth. Full tensor product smooths instead are multivariate smooths that are calculated differently and have the advantage that the covariates can have different units (Wood 2006). In this study, we have used a tensor product of the spatial coordinates and year to account for spatiotemporal autocorrelation, thus including the effect of space and time and their interaction and an isotropic smooth for the combined effect of NH₄ and NO₃.

Generalised additive models and most other regression techniques in general fit the conditional mean or expected value under given covariate combinations. In the fitting process, the assumption of a conditional error distribution, e.g. the Gaussian distribution, is needed to estimate the expected value and minimizing model deviance. Although there are plenty of different conditional distributions available to be used within GAMs, it is often difficult to find one that leads to a sufficiently good model (the appropriateness of the conditional distribution can be checked, e.g. via Q-Q plots of the often transformed residuals). Additive quantile regression models (qGAMs) represent a nonparametric alternative to this and allow to fit any conditional quantile of interest instead of the conditional mean (Fasiolo et al. 2020). This means that no choice of a conditional error distribution needs to be made. Fitting is achieved with a smooth generalization of the “pinball” loss in a Bayesian framework (Bissiri et al. 2016; Fasiolo et al. 2020). To check the appropriateness of the model, alternative methods have been created as, e.g. Q-Q plots are not useful here. If the conditional median is fitted for example, 50% of the observed values should be above the fit. Confidence intervals based on the binomial distribution are created across bins of n observations to test for significant deviations (Fasiolo et al. 2020).