How do invasive trees impact shrub layer diversity and productivity in temperate forests?

Abstract

• Key message

Invasive tree species alter taxonomic diversity and functioning of forest shrub layers: Prunus serotina increases shrub layer biomass two to three times but decreases its biodiversity, Robinia pseudoacacia slightly increases shrub layer biomass and has no effect on its biodiversity, while Quercus rubra both biomass and biodiversity of the shrub layer.

• Context

Although the impact of invasive trees on understory biodiversity is known, very little data exist about their influence on shrub layer biodiversity and productivity.

• Aims

To assess impacts of Prunus serotina Ehrh., Quercus rubra L., and Robinia pseudoacacia L. on shrub layer aboveground biomass, species composition, and alpha diversity.

• Methods

We measured stand structures in a set of 168 study plots established in Wielkopolski National Park (W Poland), and we compared biomass and diversity metrics using generalized mixed-effects linear models.

• Results

We found the lowest aboveground biomass of shrub layers in Q. rubra forests. P. sylvestris forests invaded by P. serotina had two to three times higher aboveground biomass than non-invaded forests. R. pseudoacacia forests had 27.8% higher shrub layer biomass than Quercus-Acer-Tilia forests. We found negative impacts of Q. rubra and negligible impacts of R. pseudoacacia on shrub layer alpha diversity metrics. However, the effect of Q. rubra was similar to native F. sylvatica. P. serotina negatively affected functional diversity, but its effects were lower in rich P. sylvestris forests than in poor P. sylvestris forests.

• Conclusion

The introduction of alien tree species alters ecosystem services and species diversity of shrub layers. The direction and magnitude of these alterations are alien species-specific and context-dependent. Therefore, their management should account for their impacts.

1 Introduction

Shrubs play an important role in forest ecosystems. Due to functional distinctiveness and ability to grow beneath forest canopies, shrubs utilize the available light between understory and canopy layers. Numerous shrubs provide ecosystem services by increasing ecosystem biomass and soil erosion control. Most shrub species provide easily-decomposable litter, influencing nutrient cycling (e.g., Horodecki and Jagodziński 2017). Shrubs also provide food for both large (e.g., Bodziarczyk et al. 2017) and small herbivores (e.g., Karolewski et al. 2013, 2020).

Invasive tree species significantly alter ecosystem services (Castro‐Díez et al. 2019) and biodiversity (Terwei et al. 2016; Šibíková et al. 2019). Such effects might be both positive, e.g., increase of carbon stock due to increased nitrogen availability (Rice et al. 2004), or negative, e.g., decreasing growth of native trees by the more effective acquisition of soil resources by invasive species (Hartman and McCarthy 2007). Moreover, invasive trees can affect biodiversity at various trophic levels (Hejda et al. 2017).

Although tree biomass pools in forest ecosystems are well-recognized, shrub biomass is rarely estimated (Woodbury et al. 2007). Despite the development of allometric models (Conti et al. 2019), forest inventories still provide only scarce information about shrub layer species composition and abundance. For that reason few studies estimated the biomass of shrub layers in forest complexes or at a national scale (Škėma et al. 2015). Lack of measurements has also prevented studies on impacts of invasive trees on shrub layers, with some exceptions (Slabejová et al. 2019; Mikulová et al. 2019). Therefore, we aimed to quantify the impact of three invasive woody species: Prunus serotina Ehrh., Quercus rubra L., and Robinia pseudoacacia L. on shrub layer biomass and functional, phylogenetic, and taxonomic diversity. We hypothesized that (1) the impact of P. serotina will differ from Q. rubra and R. pseudoacacia, as P. serotina occupies the shrub layer niche, in contrast to the canopy-dominants Q. rubra and R. pseudoacacia, and (2) the invasive species studied will decrease both biomass and taxonomic, phylogenetic, and functional alpha diversity of shrub layers.

2 Material and methods

2.1 Study design

