The long‐term effects of invasive earthworms on plant community composition and diversity in a hardwood forest in northern Minnesota

Abstract Nonnative European earthworms are invading hardwood forests of the Chippewa National Forest, MN. While effects on plant communities at the leading edge of invasion have been studied, little is known about longer‐term effects of invasive earthworms. We applied a model using historic O‐horizon soil thickness and a chronosequence approach to classify 41 hardwood sites in the Chippewa National Forest as “long‐term wormed” (wormed >2 decades), “short‐term wormed” or “unwormed/lightly wormed.” Graminoids, especially Carex pensylvanica, had the greatest mean percent cover in sites that had been wormed for over two decades. The families with the greatest negative change in mean percent cover after over two decades of earthworm invasion were Asteraceae, Violaceae, and Sapindaceae (specifically Acer species). Across all diversity metrics measured, long‐term wormed sites had the lowest understory plant species diversity, short‐term wormed sites had intermediate diversity, and unwormed/lightly wormed sites exhibited the highest diversity. Long‐term wormed sites had the lowest mean species richness across all sample scales (1–1024 m2). The greatest within‐group compositional dissimilarity occurred at sites that had been wormed for over two decades, suggesting that sites that had been wormed for over two decades have not reached a compositionally similar end‐state “wormed” community type. Our study suggests that understory diversity will decrease as hardwood forest stands become wormed over time. While our results support other findings that exotic earthworm invasion is associated with lower understory plant diversity in hardwood forests, our study was the first to use space‐for‐time substitution to document the effects after multiple decades of earthworm invasion.


| INTRODUC TI ON
European earthworms can have undesirable impacts on hardwood forests of the Great Lakes region in North America that were once devoid of earthworms due to glaciation (James, 1995;Resner et al., 2015). The establishment of invasive earthworms in these forests has been associated with a wide range of repercussions for hardwood forest ecosystems, including changes to ecosystem function (Dobson et al., 2017), ecosystem services (Frelich et al., 2019), and biodiversity (Ferlian et al., 2018). Recent studies suggest that invasive earthworms can affect hardwood forest plant communities both indirectly and directly. Indirectly, earthworms can alter soil structure and composition (Knowles et al., 2016), disrupt microbial networks (Dempsey et al., 2013), alter nutrient availability (Ferlian et al., 2020), and affect hydrological cycles (Catford, 2017). Directly, earthworms can ingest seeds (Cassin & Kotanen, 2016), consume roots and leaves (Curry & Schmidt, 2007), and modify seed bank composition .
Although studies have assessed the initial effects of invasive earthworms on plant community composition across leading edges of earthworm invasion in hardwood forests of the Great Lakes region Holdsworth et al., 2007), we are not aware of any studies that have used permanent plots to assess the effects of invasive earthworms on plant community composition in hardwood forests over multiple decades. However, one challenge of using permanent plots to study the long-term effects of invasive earthworms is that study sites that have historic data on soils and plant communities generally lack information about historic earthworm presence. Here, we addressed this problem by creating a generalized linear model using the thickness of the O horizon to classify plots as wormed or unwormed when they were first sampled (Alexander et al., 2020).
Although studies commonly use earthworm abundance or biomass to approximate earthworm impacts (Mahon & Crist, 2019;McCay & Scull, 2019), these two measures typically increase throughout the summer months. Because our sampling took place over multiple months, earthworm abundance and biomass were not appropriate for inclusion in the GLM model described in Alexander et al., 2020. In this paper, we used permanent plots because it allowed us to classify plots by history of worm invasion. However, differences among observers meant that we could not make direct comparisons using historic vegetation data. Therefore, we used a space-for-time substitution approach to investigate differences among understory plant communities in plots that differed in the duration of worm invasion. To the best of our knowledge, our work is the first to analyze the effects of invasive earthworms on plant community composition and richness over multiple decades in the Great Lakes region. We sought to answer the following questions: (1) which plant species and families are the relative winners and losers as earthworms become established, and are there patterns in plant functional traits, such as mycorrhizal status or life form, that predict plant success?
(2) What are the effects of invasive earthworms on understory plant diversity, and do these effects vary with scale? (3) What are the effects of invasive earthworms on plant community composition? 2 | ME THODS

