Contrasting protist communities (Cercozoa: Rhizaria) in pristine and earthworm-invaded North American deciduous forests

Earthworms are considered ecosystem engineers due to their fundamental impact on soil structure, soil processes and on other soil biota. An invasion of non-native earthworm species has altered soils of North America since European settlement, a process currently expanding into still earthworm-free forest ecosystems due to continuous spread and increasing soil temperatures owing to climate change. Although earthworms are known to modify soil microbial diversity and activity, it is as yet unclear how eukaryote consumers in soil microbial food webs will be affected. Here, we investigated how earthworm invasion affects the diversity of Cercozoa, one of the most dominant protist taxa in soils. Although the composition of the native cercozoan community clearly shifted in response to earthworm invasion, the communities of the different forests showed distinct responses. We identified 39 operational taxonomic units (OTUs) exclusively indicating earthworm invasion, hinting at an earthworm-associated community of Cercozoa. In particular, Woronina pythii, a hyper-parasite of plant-parasitic Oomycota in American forests, increased strongly in the presence of invasive earthworms, indicating an influence of invasive earthworms on oomycete communities and potentially on forest health, which requires further research.


Introduction
Large parts of North American temperate forests have been devoid of native earthworm species since the last glacial period (Hendrix and Bohlen 2002;Bohlen et al. 2004;Callaham et al. 2006;Frelich et al. 2019). Although a variety of non-native earthworm species have been repeatedly introduced since the arrival of European settlers, extensive areas of the northern regions still remained earthworm-free due to their severe climatic conditions with harsh winters (Callaham et al. 2006). The increase of winter soil temperatures in conjunction with the current climate change led to an invasion of native forest ecosystems at higher altitudes and higher latitudes by non-native earthworm species (Fisichelli et al. 2013). Moreover, manifold human activities facilitate the distribution and spread of non-native earthworms in northern North America (Holdsworth et al. 2007;Cameron et al. 2007).
Earthworms are known as ecosystem engineers which profoundly shape soil ecosystems (Knollenberg et al. 1985;Lavelle et al. 1997;Butenschön et al. 2009;Eisenhauer 2010). Earthworms process large amounts of plant litter by fragmentation and bioturbation, with major consequences for the cycling of carbon and nutrients in forest ecosystems (Scheu 1987(Scheu , 1994Scheu and Wolters 1991;Fahey et al. 2013). Earthworm burrow walls and casts, i.e. the drilosphere, are soil hotspots with increased availability of C and N (Brown 1995;Eisenhauer 2010) and with much higher microbial activity and abundance than surrounding soil (Maraun et al. 1999;Tiunov et al. 2001;Tao et al. 2009;Andriuzzi et al. 2016). Subsequent shifts of soil microbial food webs lead to an altered soil animal and microbial community and may even affect plant community composition (Eisenhauer et al. 2007;Nuzzo et al. 2009;Craven et al. 2017;Ferlian et al. 2018). However, despite evidence for significant shifts in soil biodiversity and community composition in response to earthworm invasion (Ferlian et al. 2018), there have been no studies on the fate of soil protists.
There are indications of a tight trophic relationship of earthworms with protists (i.e. unicellular eukaryotes) (Bonkowski and Schaefer 1997;Scheu and Setälä 2002). Thus, indirect effects of earthworms, via casting and burrowing, and direct impacts through consumption of protist communities are expected with potential ripple effects on other trophic levels . Although soil protists have been predominantly considered as consumers of bacteria (Clarholm 1981(Clarholm , 1985Glücksman et al. 2010;Trap et al. 2016;Flues et al. 2017), they are increasingly recognised as consumers of fungi in soil ecosystems (Geisen et al. 2015a;Dumack et al. 2016Dumack et al. , 2018. Due to their high abundance, high turnover and assimilation rates, they occupy a key position for the energy transfer from microorganisms to higher trophic levels in forest soil food webs (Schaefer and Schauermann 1990;Schröter et al. 2003;Moore et al. 2005). Moreover, there is accumulating evidence for the significant role of plant parasitic protists to play crucial roles in many ecosystems (Delgado-Baquerizo et al. 2020;Hittorf et al. 2020).
Even with the rise of molecular methods in the assessment of microbial communities, the molecular assessment of protists remains difficult. Protists are polyphyletic and highly diverse; accordingly, there is no universal protocol to assess all protists. Primerbased metabarcoding is inevitably selective, as there is no general primer that enables the assessment of all protist groups (Lentendu et al. 2014;Hirakata et al. 2019). In addition, protists vary in their genome sizes and ribosomal copy numbers, which obscures results of metabarcoding surveys (Amend et al. 2010;Weber and Pawlowski 2013;Lentendu et al. 2014;Sibbald and Archibald 2017;Hirakata et al. 2019). However, despite these difficulties, there is a consensus that the protistan taxa Cercozoa and Amoebozoa represent the dominant protistan taxa in soil (Urich et al. 2008;Geisen et al. 2015b;Grossmann et al. 2016;Bonkowski et al. 2019) and forest litter (Voss et al. 2019). Specific metabarcoding protocols for the comprehensive assessment of the total diversity of one of these two taxa, the Cercozoa, have been developed (Bugge Harder et al. 2016;Fiore-Donno et al. 2018Degrune et al. 2019). The Cercozoa comprise a vast variety of different morphotypes, including flagellates, amoeboflagellates, naked amoebae, and testate flagellates as well as amoebae with diverse functional roles, such as bacterivores, fungivores, eukaryvores, and omnivores or algivores (Seppey et al. 2017;Bonkowski et al. 2019;Dumack et al. 2019). This means that protists do not only have many important ecosystem effects, they can also serve as crucial bioindicators for soil community composition and health (Geisen et al. 2020).
Following the protocols published by Fiore-Donno et al. (2018), we assessed the effects of earthworm invasion on cercozoan communities of three North American boreal and temperate-boreal forests (in Minnesota, USA, and Alberta, Canada). These sites have been invaded by anecic (Lumbricus terrestris), endogeic (Octolasion tyrtaeum) and epigeic (Dendrobaena octaedra) earthworms (Frelich et al. 2006;Eisenhauer et al. 2007). We hypothesised that (1) earthworms affect the composition of cercozoan communities in soil, and (2) that bioturbation by earthworms homogenises cercozoan communities across the upper soil profile.

Study site
The study sites were chosen as they have been under investigation for two decades and accordingly show a gradient of differing abiotic or biotic conditions along the earthworm invasion history (Eisenhauer et al. 2019;Jochum et al. 2021;Supplementary Table 1).
The first study site is located in a boreal forest north of Barrier Lake in Kananaskis Valley in the Rocky Mountains, south west Alberta, Canada (hereafter called Barrier Lake North; 51°02 0 N, 115°03 0 W). The site is at 1510 m a.s.l. on a slope facing south. The second study site is located south of Barrier Lake in Kananaskis Valley (hereafter called Barrier Lake South, 51°01 0 N, 115°04 0 W) at 1380 m a.s.l. on a slope facing west about 2.8 km from the first site. Both sites are characterised by a continental climate with cold, long winters and dry, short summers. The soil is usually frozen between November and March. The mean annual precipitation and temperature are 638 mm and 3.4°C, respectively. The soil type is classified as Orthic Grey Luvisol. The soil of the first site has a pH of 6.1 and an organic layer thickness of about 7.3 cm. The soil of the second site has a pH of 5.8 and an organic layer thickness of around 5.5 cm. The tree communities of both sites are dominated by Populus balsamifera L. and Populus tremuloides Michx. interspersed with single Picea glauca (Moench) Voss. The understorey vegetation is dense, species-rich and dominated by Rosa L. and Aster L. species, Viola canadensis L., Epilobium angustifolium L. Holub, and Delphinium glaucum S.Watson. The sites have moving earthworm invasion fronts as documented before (e.g. Eisenhauer et al., 2007;Straube et al., 2009) and confirmed by earthworm sampling conducted on all plots of the present study. The part of each site uphill from the invasion front was designated as uninvaded area with no or only low biomass (single individuals) of epigeic, endogeic and anecic earthworm species (hereafter called uninvaded area), whereas the part of each site downhill from the invasion front was designated as invaded area, where biomasses of epigeic, endogeic, and anecic earthworm species were high (hereafter called invaded area).
The third study site is located in a temperate-boreal forest, part of the Saint John's University, 120 km north west of Minneapolis, Minnesota, USA (hereafter called Saint John's, 45°34 0 N, 94°23 0 W). The site is flat at 350 m a.s.l. It is characterised by a humid continental, cold temperate climate. The mean annual precipitation and temperature are 701 mm and 6.2°C, respectively. The soil type is classified as Glossic Eutroboralf. The soil has a pH of 4.9 and an organic layer thickness of about 1.6 cm. The tree community of the site is dominated by Acer saccharum Marshall. The understorey vegetation is sparse and dominated by Urtica dioica L., Carex pennsylvanica Lam., Amphicarpaea bracteata (L.) Fernald, and Eurybia macrophylla (L.) Cass. The uninvaded and invaded areas of the site are around 400 m apart separated by a small lake. Areas of earthworm presence and absence/ low biomasses were confirmed by earthworm sampling conducted on the plots of the present study.

Sampling
In August 2016, we set up 20 plots of 1 9 1 m per forest site with an additional adjacent subplot of 0.3 9 1 m, ten in the uninvaded, ten in the invaded area, respectively. The plots were chosen randomly with a minimum distance between plots of 15 m. We distributed ten plots in the uninvaded and ten plots in the invaded area without any pairing of plots. The distance between the closest invaded and uninvaded plots was around 400 m in Alberta, whereas the distance between the closest sites was roughly 3 km. In Minnesota, the distance between uninvaded and invaded areas was roughly 2.5 km. Each plot was separated into four quarters which were assigned to different sample types, such as vegetation, plant litter, soil invertebrate macrofauna, and bulk soil. Soil samples for DNA sequencing analysis were taken at a separate adjacent subplot outside the four quarters to avoid disturbance during sampling. Four soil cores with a diameter of 2 cm were taken per plot. All soil cores were subdivided into 0-5 and 5-10 cm soil layers, homogenised per plot and soil layer using a 2-mm sieve, and transferred into 15 ml tubes. Samples were stored on dry ice in the field and for transport and later at -80°C in the lab until further processing. All tools were washed and sterilised with ethanol before and between sampling.

DNA Extraction
Soil microbial genomic DNA was extracted from each of the soil samples using a PowerSoil DNA Isolation Kit (MO BIO Laboratories Inc., Carlsbad, CA, United States) according to the slightly modified manufacturer's instructions. We used 0.4 g instead of 0.25 g of soil for the extractions and eluted in 100 ll of the elution buffer. DNA yields from each sample were checked with a NanoDrop ND-8000 spectrophotometer (Thermo Fisher Scientific, Dreieich, Germany), and the extracts were stored at -20°C for further analysis. The extracted DNAs were diluted by the factor of ten for subsequent PCRs.
PCRs and sequence processing PCRs were conducted in two steps. In the first PCR, the forward primers S616F_Cerco and S616F_Eocer were mixed in the proportions of 80% and 20%, and used with the reverse primer S963R_Cerco (Fiore-Donno et al. 2018). One ll of ten times diluted DNA were used as a template for the first PCR and 1 ll of the resulting amplicons were used as a template for a following semi-nested PCR. We employed the following final concentrations: Dream Taq polymerase (Thermo Fisher Scientific, Dreieich, Germany) 0.01 units, Thermo Scientific Dream Taq Green Buffer, dNTPs 0.2 mM and primers 1 lM. The conditions were set to an initial denaturation step at 95°C for 2 min, 24 cycles at 95°C for 30 s, 50°C for 30 s, 72°C for 30 s; and a final elongation step at 72°C for 5 min. The second PCR was conducted with barcoded primers (see Fiore-Donno et al. 2018, combinations used are given in Supplementary Table 2). All PCRs were conducted twice to reduce the possible artificial dominance of few amplicons by PCR competition, and then pooled.
A mock community with known species richness of diverse cultivated cercozoan taxa was run in parallel to assist the fine-tuning of the bioinformatics pipeline as described in Fiore-Donno et al. (2018). The amplicons were checked by electrophoresis and 25 ll of each pooled PCR product were purified and normalised using SequalPrep Normalization Plate Kit (Invitrogen GmbH, Karlsruhe, Germany). We then pooled the samples and the mock community and proceeded for a single library preparation. Library preparation and paired-end MiSeq sequencing with the MiSeqv3 2 9 300 bp kit were carried out by the Cologne Center for Genomics (CCG).
Paired reads were assembled using MOTHUR v.3.9 (Schloss et al. 2009) allowing one difference in the primers, no difference in the barcodes, no ambiguities, no mismatches higher than three and removing assembled sequences with an overlap \ 200 bp. Reads were sorted into samples according to the barcodes. The mock community functioned as an internal standard for the analyses and saturation was reached. We could retrieve the expected ten OTUs of the mock community at 97% similarity when deleting OTUs represented by\ 1000 sequences, which equals a relative abundance of less than 3.5%, and set this as a cutoff threshold for OTU delineation (see Fiore-Donno et al. 2018). The quality check and removal/cutting of low-quality reads were conducted with the default parameters. Using BLAST ? with an e-value of 1e -50 and keeping only the best hit, sequences were identified a self-updated version of the PR2 database and non-cercozoan sequences were removed (Guillou et al. 2013). Chimeras were identified using UCHIME as implemented in MOTHUR with a penalty for opening gaps of -5 and a template for aligning operational taxonomic units. Sequences were clustered using VSEARCH v.1, with abundancebased greedy clustering (agc) and a similarity threshold of 97% as indicated by analysing the mock community. Sequencing data of the protists were deposited at the European Nucleotide Archive under the accession number ERS2039495 (SAMEA104421553).

Statistical analysis
All statistical analyses and data visualizations were conducted in R version 3.5.2 with the vegan package (Oksanen et al. 2014). Based on the obtained sequence reads for each OTU, we calculated the Chao1 index as a measure of OTU richness, evenness, Shannon index as a measure of alpha diversity, and Bray-Curtis dissimilarities for beta diversity. Nested Permutational Multivariate Analysis of Variance (Permanova) was employed to test for differences in Cercozoa community composition across treatments. The herein-defined treatments were thus (A) the earthworm invasion, which was treated as a factor with two levels and the strata (B) sampling site (forest) and (C) soil depth. The site effect explained most of the variation in the dataset. To cope with the large influence of the site, we computed a partial dbRDA with the capscale function to partial out the site effect. An indicative value analysis (Dufrêne and Legendre 1997) was performed with the indicspecies package (De Cáceres and Legendre 2009) to identify indicator taxa in the respective treatments. The taxonomic network was made with the igraph package and based on the taxonomic assignment inferred from the PR2 database and the therein used hierarchy (Csardi and Nepusz 2006;Guillou et al. 2013). The code for the statistical analysis is available online on Github (https://github.com/Kenneth-Dumack/Ecoworm).

Results
We obtained in total 555 cercozoan OTUs from initially 8,748,548 sequences of 120 samples (3 forests 9 2 invasion status 9 2 soil depths 9 10 replicates; Fig. 1; Supplementary Table 3). Overall, the bacterivorous Sarcomonadea were the most diverse and most abundant taxon in this study. They represented about 45% of all OTUs (Fig. 1), with five OTUs ranking among the ten most frequently sequenced OTUs. The second most abundant OTU was Rhogostoma cylindrica, an omnivorous amoeba species. Moreover, we found a notably high number of sequence reads assigned to the bacterivorous Thaumatomonas (third most abundant OTU) and Spongomonadidae (eighth-most abundant OTU), and the oomycete hyperparasite Woronina pythii (ninth most abundant OTU).
Among the three forest sites, differences in cercozoan communities in terms of evenness, alpha and beta diversity could be detected (Permanova, richness: df [2,108] , F = 6.46, p = 0.1343, evenness: df [2,108] , F = 5.523, p = 0.0012, alpha diversity: df [2,108] , F = 1.95, p = 0.0007, beta diversity: df [2,108] , F = 12.65, p = 0.0001; Fig. 2). However, the differences in the three measures were mostly due to the distinctness of Saint John's in comparison to Barrier Lake North and South (Fig. 2). In contrast to both Barrier Lake sites, Saint John's exhibited a high difference in alpha diversity and evenness between soil layers (Fig. 2). In general, the the sampling site affected the diversity indices (alpha diversity, evenness, and beta diversity, see Supplementary Table 4). Although all three sampling sites accommodated the same taxa with no difference in richness (df [2,108] , F = 6.46, p = 0.13), the taxa responded differently to changes in soil depth and earthworm invasion, which is shown by the numerous significant interactions between these to treatments. For instance, there was a significant interaction between site and soil depth for alpha and beta diversity as well as evenness (Permanova, richness: df [2,108] , F = 1.21, p = 0.2952, evenness: df [2,108] , F = 5.09, p = 0.0078, alpha diversity: df [2,108] , F = 5.31, p = 0.0044, beta diversity: df [2,108] , F = 1.91, p = 0.0002; Supplementary Table 4).
The comparison of the influence of the treatments on cercozoan OTUs showed that in total 124 taxa were indicative for earthworm invaded soils, whereas 76 taxa indicated uninvaded soils. However, when taking into account the other treatments (site and depth), there were 39 taxa that exclusively indicated earthworm invasion (Supplementary Table 4). Interestingly, the abundance of Woronina pythii was dependent on earthworm invasion, showing a significantly higher relative abundance, with a mean increase about 20 times, in earthworm-invaded sites ( Supplementary Fig. 3).

Discussion
Earthworms affect soil microbial biomass and protist abundance directly via feeding (Bonkowski and Schaefer 1997) and indirectly via bioturbation, casting and burrowing (Maraun et al. 1999;Tiunov et al. 2001;Bonkowski et al. 2001;Tao et al. 2009). More recently, Andriuzzi et al. (2016) reported an altered community composition of protists in earthworm burrows. We investigated and compared how earthworm invasion affects cercozoan communities in northern North American forest soils.

Cercozoan communities differ among forests
The cercozoan communities of the three investigated sites were vastly different, and accordingly, cercozoan communities responded also differently to earthworm invasion. Differences in the cercozoan communities of our sampling sites, and thus also in response to earthworm invasion, may be explained by different biotic and abiotic conditions. Barrier Lake North and Fig. 1 Taxonomy-based network-like presentation of the diversity of Cercozoa in North American forest soils. Terminal nodes correspond to OTUs, and their respective size corresponds to the relative abundance of sequence reads in the dataset. Edges represent the taxonomic ranks from order to OTU level. The ten most abundant OTUs are labelled in grey according to their closest hit in the PR2 database South were dominated by poplar trees and shared a similar soil type, while soil type and tree species differed profoundly in Saint John's maple forest. Fiore-Donno et al. (2020) could assign differences in cercozoan communities to both, differences between forest tree species and soil types, but could not disentangle both factors, because specific forest tree species occurred on distinct soils types. Accordingly, the differences in tree species and soil type of Barrier Lake North and South and Saint John's most likely explain the strong differences between cercozoan communities between the sampling sites. Moreover, the two Barrier Lake forests and the Saint John's forest also considerably differed in other soil abiotic properties. The former two had comparably higher soil water content, carbon and nitrogen concentrations, and thickness of the humus layer (Jochum et al. 2021). The understory vegetation was more developed and comprised different species in the Barrier Lake forests compared to Saint John's forest. And in addition, the forests may have differed in their invasion stage and, thus, in the proportions of earthworm ecological groups that are known to differently affect soil biota (Ferlian et al. 2018). Biotic and abiotic factors were shown to have strong effects on the community composition of protists (Tedersoo et al. 2016;de Araujo et al. 2018) and likely contribute to patterns of spatial separation in protist communities over large geographic distances (Venter et al. 2017;Singer et al. 2018Singer et al. , 2019. These general differences in cercozoan communities across forests likely explain their dissimilar responses to earthworm invasion.

Cercozoan communities are not uniformly altered by earthworm invasion
The cercozoan community composition changed with earthworm invasion. Although the present species were similar in all three forest sites the changes induced or accompanied by earthworm invasion were  Fig. 2 Chao1 richness, Shannon alpha diversity, and Sheldon evenness of cercozoan OTUs in earthworm invaded forests in Northern American are shown not similar, which was shown by the numerous interaction effects in our statistical analysis. A reasonable interpretation for the observed patterns is therefore that earthworms and an interplay of environmental and potentially other biotic interactions affect the cercozoan community. Even the quite similar sampling sites, like Barrier Lake North and Barrier Lake South, showed contrasting results. For instance, Barrier Lake South shows no changes in alpha diversity and evenness, although Barrier Lake North does show a decrease in both (Fig. 2). We conclude that the complexity of how soil microbial communities are affected needs further investigation taking other groups of soil microorganisms into account.
Earthworm invasion alters protist communities and potential consequences for plant communities Sarcomonadea and an OTU assigned to Rhogostoma cylindrica (Thecofilosea) generally dominated the cercozoan communities, despite strong differences of community composition between Saint John's and Barrier Lake North and South. These taxa are common to be diverse and abundant in soil (Dumack et al. 2017;Fiore-Donno et al. 2018Degrune et al. 2019;Jauss et al. 2020;Ö ztoprak et al. 2020). However, the high abundance of OTUs, which were assigned to Thaumatomonas, Spongomonadidae and Woronina pythii was unexpected, since they occur only rarely or were found to be entirely absent in forests Voss et al., 2019; Fig. 1).
We found a number of indicator taxa for earthworm invasion ( Fig. 4; Supplementary Table 5). Among the 124 taxa that indicated earthworm invasion 39 taxa did not indicate any other treatment. Among those were  numerous OTUs referring to small bacterivorous flagellates of the Cercomonadida and Glissomonadida that might profit from earthworm invasion by increased dispersal or increased availability of carbon compounds. Interestingly a surprisingly high number (seven) of undetermined OTUs was found among the indicator taxa for earthworm invasion. Cercozoan OTUs that can not even be assigned to an order are common but usually lowly abundant and it might be presumed that especially these OTUs found here have not yet been characterized in the recent years of taxonomic revisions of Cercozoa as they might be a part of the earthworm microbiome that needs further investigation. However, it has to be mentioned that relationships found in an observational study are only correlative, and we cannot entirely rule out that the specific biotic and abiotic environment may have driven the distribution of earthworms. To keep such effects at a minimum, we investigated forests that have been studied for almost three decades using transect sampling methods that allow for a small-scale localization of invasion fronts as well as identification of potential dynamics of invasion fronts (Scheu 1994;Straube et al. 2009;Schlaghamerský et al. 2014;Eisenhauer et al. 2019). This is further supported by two recent meta-analyses investigating the effects of earthworm invasion on soil biotic and abiotic parameters, respectively (involving the present forests as well; Ferlian et al. 2018Ferlian et al. , 2020. Both studies analysed the influence of study type and found marginal differences between lab and (observational and experimental) field studies, indeed suggesting an effect of earthworms on soil biotic and abiotic parameters. Especially noteworthy is the high number of sequence reads assigned to W. pythii in our data set. Fiore-Donno et al. (2020) did not find this species in any of 150 German forest sites, while W. pythii was quite abundant in nearby grasslands. In the present study, we found it to occur in very low abundance in uninvaded sites, but with strongly increased relative abundance in invaded sites. W. pythii is a hyperparasite of ooymcetes, which are among the most devastating tree pathogens on the American continent (Jung et al. 2005;Garbelotto and Hayden 2012;Vargas and Negron-Ortiz 2013;Guajardo et al. 2019 Data availability Sequencing data of the protists were deposited at the European Nucleotide Archive under the accession number ERS2039495 (SAMEA104421553).
Code availability The code for the statistical analysis is available online on Github (https://github.com/Kenneth-Dumack).

Declarations
Conflict of interest Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.