We conducted our study in Wielkopolski National Park (WNP; W Poland; 52° 16′ N, 16° 48′ E), using a set of study plots designed to assess the spread of natural regeneration of the studied species, and described in detail in previous studies (Dyderski and Jagodziński 2019a). The climate in WNP is temperate, with a mean annual temperature of 8.4 °C and mean annual precipitation of 521 mm, for the years 1951–2010. Forests of WNP were strongly transformed by former forest management, replacing mixed and broadleaved forests with monocultures of Scots pine. Before the establishment of the national park in 1957, numerous alien trees and shrubs had been introduced in WNP, resulting in the highest number of alien woody species among all national parks in Poland (Gazda and Szwagrzyk 2016). In WNP we arranged a set of 168 study plots with an area of 150–2000 m² (mean 667 ± 26 m²). The area of study plots was limited by canopy homogeneity—we aimed to cover the largest possible area within the stands studied, to stabilize the stand structure measurements. We distributed study plots across 21 blocks (Fig. 1), covering one or two pairs of 100 m² formerly established square-plots for natural regeneration assessment (Dyderski and Jagodziński 2019a). Each block was established with a central plot within a monoculture of the invasive species studied (in the case of P. serotina–P. sylvestris monocultures with understories dominated by P. serotina). The plots used in this study covered two or four regeneration plots (pairs of square-plots), depending on the stand homogeneity: when two pairs were located in a homogenous stand, we established a single large plot, otherwise—two separate plots (Dyderski and Jagodziński 2019a). We excluded three study plots located in non-forest vegetation.

Fig. 1

Scheme for a block of experimental plots in the field. Each of the 21 blocks covers a set of 18 square plots (100 m²), established for natural regeneration assessment (see Dyderski and Jagodziński 2018 for details), marked as grey squares. These plots were covered by stand structure plots (bold rectangles), established to cover homogenous forests. For that reason, a single study plot can cover two or four regeneration plots

We divided study plots into nine categories (Table 1), according to tree stand species composition and soil fertility, as described in detail in previous studies (Dyderski and Jagodziński 2020b). The division follows the phytosociological variability of invaded and uninvaded vegetation: Fagus sylvatica-dominated forest refers to Deschampsio flexuosae-Fagetum, an acidophilous beech forest; Quercus petraea-dominated forest refers to Calamagrostio arundinaceae-Quercetum petraeae, an acidophilous oak forest; and Quercus-Acer-Tilia forest refers to Galio sylvatici-Carpinetum, a fertile broadleaved forest. Pinus sylvestris-dominated forests represented two groups: poor (occupying mainly mesic sites of Leucobryo-Pinetum and Calamagrostio arundinaceae-Quercetum petraeae on podzols and brunic soils), and rich (on nutrient-rich luvisols and cambisols soils, which replaced Galio sylvatici-Carpinetum). In each of these two P. sylvestris groups we distinguished a variant invaded by P. serotina, which spontaneously colonized both types of forests. We assumed plots with more than 500 ind. ha⁻¹ of P. serotina individuals as invaded.

Table 1 Overview of soil characteristics, biomass, and species composition of forest types included in the study. Notation for soil pH, C:N ratio, shrub layer, and overstory biomass shows mean ± SE (min-max). Type with “-Ps” indicates Prunus serotina invasion

2.2 Data collection

We recorded all trees and shrubs with heights above 1.3 m during August 2014 (9 blocks) and 2015 (12 blocks) within each stand structure plot (Fig. 1). We measured all living individuals with a diameter at breast height (DBH) ≥ 5 cm including bark. For individuals with DBH < 5 cm we recorded the number of individuals. We identified all individuals at the species level in the field, following the Global Biodiversity Information Facility taxonomic backbone (GBIF 2019) for species nomenclature. For calculations of basal area (m² ha⁻¹) and biomass (see details below) we assumed the DBH of trees not measured to be 2.5 cm, as this is the mid-point of the non-measured interval (0–4.9 cm). Despite decreasing calculation accuracy in comparison with overstory trees reaching large diameters, we considered such bias would have low significance for the total results. We assumed the shrub layer as trees and shrubs reaching up to 1/2 of the height of the main canopy, classifying them in the field, during measurements. We used that classification to account for functional differences between forest strata rather than applying a particular threshold of DBH or height. During measurements we measured the slope of the plot using a clinometer. We also used data about soil type and soil C:N ratio from earlier studies (Dyderski and Jagodziński 2019b).

We calculated the aboveground biomass of each recorded individual using species- or genus-specific allometric models. When tree dimensions exceeded the maximum diameter of sample trees from the dataset used for a particular allometric model by > 20%, or specific models for particular species were not available, we used the general model for broadleaved trees (Forrester et al. 2017). For some species not reaching DBH > 5 cm we used species-specific models based on root collar diameter, assuming root collar diameter to be 2.5 cm. For details see the particular biomass models in the dataset (Dyderski and Jagodziński 2020a). We used the biomass of each species in each plot to calculate biomass proportions and for ordination analysis.