| Study area
The Chippewa National Forest, located in north-central Minnesota, contains permanent vegetation plots that were initially established and inventoried between 1989 and 1996 ( Figure 1). Although these plots have historic data on soils and plant communities, no earthworm data were collected when they were initially sampled. In recent decades, nonnative European earthworms have been invading hardwood forest soils in the Chippewa National Forest through their widespread use as fishing bait and through transportation on logging roads and via horticultural soils (Hendrix & Bohlen, 2002).
Sampling took place within the three counties encompassing the forest: Cass, Beltrami, and Itasca. The Chippewa National Forest has a humid continental climate, with warm summers and cold winters. The mean monthly minimum temperature in the nearby town of Cass Lake is −20.8°C in January and the mean monthly maximum temperature is 19.3°C in July (Arguez et al., 2010). The mean monthly precipitation in Cass Lake is 2.2 inches (Arguez et al., 2010).
The permanent plots selected for this study were limited to Mesic hardwood forests. Soils in this area are primarily well-drained, loamy haplic glossudalfs (Richardson, 1997).
Historic permanent vegetation plots (relevés) were first established and sampled from early June to early September from 1989 to 1996 as part of a national effort within the U.S. Forest Service to map terrestrial ecosystems at several different scales (Cleland et al., 1997). On the Chippewa National Forest, this effort was completed in cooperation with the Minnesota County Biological Survey of the Minnesota Department of Natural Resources (Aaseng et al., 2011). In this study, a 'relevé' refers to a sampling unit and procedure modified from that first introduced by Braun-Blanquet in the early 1900's (Braun-Blanquet, 1932). During the initial 1989 to 1996 sampling period, representative samples of mature, homogeneous vegetation were targeted for plant community classification and vegetation characterization studies.
A subset of the relevés established in the Chippewa National Forest between 1989 and 1996 was inventoried between mid-May and mid-August 2017. In the selection process, we sought to reduce potential sources of variation other than worm impacts. We used four criteria: (1) we sampled only Northern Mesic Hardwood Forest or Northern Rich Mesic Hardwood Forest, (2) all plots had to have been unaffected by recent management (plots harvested between the two sampling periods were excluded and the age of re-inventoried forest stands thus ranged from 50 to 137 years, with an average stand age of 90 years), (3) to distribute the plots throughout the study area, all plots were at a minimum of 10 chains (approx. 200 m) apart, and (4) surveyor, with priority given to relevés which were sampled by a surveyor with known soil science expertise. Within these criteria, we also sought out plots with a high likelihood of remaining unwormed because such sites have become increasingly rare. Thus, our study was not designed to assess the percent of relevés that were wormed or unwormed as an indicator of the overall rate of earthworm invasion or spread in the Chippewa National Forest. At one sampling location, a species-area plot was not completed and therefore our analysis consisted of 41 relevés and 40 species-area plots.

| Sampling methodology
The modified relevé method as applied on the Chippewa National Forest involved establishing a 10 × 10-m (100 m 2 ) plot for subsequent soil and vegetation sampling. Relevés were re-established by relocating magnets placed at the center of the relevé during the 1989 to 1996 sampling period using a metal detector. Methods for sampling soil horizons during the initial sampling period were replicated during the 2017 sampling season. Soil sampling involved digging a small soil pit (roughly 0.5 m 2 ) and sampling soils to a depth of 30 cm. Thickness of the A-horizon and О-horizon were recorded in centimeters during both sampling periods, as was the texture and color of the A-horizon. This soil sample was taken within a few meters outside of the permanent relevé in a random direction to avoid damaging vegetation inside the relevé. Inventorying plants in the 10 × 10-m relevés involved recording all vascular plants and estimating their abundance using the Braun-Blanquet cover/abundance scale ( Table 1). Although the Minnesota Department of Natural Resources relevé protocol includes information about vegetation at all height strata, we limit our analyses to understory vegetation less than 2 m.
In addition to inventorying permanent 10 × 10-m relevés, new 32 × 32-m nested species-area plots were established at permanent relevé locations in 2017 (n = 40, Figure 2). Species-area plots were established as part of a statewide effort by the Minnesota Department of Natural Resources to understand species richness at multiple scales. The direction of these species-area plots was expanded from the permanent 10 × 10-m relevé was generally random, although we attempted to select a direction where the community  Prior to analyses, some taxa for which field identification to species was particularly challenging were grouped together at the generic (e.g.,: Amelanchier) or subgeneric (e.g.,: Viola pubescens/ canadensis) level. All statistical analysis was conducted in R Studio (R Core Team, 2020). For all understory plant species or plant families that occurred in greater than five relevés, we created a table of mean percent cover for unwormed/lightly wormed, short-term wormed, and long-term wormed relevés. We used understory Küchler life form cover data from the 10 × 10-m relevés to assess whether life forms (woody evergreen, woody deciduous, forbs, and graminoid plant species) explained understory differences among the three worm groupings.

| Statistical analyses
Mean percent cover was used to assess differences in species, family, and Kuchler life forms among unwormed/lightly wormed, short-term wormed, and long-term wormed plant communities in 10 × 10-m relevés. The Kruskal-Wallis H test was used to test the significance of mean percent cover differences among earthworm categories.
Four biodiversity indices (Shannon-Wiener, inverse Simpson's index, Pielou's evenness, and species richness) were used to assess differences in diversity among unwormed/lightly wormed, short-term wormed, and long-term wormed plant communities in 10 × 10-m relevés. Shannon-Wiener diversity is more sensitive to both species richness and species evenness than Simpson's index.
Because Simpson's index decreases as diversity increases, we reported the inverse Simpson's index (1/D). Species richness was calculated as the number of understory plant species per relevé, or the number of total species per species-area plot. A one-way ANOVA with pairwise comparisons was conducted to test the significance of mean Shannon-Wiener diversity, Simpson's index, and species richness among unwormed/lightly wormed, short-term wormed, and long-term wormed relevés.
We created a species accumulation curve for each of our three earthworm categories: unwormed/lightly wormed, short-term wormed, and long-term wormed species-area plots. The curves describe the cumulative count of species per increase in area sampled for species-area plot data. We used a logarithmic linear model to plot the line of best fit with the "ggplot2" package in R (Wickham, 2011).
Nonmetric multidimensional scaling (NMS) was used to ordinate plant species and sites and to explore gradients in plant species composition in relation to earthworm presence in the Chippewa National Forest. NMS is an ordination tool with several unique properties, including being the only ordination technique where sample separation is directly linked to sample dissimilarity. NMS does not assume linearity, it does not presume an underlying model of species response, and it does not assume an inherent dimensionality of the data. In NMS, the user chooses a set number of axes for ordination prior to analysis, and the data are subsequently fitted to these dimensions (Kruskal 1964). The NMS ordination was created using understory quantitative vegetation data collected in 2017. We created an NMS ordination using the 2017 10 × 10-m relevé data with the "vegan" package in R (Oksanen et al., 2020). Prior to ordination, data were transformed using a Wisconsin double standardization and square root transformation. A Bray-Curtis dissimilarity matrix was used for NMS analysis. After selecting an appropriate dimensionality for the NMS analysis, a 30 iteration NMS ordination with random starts was run and the iteration with the lowest stress level was selected. The NMS configuration was rotated to orient axes that explained the most variance as the primary axes. Finally, sites were added as points on the NMS plot and categorized by worm grouping.
The "envfit" function in the "vegan" package was used to assess correlations between a subset of environmental variables, earthworm assemblages, and NMS axes. Compositional differences between unwormed/lightly wormed, short-term wormed, and long-term wormed plant communities were tested using the permutational multivariate analysis of variance method (perMANOVA) described in Anderson (2001). This was accomplished using the "adonis2" function of the "vegan" package in R (Oksanen et al. 2018). The perMANOVA was calculated using a Bray-Curtis dissimilarity matrix applied to the 2017 understory relevé and species-area plant data and 100,000 permutations. A post hoc pairwise comparison of the perMANOVA results was conducted with a false discovery rate (FDR) correction and 100,000 permutations using the "pairwise. adonis" function in R (Martinez Arbizu, 2017).