We analyzed three aspects of shrub layer alpha diversity—taxonomic, phylogenetic, and functional. We used species richness and Shannon’s diversity index as metrics of taxonomic alpha diversity. We calculated them using the vegan package (Oksanen et al. 2018). To assess phylogenetic diversity we used a phylogenetic megatree included in the V.phylo.maker package (Jin and Qian 2019) to construct a tree of species occurring in shrub layers of study plots. We calculated Faith’s phylogenetic diversity (i.e., sum of phylogenetic tree branch lengths, representing all species present in the community) and mean pairwise phylogenetic distance between species within the community, using the PhyloMeasures package (Tsirogiannis and Sandel 2016). We obtained functional traits from LEDA (Kleyer et al. 2008), BIEN (Enquist et al. 2016), BiolFlor (Klotz et al. 2002), and Pladias (Wild et al. 2019) databases: pollination mode, flowering start, and duration, specific leaf area, seed mass, height, and wood density. We obtained complete information about flowering, pollination, and seed mass traits, but SLA was available only for 88.9%, height—for 98.1%, and wood density—for 61.1% of species. We imputed missing data (see Pyšek et al. 2015 for rationale) using random forest imputation (Penone et al. 2014), implemented in the missForest package (Stekhoven and Bühlmann 2012). This method estimated missing values using known trait values and the first ten phylogenetic eigenvectors (Diniz-Filho et al. 1998), obtained using the PVR package (Santos 2018) and covering 71.7% of the variation in phylogenetic distances among species. Normalized root-mean-squared error of imputed traits was 0.0314. We calculated two indices of functional diversity: functional richness, expressing the quantity of plant functional types present in a community, and functional dispersion, expressing the size of the community species traits hypervolume within the functional trait space (Mason et al. 2005; Laliberté and Legendre 2010), using the FD package (Laliberté et al. 2014).

2.3 Data analysis

We analyzed data using R software (v. 3.5.3; R Core Team 2019). We assessed the species composition of shrub layers within forest types using non-metric multidimensional scaling (NMDS), implemented in the vegan package (Oksanen et al. 2018). We used species biomass within study plots as an abundance metric and study plots as analytic units (n = 168). Species biomass was standardized using log-transformation, and NMDS was based on Bray-Curtis distance matrix. We tested the significance of differences among forest types using adonis–permutation-based multivariate ANOVA, also accounting for a block design in permutations. We assessed differences among forest types studied using linear (LMMs) and generalized linear mixed-effects models (GLMMs), implemented in the lme4 package (Bates et al. 2015). We accounted for spatial dependence within blocks by including block ID as random intercepts in the models. In the initial models we included forest type, overstory aboveground biomass, stand age (data from management plans), slope, soil type, soil C:N ratio, and soil pH, and then we reduced models to decrease Akaike’s Information Criterion, corrected for small sample size (AICc). We also reported AICc₀–AICc of null models (intercept-only) and AICc_full–AICc of the full model (including all variables) to show how final models were improved. After inspections of histograms and model QQ plots we used a log-normal distribution of shrub layer aboveground biomass, and normal distributions of other alpha diversity metrics. Due to the discrete character of species richness we assumed a Poisson distribution, after ensuring that the model was not overdispersed. We inspected residual distributions, impacts of outliers on results, and dispersion using diagnostic tests implemented in the DHARMa package (Hartig 2020). We also reported marginal (R²_m) and conditional (R²_c) coefficients of determination, expressing the proportion of variance explained by fixed effects only, and both fixed and random effects, respectively (Nakagawa and Schielzeth 2013), implemented in the MuMIn package (Bartoń 2017). In the case of functional richness we excluded plots where this parameter was unavailable due to low species richness (less than three species). To assess marginal effects of forest type (i.e.. assuming mean values of other predictors and excluding random effects), we calculated marginal means using the emmeans package (Lenth 2019) and we checked differences among variants using Tukey posteriori tests with multiple hypothesis testing adjustments. While interpreting the results, we followed the guidelines of the American Statisticians Association (Wasserstein and Lazar 2016) and focused more on effect sizes than on p-values of pairwise comparisons. This allowed us to make conclusions based more on biological effects rather than on mere statistical measures undermined by small sample sizes.

3 Results

3.1 Shrub layer aboveground biomass