| RE SULTS
Locations of relevés resampled in 2017 are depicted in Figure 1. (endogeic), Dendrobaena octaedra (epigeic), Eiseniella tetraedra (epigeic), and Dendrodrilus rubidus (epigeic). We performed a one-way ANOVA test to ensure that observed differences in plant communities among unwormed/lightly wormed, short-term wormed, and longterm wormed relevés could not simply be ascribed to differences in canopy cover among the three worm categories. We found no statistically significant differences in canopy cover among the three worm categories [F(2, 38) = 1.07, p = 0.352].

| Winners and losers
Woody deciduous species had similar mean percent cover for unwormed/lightly wormed, short-term wormed, and long-term wormed relevés (Appendix 1). Woody evergreen species had low mean percent cover at all relevés, with the greatest mean cover at unwormed/ lightly wormed sites and the smallest mean percent cover at longterm wormed sites. The largest differences in mean percent cover occurred with graminoid species, with over 30% lower mean percent cover of graminoids in unwormed/lightly wormed (11.7%) compared to long-term wormed (45.0%) relevés. Forbs had the highest mean percent cover in short-term wormed relevés, and the lowest mean percent cover in long-term wormed relevés.
The families Cyperaceae and Sapindaceae occurred in all sampled relevés ( Table 2). The greatest relative change in mean percent cover among families occurred for the Cyperaceae family, which had a mean cover of 7.4% in unwormed/lightly wormed relevés and a mean cover of 40.0% in long-term wormed relevés.

| Diversity
Across all measures of diversity (Shannon-Wiener, inverse Simpson, Pielou's evenness, and species richness), unwormed/lightly wormed sites had the highest diversity values and long-term wormed sites experienced the lowest diversity values (Figure 3). A oneway ANOVA test indicated statistically significant differences in

| Species richness and scale
Species accumulation curves indicated that long-term wormed plots had the lowest mean richness at all measured scales (1-1024 m 2 , Figure 4). Short-term wormed plots had similar low mean richness as long-term wormed plots at the smallest scale, but at the largest scale measured (1024 m 2 ), short-term wormed plots had mean richness similar to unwormed/lightly wormed plots (Figure 4). At 1024 m 2 , unwormed/lightly wormed sites had the greatest mean total number of species (60), when compared to short-term wormed (58), and longterm wormed (50).

| Plant community composition
An optimal NMS solution with two dimensions and a stress level of 0.17 was selected for the 2017 relevé understory vegetation data ( Figure 5). Environmental variables that were significantly correlated with NMS axes include Easting, Northing, elevation, 1989Northing, elevation, -1996 A-horizon thickness, and 2017 A-horizon thickness ( Table 3).
Earthworm-related variables that were significantly correlated with NMS axes include anecic earthworm presence, epi-endogeic/anecic earthworm presence (juvenile Lumbricus species), and sum of earthworm eco-groups. As sites became more wormed, they appeared to cluster toward the right side of the NMS ordination, which was significantly associated with Northing, 1989-1996 A-horizon thickness, 2017 A-horizon thickness, anecic earthworm presence, epi-endogeic/ anecic earthworm presence, and sum of earthworm eco-groups.
TA B L E 2 Percent cover of understory plants by family and worm category for unwormed/lightly wormed, short-term wormed, and longterm wormed relevés. The "spp" column represents the number of species (or species complex) per family and 100,000 permutations revealed significant compositional differences between short-term and long-term wormed sites (Table 4), as well as between unwormed/lightly wormed and long-term wormed sites (Table 4). Compositional differences between unwormed/ lightly wormed and short-term wormed sites were not statistically significant (Table 4).

| DISCUSS ION
Taken together, these results indicate that earthworm invasion is associated with reduced species diversity and altered composition in forest understories and, further, that worming effects increase with the duration of earthworm invasion such that even sites wormed for more than two decades are still experiencing species loss and composition shifts.

| Winners and losers
Differences in abundance of understory plant species after earthworm invasion followed similar distributional patterns as was reported by Hale et al. (2006). Of the relevés sampled in 2017, Carex pensylvanica had a sevenfold difference in cover among the wormed categories. Unlike the majority of native understory plants that are present in sugar maple--basswood forests of the Great Lakes region, Carex pensylvanica is non-mycorrhizal (Brundrett & Kendrick, 1988). Nonmycorrhizal plant species may receive a relative competitive advantage after earthworms invade and disrupt existing mycorrhizal networks.
Arisaema triphyllumis a native forb that possesses calcium oxalate crystals in its foliage and roots. These crystals may deter above and belowground herbivores, including earthworms. Hale et al. (2006)  and Acer spicatum. Previous studies found that earthworm invasion can have negative effects on Acer saccharum colonization (Cassin & Kotanen, 2016;Lawrence et al., 2003), although mean percent cover of understory Acer saccharum was greater in long-term wormed relevés compared to unwormed/lightly wormed relevés in this study. Both Acer saccharum and Acer rubrum are known to form symbioses with arbuscular mycorrhizae (Brundrett and Kendrick 1988;Wiseman & Wells, 2005). Overall, the Acer genus had slightly greater mean cover in short-term wormed relevés than unwormed/ lightly wormed relevés, but in long-term wormed relevés the Acer genus had mean understory cover values less than unwormed/ lightly wormed relevés, indicating members of the Acer genus are "losers" as earthworms become established over multiple decades.
Mean total cover of the Araliaceae family was lowest in relevés that had been wormed for over two decades. Aralia nudicaulis experienced a 45.9% decrease in mean cover in long-term worm relevés (2.0% cover) when compared to unwormed/lightly wormed releves (3.7%).
Both Aralia nudicaulis and Aralia racemosa are obligately mycorrhizal (Lawrence et al., 2003). In a mesocosm experiment by Hale et al. (2008), Aralia racemosa did not experience significantly higher mortality when exposed to earthworms over a period of 13-18 weeks, perhaps echoing our results that decreases in Aralia racemosa cover associated with invasive earthworms must be observed over a longer time frame.
The effects of invasive earthworms on mean cover of understory vegetation varied by Küchler life form (woody deciduous, woody evergreen, graminoids, and forbs). Mean woody deciduous cover was lowest in unwormed/lightly wormed sites, followed by long-term wormed sites. Mean woody deciduous cover was highest in shortterm wormed sites. This phenomenon may be a result of a flush of readily available nutrients available for plant uptake immediately after earthworm invasion, followed by a longer-term decrease in nutrients as these nutrients become prone to leaching and outwash (Resner et al. 2015). Woody evergreen species had relatively low dominance in all categories, likely because relevés selected for resampling were classified as sugar-maple basswood forests. Mean cover of forbs was intermediate at unwormed/lightly wormed sites, greatest at short-term wormed sites, and reached its lowest value at longterm wormed sites. Forbs may experience an initial flush of growth following a recent earthworm invasion, followed by a longer-term decrease in total cover (Hale et al., 2008). Like the growth observed in short-term wormed understory woody deciduous vegetation, the initial pulse of growth observed in short-term wormed forbs could be due to a buildup of readily available nutrients resulting from mineralization of the O-horizon (Dobson et al., 2017). The two plant families with the largest total mean negative change in cover from unwormed/lightly wormed to long-term wormed relevés were both exclusively represented by forbs in our dataset: Asteraceae (43% decrease in mean cover) and Violaceae (48% decrease in mean cover).
Mean total cover of graminoids experienced the greatest degree of change as earthworms became established over time. The greater mean cover values of graminoids that we observed after short-term and long-term earthworm invasion mirrors that reported by Craven Cyperaceae (+32% mean cover, a 440% increase), but still showed a net positive change in total mean cover (+1% mean cover, a 28% increase) in unwormed/lightly wormed relevés compared to relevés that had been wormed for over two decades. Graminoid species that possess the following traits may experience a competitive advantage after earthworms become established: tolerance of drought and root herbivory, the presence of bud banks, and non-obligate mycorrhizal associations (Craven et al., 2017).

| Diversity
Similar to findings by Hale et al. (2006), Holdsworth et al. (2007, and Hopfensperger et al. (2011), our study found that understory plant diversity was lowest in sites that had experienced earthworm invasion. Across all diversity metrics measured (Shannon-Wiener, inverse Simpson's, Pielou's evenness, and species richness), longterm wormed relevés had lower understory plant species diversity than unwormed/lightly wormed and short-term wormed relevés.
Unwormed/lightly wormed relevés exhibited the highest average diversity (Shannon-Wiener, inverse Simpson's, Pielou's evenness, and species richness) at the 10 × 10-m scale. Species evenness was lowest in sites that had been wormed for over two decades, likely reflecting a greater dominance of Carex pensylvanica monocultures.

| Species richness and scale
We were not able to find published information about the effects of invasive earthworms on understory plant species richness at multiple scales. In our study, we found that at the smallest scale (1 m 2 ) unwormed/lightly wormed species-area plots had the greatest mean species richness. At this scale, short-term and long-term wormed sites had similar numbers of mean species richness. At intermediate scales, unwormed/lightly wormed species-area plots had the highest mean species richness, followed by short-term wormed sites, and long-term wormed sites had the lowest mean species richness. At the largest scale, unwormed/lightly wormed and short-term sites had similar mean species richness, although unwormed/lightly wormed mean species richness was highest. The lower mean species richness in long-term wormed sites became more pronounced at the largest scale (1024 m 2 ). Our results indicate that sites that have been wormed for multiple decades have, on average, lower species richness at both small (1 m 2 ) and large (1024 m 2 ) scales. The full extent of reduced species richness associated with earthworm invasion may not be apparent until earthworms have been established for multiple decades.
F I G U R E 5 NMS plot for unwormed/ lightly wormed, short-term wormed, and long-term wormed plant communities from 10 × 10-m relevé data collected in hardwood forests of the Chippewa National Forest in 2017. Larger circles represent the centroid for each community type. Smaller circles represent individual species-area plots. Stress for the NMS ordination was 0.17

| Plant community composition
Sites that were classified as unwormed/lightly wormed (withingroup) had the lowest dissimilarity, indicating these sites were the most compositionally similar. The second least compositionally dissimilar categories were unwormed/lightly wormed and short-term wormed sites (between-group), indicating that these sites also shared a large degree of compositional similarity. Sites that had been wormed for the longest period had the greatest rates of dissimilarity between short-term wormed (between-group) and unwormed/lightly wormed (between-group) relevés. This may indicate that the greatest compositional changes in plant communities that have been invaded by European earthworms occur after two decades or more of invasion. If greater within-group compositional dissimilarity indicates a community currently in transition, then it appears that both short-term wormed and long-term wormed sites are currently in a state of transition. Compositional dissimilarity was greatest at sites that had been wormed for over two decades, and sites that had been wormed for at least two decades did not appear to have reached a compositionally similar endstate "wormed" community type.
Anecic earthworms feed on both leaf litter and mineral soil and create extensive vertical burrows in the soil. The general sequence of earthworm invasion in the northern Midwest appears to follow the order of Dendrobaena (epigeic) > Lumbricus juveniles (anecic/epiendogeic) > Aporrectodea (endogeic) > Lumbricus terrestris (anecic) (Holdsworth et al., 2007). Epigeic earthworms invade first, typically followed by anecic, endogeic, and epi-endogeic earthworm species, subsequently followed by a high biomass of mature anecic earthworms (Loss et al., 2013). Therefore, sites with anecic earthworms present (located in the bottom right corner of the ordination diagram) are more characteristic of late-stage earthworm invasion. unwormed or had only leaf-litter dwelling earthworms present in 1/3 samples. It is essential that baseline data on plant community composition of unwormed sites is collected before these areas disappear, especially given their apparent diversity when compared to sites invaded by earthworms. In the Chippewa National Forest, unwormed sites may persist in isolated areas that function as habitat islands, such as forest stands surrounded by peatlands. These few remaining unwormed areas should be targeted for baseline studies.
If results of this study are indicative of future trends as earthworms become established, it can be expected that understory diversity will decrease as hardwood forest stands become wormed over time. There may be an initial flush of understory diversity immediately post earthworm invasion as nutrients become readily available, however this will likely be followed by a long-term de- The Great Lakes Worm Watch provided beneficial training in earthworm sampling and identification.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

CO D E AVA I L A B I LIT Y
Raw data and code are available online at github.com/jinny alexa nder/Earth worm_plant_data

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in GitHub at https://github.com/jinnyalexander/ Earthworm_plant_data.