Despite testing stand age, slope, soil variables, and overstory aboveground biomass, the final model of aboveground biomass included only forest type (AICc = 409.1, AICc_full = 489.1, AICc₀ = 490.8, R²_m = 0.337, R²_c = 0.529; Table 2; Fig. 2). We found the lowest mean aboveground biomass of shrub layers in Q. rubra forests (0.52 ± 0.14 Mg ha⁻¹), which was about one-third of that in F. sylvatica forests (1.43 ± 0.43 Mg ha⁻¹) and one-fourth of that in Q. petraea forests (2.17 ± 0.54 Mg ha⁻¹). In both types of P. sylvestris forests we found approximately 2.5-times higher aboveground shrub biomass in P. serotina invaded than non-invaded forests (8.11 ± 1.87 vs. 3.50 ± 0.59 Mg ha⁻¹ in rich and 5.59 ± 1.23 vs 2.09 ± 0.78 Mg ha⁻¹ in poor P. sylvestris forests); however, these differences were statistically insignificant in a pairwise comparison. R. pseudoacacia and Quercus-Acer-Tilia forests (4.85 ± 1.41 and 4.57 ± 1.10 Mg ha⁻¹, respectively) had higher shrub layer biomass than Q. rubra and F. sylvatica forests. Their shrub layer biomass was lower than invaded rich and poor P. sylvestris forests, but the differences were not statistically significant.

Fig. 2

Mean (+ SE) shrub layer aboveground biomass estimated using GLMM assuming log-normal distribution. Letters denote variants that are not different at p = 0.05, according to Tukey posteriori test. For model details see Table 2

3.2 Shrub layer species composition

The shrub layers in forest types studied were dominated by the young regeneration of trees rather than shrub species (Table 3). Ordination (Fig. 3) revealed a gradient of species composition from fertile forest types on the left side of NMDS space (Robinia, Quercus-Acer-Tilia, and rich P. sylvestris forests) to nutrient-poor forest types on the right side (poor P. sylvestris and Q. petraea forests). Forest types differed significantly in species composition (ADONIS: d.f.  = 8, F = 10.01, p = 0.001) and explained 33.9% of the variation in species composition. Shrub layers dominated by shade-tolerant species, associated with high canopy cover, were grouped at the bottom of NDMS space. We also found that few species scores indicated their ecological optima in shrub layers outside forest types where they are dominants: R. pseudoacacia and Acer pseudoplatanus scores were shifted towards P. sylvestris forests (Fig. 4). Among invasive species studied P. serotina comprised a majority of biomass in most forest types. However, in non-invaded stands, its biomass was lower. We found only small proportions of Q. rubra in shrub layers in forest types other than those dominated by Q. rubra.

Fig. 3

Results of non-metric multidimensional scaling (NMDS) of shrub layer aboveground biomass proportions among species in each study plot (n = 168), colored by forest type. NMDS stress = 0.1638. Labels represent species scores occurring in at least five study plots (alien species marked by red labels), abbreviated to the first three letters of each part of a binomial name (e.g., Fag.syl = Fagus sylvatica)

Fig. 4

The proportions of main species biomass in shrub layers among distinguished forest types. Boxplots indicate the first and third quartiles, the line inside each box indicates the median, and whiskers indicate the minimum and maximum without outstanding observations (> 1.5 interquartile range). Points represent particular observations. For clarity we drew only species occurring in at least 20 study plots and observations from plots with proportions > 5%, sorted by overall proportion. Invasive species studied were marked by red color

3.3 Shrub layer alpha diversity

We found differences in alpha diversity among forest types studied for all biodiversity indices except Shannon’s diversity index and functional richness (Fig. 5; Table 4). All diversity metrics depended on the plot area; however, its effect size was low. Phylogenetic diversity indices were controlled by forest type functional richness depended on soil type and forest type. We found the lowest taxonomic diversity in F. sylvatica and Q. rubra forests, while the highest—in all P. sylvestris forest types, R. pseudoacacia, and Quercus-Acer-Tilia forests. We found the highest phylogenetic diversity and mean pairwise distance in non-invaded poor P. sylvestris forests, both types of rich P. sylvestris forest and Quercus-Acer-Tilia forests, while the lowest—in F. sylvatica forests. We found the lowest functional richness in Q. rubra and F. sylvatica forests, and the highest—in Quercus-Acer-Tilia, Q. petraea, and non-invaded poor P. sylvestris forests. The functional richness of Q. petraea forests was twice that of Q. rubra forests. We found the lowest functional dispersion in Fagus and invaded rich P. sylvestris forests, while the highest was in Quercus-Acer-Tilia forests.

Fig. 5

Mean (+ SE) values of shrub layer alpha biodiversity indices estimated using GLMMs. Letters denote variants which are not different at p = 0.05, according to Tukey posteriori test, n.s. the final model did not reveal an impact of forest type. For model details see Table 4

4 Discussion

We found differences in shrub layer biomass among forest types studied. In most cases we found higher biomass than the averages from Lithuania fertile (0.896 ± 0.029 Mg ha⁻¹) and very fertile (1.384 ± 0.046 Mg ha⁻¹) sites (Škėma et al. 2015) and S Poland (1.540 ± 0.011 Mg ha⁻¹; Orzeł 2015). Our results are in line with findings from Lithuania where more fertile forest types had higher shrub biomass (Škėma et al. 2015). This contradicts the findings of Orzeł (2015), who found higher shrub layer biomass in the less fertile forest types.

Our study revealed that P. serotina and R. pseudoacacia invaded forests had higher shrub layer biomass than native forest types, except Q. petraea and F. sylvatica forests. The latter types are known for either low soil fertility, unfavorable for numerous shrub species (Dyderski and Jagodziński 2020b), or limited light availability, similarly to Q. rubra (Jagodziński et al. 2018). Negative impacts of Q. rubra might be compared with other tree species with a high leaf area. For example, a removal experiment in Picea abies monocultures revealed that removing 50% of tree basal area resulted in twice higher shrub species richness (Heinrichs and Schmidt 2009). Also, alien Pinus nigra, transmitting more light beneath the canopy layer than native species, had twice the shrub layer cover (Mikulová et al. 2019).

This is in line with the findings of Ostrom (1983), who found almost 2.5-fold higher shrub biomass in light-transmitting Larix laricina forests than shade-casting Thuja occidentalis forests. The negative impacts of Q. rubra on both shrub layer biodiversity and biomass is in line with studies concerning herb layer biodiversity (Marozas et al. 2009; Chmura 2013; Dyderski and Jagodziński 2020b) and dominant species biomass (Woziwoda et al. 2019). Similarly, the mechanism of increased light availability might be connected with the positive impacts of R. pseudoacacia on shrub layer biomass. Moreover, due to the ability to absorb nitrogen from the atmosphere (Rice et al. 2004), R. pseudoacacia forests allow for higher growth rates of understory plants, both woody and non-woody. This can be also indicated by high overstory biomass in R. pseudoacacia forest (Table 1). Therefore, previous studies reported positive (e.g., Gentili et al. 2019), negative (e.g., Slabejová et al. 2019), or no effects (e.g., Hejda et al. 2017) of R. pseudoacacia on understory vegetation.

The main differences among invasive species studied are connected with their life forms—P. serotina usually dominates shrub layers, contributing to higher shrub layer biomass. It is connected with the ability for fast height increments after release from growth suppression (Closset-Kopp et al. 2007), allowing it to increase biomass up to 24,000 times in eight years (Jagodziński et al. 2019). Although R. pseudoacacia is also able to dominate the shrub layer (Slabejová et al. 2019), we did not find such an effect, due to limitations on the survival of natural regeneration (Dyderski and Jagodziński 2019b). We found lower effects of P. serotina on shrub layer biomass in rich P. sylvestris forests than in poor P. sylvestris forests. This may be connected with higher soil fertility, hosting more native species of shrubs, mainly typical of broadleaved forests (Zerbe and Wirth 2006). P. serotina successfully colonized both habitat types and outcompeted native species. Different impacts of invasive species in low and high soil fertility are widely known in invasion ecology as context-dependence (Sapsford et al. 2020). This highlights the need to account for habitat-specificity in management plans for invasive trees.

5 Conclusion

Our study demonstrated how invasive tree species influenced productivity and biodiversity in temperate forests. Depending on forest management and conservation aims, removal of invasive trees might lead to decreasing ecosystem biomass pools but allow for regeneration of native biodiversity. However, impacts are species- and context-dependent, therefore decision-making about the introduction or eradication of invasive tree species requires accounting for a wide range of impact assessments. Moreover, we also provided new data about the primary production and carbon sequestration of the shrub layer.

6 Disclaimer

The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Appendix

Table 2 Impacts of forest type on shrub layer aboveground biomass, estimated using LMM assuming log-normal distribution. AICc = 409.1, AICc_full = 489.1, AICc₀ = 490.8. Significant variables were italicized

Table 3 Overview of species frequency (%) within forest types, divided into shrubs and regeneration of tree species

Table 4 Models of shrub layer alpha diversity, estimated using Gaussian distribution GLMMs (for species richness – Poisson distribution and z statistic instead of t). Significant variables were bolded. AICc – Akaike Information Criterion, AICc₀ – AIC of null model (intercept-only), AICc_full – AIC of full model (including all hypothesized variables), R²_m – marginal coefficient of determination, R²_C – conditional coefficient of determination; *for Shannon’s index model block random effect SD = 0.43, AICc = 198.4, AICc_full = 275.0, AICc₀ = 199.5, R²_m = 0.005, R²_c = 0.577