DNA metabarcoding reveals the impact of Cu2+ on soil cercozoan diversity

Although


Introduction
Copper (Cu 2+ ) is an essential element in plants and animals and many microorganisms (Adriano, 2005;Stern et al., 2017); however, large amounts of the heavy metal may pose an issue for agricultural soils (Ballabio et al., 2018).Cu 2+ enters the soil through anthropogenic activities such as pig slurry (Christie and Beattie, 1989), fertilizers (Baker and Senft, 1995), pesticides (Vasseur et al., 1988), and waste from industrial and wastewater facilities (Tabatabai, 1976;Giller et al., 1998).In Europe for example, the mean Cu 2+ concentration in soil samples is approximately 16.9 mg kg − 1 , but this value is much greater in some agricultural soils, especially in vineyards of southern Europe that have a mean Cu 2+ concentration of approximately 49.3 mg kg − 1 (Ballabio et al., 2018).The total Cu 2+ concentration in soil varies depending on the soil type, mineral properties and how the soil was formed, however, the optimal Cu 2+ concentration range is possibly between 5 and 30 mg kg − 1 as lower Cu 2+ concentrations can potentially have negative impacts on crops (Adrees et al., 2015;Cornu et al., 2017;Glibota et al., 2019).Nonetheless, elevated concentrations of Cu 2+ can affect various organisms in soil food webs from microorganisms (Nunes et al., 2016) to earthworms (Van Zwieten et al., 2004).Carley et al. (2020) in a study focusing on Cu 2+ nanoparticles, which supposedly have less environmental impact than conventional pesticides, found that exposure to Cu 2+ nanoparticles affected a variety of taxonomic groups.In particular, soil microorganisms are more sensitive to heavy metal stress compared to fauna and flora that inhabit the same soil (Giller et al., 1998;Carley et al., 2020).
Protozoa encompass a taxonomically diverse group of unicellular eukaryotes that occur in large numbers in soil, ranging on average from 10 3 to 10 6 cells per gram of soil (Adl and Coleman, 2005;Ekelund and Rønn, 1994).Protozoa are also ecologically diverse and have a wide range of functional roles in the ecosystem that may contrast or complement one another (Pawlowski et al., 2012;Dumack et al., 2020).In particular, soil Protozoa are important in decomposition pathways as they consume and regulate a wide variety of microorganisms (Clarholm, 1981;Ekelund and Rønn, 1994;De Ruiter et al., 1995;Ekelund, 1998;Adl and Gupta, 2006), and stimulate microbial activity and nutrient cycling (De Ruiter et al., 1995;Bonkowski and Clarholm, 2012).Hence, they have an important role in plant primary production (Adl, 2003;Jassey et al., 2015;Geisen et al., 2018).
Research on soil Protozoa, however, is still scarce compared to other soil organisms, such as bacteria that have been investigated in greater detail; many protozoan lineages have yet to be investigated (Dassen et al., 2017;Islam et al., 2020).A major reason for the lack of research on soil Protozoa is methodological challenges that have impeded earlier research.Recent advances in molecular biology though have made it more achievable (Geisen et al., 2017), primarily through high throughput sequencing of amplified marker genes, known as metabarcoding (Pawlowski et al., 2012;Geisen et al., 2017).As opposed to traditional (non-molecular) approaches, the high resolution of molecular techniques now allows for the detection of minor changes in the diversity, which could potentially impact soil functioning.Understanding the diversity and community structure of soil Protozoa, and how abiotic factors may affect community composition, are still topics that need in-depth investigation (Geisen et al., 2017).With advancements in molecular methods such as the accessibility of new primers and improved sequencing platforms, protozoan communities can now be more readily assessed and analysed (Harder et al., 2016;Fiore-Donno et al., 2018).
We focused on a particular group of Protozoa, the Cercozoa, and on a particular heavy metal, Cu 2+ .Cercozoa, a major protozoan phylum, encompass a morphologically, functionally and ecologically diverse group adapted to both marine and terrestrial habitats (Bass et al., 2009;Ploch et al., 2016), however, the environmental diversity of Cercozoa is relatively unexplored (Fiore-Donno et al., 2018).Cercozoans are largely bacterial feeders (Oliverio et al., 2020), known to regulate the flow of resources to microorganisms at lower trophic levels (Delgado-Baquerizo et al., 2020;Jiao et al., 2021).In recent years, communities of soil Cercozoa have been investigated (Fiore-Donno et al., 2019); however, it is not known exactly how a pollutant can impact cercozoan populations.
We examined how Cu 2+ affects the diversity and community composition of soil Cercozoa.Using molecular techniques, this study aimed to determine how an array of Cu 2+ concentrations affect soil Cercozoa over a 43-day period.We detected the changes in the cercozoan community and taxonomic composition using deoxyribonucleic acid (DNA) metabarcoding and compared this to protozoan abundance measured by a conventional microbiological technique.We addressed the following two hypotheses: (1) Cercozoan alpha diversity correlates negatively with Cu 2+ addition and length of exposure.(2) Some cercozoan groups are more Cu 2+ sensitive and are most likely to decrease in abundance or disappear with increased Cu 2+ concentrations and length of exposure.

Respiration and protozoan measurements
Both Cu 2+ concentrations and incubation time significantly affected CO 2 production in the microcosms (Fig. 1).Respiration was largely unaffected at the lowest Cu 2+ concentrations, while respiration of the higher Cu 2+ concentrations differed considerably from the lower.
The total number of Protozoa significantly decreased with time and with increasing Cu 2+ concentrations (Fig. 2).

Cercozoan sequence data and composition
In total, we recorded 9,899,296 reads distributed among 5,345 ASVs.After taxonomic assignment, we identified 8,701,072 reads and 4,419 ASVs, which was further filtered and resampled leaving 5,201,853 (52.5 %) reads and 4,402 (82.4 %) sequences.Altogether, we identified 34 cercozoan families and 14 orders comprising 76.5 % and 84.4 % of the total cercozoan sequences respectively.Supplementary data (Tables 1-3) illustrates the total number of reads for each cercozoan taxonomic group.

Cercozoan orders
Cercozoan orders with the highest relative read abundance were Cercomonadida, Crycomonadida, Euglyphida, Glissomonadida, Limnofilida, Plasmodiophorida and Spongomonadida (Fig. 3).Glissomonadida and Cercomonadida had the highest relative read abundances across all Cu 2+ concentrations on both sampling days.Some orders such as Spongomonadida and Limnofilida occurred scarcely at the higher Cu 2+ concentrations, albeit with some recovery at day 43.Supplementary Table 1 presents the total abundance of all cercozoan orders.

Cercozoan families
The most abundant families identified included Allapsidae, Cercomonadidae, Limnofilidae, Paracercomonadidae, Plasmordiophoridae, Rhogostomidae, Sandonidae and Spongomonadidae (Fig. 4).Sandonidae was the most common family compromising a majority of the relative read abundance in each sample.Some families occurred in all samples (Supplementary Table 2), but in some samples they were not abundant enough to be represented in Fig. 4. Of these, notably Limnofilidae, Spongomonadidae, Paracercomonadidae and Rhogostomidae that occurred in low numbers at the highest concentrations, albeit with some recovery on day 43.In particular, Plasmordiophoridae was scarce at the highest concentrations (625 to 10,000 mg Cu 2+ kg − 1 dw soil).Overall, the total abundance of all families decreased substantially after the 10,000 mg Cu 2+ kg − 1 dw soil treatment and had not recovered by day 21, although both the abundance and number of families had started to recover by day 43 (Supplementary Table 2).

Cercozoan genera
For a large proportion of the ASVs, we could not assign them to a named genus, and many of the genera that we recorded had abundances so low that the genus distribution showed few distinct patterns Fig. 1.Soil respiration (accumulated CO 2 ) in microcosms amended with different Cu 2+ concentrations and measured on day 21 (red dots) and day (blue dots) after experimental start.The solid lines show fits of the data to a simple 3-parameter sigmoid equation, y = a/(1 + exp(-(x-x 0 )/b)), often used to model eco-toxicological data (Johansen et al., 2018).A two-way ANOVA showed that Cu 2+ (log-transformed) affected soil respiration significantly negatively (p < 0.0001) and that there was significantly more accumulated CO on day 43 than on day 21 (p = 0.0005).There was no significant interaction between Cu 2+ and day (p = 0.32).We assessed homoscedasticity and normality of the model by residual and quantile plots (Supplementary Fig. 2A).The x-axis shows how much Cu 2+ we added to the soil, the amount already present in the soil was 2.2 mg Kg − 1 , and is negligible compared to the addition.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) L. Lourenço et al. (Supplementary Table 3, Supplementary Fig. 1.).The most abundant genera were Eocercomonas, Limnofila, Paracercomonas, Polymyxa, Sandona, and Spongomonas.Both Sandona and Paracercomonas displayed a bimodal distribution with high abundances at high and low Cu 2+ concentrations (Supplementary Table 3).Some of the less abundant genera such as Arachnula, Orciraptor and Sainouron were only found in one single sample.

Richness patterns and community composition
After resampling and cleaning the data, one outlier was excluded at 10,000 mg Cu 2+ kg − 1 dw soil as there was one microcosm that upon inspection of our laboratory protocol received less Cu 2+ than intended.The calculated Chao1 values were almost identical to the observed ASV richness showing that sequencing depth was adequate to capture total richness.However, we chose to use the Chao1 values for further analysis as these are better at taking into account the low abundance ASVs as compared to observed ASV richness (Chao and Chiu, 2016).Results showed that both Cu 2+ concentration and time affected cercozoan richness significantly (p < 0.005, two-way analysis of variance (ANOVA)).In addition, the interaction between the factors had a significant effect on the sampling date and treatment (Fig. 5a).The Shannon diversity index displayed another pattern than the Chao1 (Fig. 5b).Still, a two-way ANOVA showed that Cu 2+ concentration and sampling day affected Shannon diversity (p < 0.05), whereas the interaction between the factors had no significant effect.

NMDS ordination
To analyse the community composition more thoroughly, we constructed an NMDS ordination plot to illustrate the shift in cercozoan diversity with increasing Cu 2+ concentrations (Fig. 6).Only little change took place between 0.0, 0.04 and 39.1 mg Cu 2+ kg − 1 dw soil, a slight shift occurred at 156.3 and 625 mg Cu 2+ kg − 1 dw soil, whereas on both sampling days, at the highest treatments (2,500 and 10,000 mg Cu 2+ kg − 1 dw soil), the shift was most pronounced.We noticed that with a higher Cu 2+ concentration, the development in individual microcosms becomes more stochastic, as the triplicates of higher Cu 2+ concentrations showed more dispersion compared to lower concentrations.Interestingly, at 2,500 mg Cu 2+ kg − 1 dw soil the community on day had a greater change in community diversity compared to the community found on day 21.While at 10,000 mg Cu 2+ kg − 1 dw soil, day showed a greater shift in community structure compared to the other concentrations.Pairwise comparisons verify the beta diversity shifts displayed in Fig. 6, as most concentrations were found to be significantly , often used to model eco-toxicological data (solid lines).A two-way ANOVA showed that Cu 2+ affected protozoan abundance significantly negatively (p < 0.0001) and that there were significantly fewer Protozoa on day 43 than on day 21 (p = 0.011).We found no significant change in Cu 2+ effect with time (p = 0.26).We assessed homoscedasticity and normality of the model by residual and quantile plots (Supplementary Fig. 2B).Please notice that due to the shortcoming of the MPN method to produce exact results for small protozoan densities, as an artefact; all numbers for the two highest concentrations were similar.However, for reasons of clarity, we placed results for day 21 slightly above those for day 43.The x-axis shows how much Cu 2+ we added to the soil, the amount already present in the soil was 2.2 mg Kg − 1 , and is negligible compared to the addition.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 3. Relative read abundance (%) of the top cercozoan orders identified across all Cu 2+ concentration treatments after 21 and 43 days in soil microcosms.Orders that comprised of less than 5 % of the total relative abundance are denoted by "Other (<5%)".

Microbial activity and MPN-based protozoan counts
Both microbial activity, as indicated by CO 2 production, and the MPN-based protozoan counts indicated a depression of microbial populations and activity with increasing Cu 2+ concentration, in line with other microcosm studies, such as Ekelund et al. (2003) in another Cu 2+ study and Holtze et al. (2003) in a mercury study.This contrasts with the situation in natural soils, where the Cu 2+ concentration seems to have only a slight effect on bacteria and Protozoa.Boon et al. (1998) attribute this to the presence of vegetation, as plant growth promotes increased biological activity.In particular, tolerant plants can provide a food base and reverse the Cu 2+ effects.Most microbial activity (ie.respiration) is caused by bacteria (Ekelund and Rønn, 1994).Therefore, the decrease of microbial activity with increasing Cu 2+ concentrations is most likely the result of i) a direct toxic effect on the bacteria combined with, ii) the   3).We used Chao1 and Shannon diversity to visualize the community richness and diversity (Hughes et al., 2001;Haegeman et al., 2012).
indirect effect from the reduction of protozoan stimulation of the bacterial activity.
Most studies on Protozoa have been conducted in aqueous solutions (Johansen et al., 2018).Often, toxicological effects on the microorganisms are stronger in liquids than in soil as heavy metals are absorbed into the soil matrix.Moreover, while microbial numbers were unaffected by a wide range of Cu 2+ concentrations, microbial composition was changed.In particular, there was a shift from gram-negative to grampositive bacteria (Ekelund et al., 2003;Wang et al., 2023).Conversely, Nunes et al. (2016) found that the prevalence of bacterial groups such as Nitrospira and Acidobacteria increased with Cu 2+ , while Proteobacteria and Actinobacteria decreased.On the species level, in particular Arthrobacter spp.are highly sensitive to Cu 2+ , while Bacillus spp.are more tolerant (Asatiani et al., 2021).These variations in results could be due the type of soil (being soil texture and organic content).For example, Actinobacteria was positively correlated with Cu 2+ in clayloam soils, but negatively correlated in silt-loam soils (Liang et al., 2017).
Cu 2+ may affect soil Protozoa in different ways.It may be a driving force in the soil system and affect the abundance and diversity of bacterial communities (Nunes et al., 2016).This will then impact the protozoan communities that rely on these bacteria.Protozoa may also be affected directly by oxidative stress.Heavy metals, such as Cu 2+ , may induce oxidative stress through the increased synthesis of reactive oxygen species.In particular, Cu 2+ generates peroxides at the nuclear and cytoplasmic levels, which can contribute to cell death.Different Protozoa have different resistance levels (Rico et al., 2009).Luu et al. (2022) examined ciliated Protozoa in loam-sandy soil and found that Cu 2+ treatment increased the number of tolerant species, while sensitive species were absent at high treatment levels.Furthermore it is possible that soil ciliates may be better able to exploit Cu 2+ levels compared to other Protozoa (Ekelund et al., 2003;Luu et al., 2022).

The most common Cercozoa in our study
The most abundant orders we found (Fig. 3) correspond with the identified sequences from grassland samples reported by Fiore-Donno et al. (2019).The majority of the Cercozoa we identified were surfaceassociated gliding or creeping bacterivorous flagellates, Glissomonadida (Howe et al., 2011), Cercomonadida (Bass et al., 2009), and Cryomonadida (Cavalier-Smith and Chao, 2003), all well adapted to, and commonly occurring in soil.It may appear more surprising that representatives from Spongomonadida occurred commonly in the samples.The known members of this group spend a majority of their life in colonies (Dumack et al., 2020), which is a lifestyle that does not agree with our microcosm systems.Still, other authors have reported Spongomonadida-sequences from soil samples, for example Fiore-Donno et al. ( 2019), in addition to yet undescribed free-living species, such as the discovery of non-sedentary bicosoecids in soil samples (Harder et al., 2014).
Another seemingly surprising finding was that members of Plasmodiophorida occurred commonly, although only among the most common in two samples (Supplementary Table 1, Fig. 3).They displayed an overall decreasing trend to increasing Cu 2+ concentration.Some Plasmodiophorida are economically important plant pathogens, which cause root diseases in potatoes and watercress, for example, and have the ability to transmit a range of plant viruses (Down et al., 2002;Kanyuka et al., 2003;Falloon et al., 2016).The occurrence of Plasmodiophorida in our soil may be due to their ability to form resting spores.However, some may parasitize algae and oomycetes (Braselton, 1995;Adl et al., 2012;Neuhauser et al., 2014).Plasmodiophorida ASVs in our study could be assigned to Polymyxa graminis, Spongospora nasturtii, Spongospora subterranean and Woronina pythii, where the latter is known to infect oomycetes (Dylewski and Miller, 1983).We notice that in previous studies of terrestrial Cercozoa (Harder et al., 2016;Fiore-Donno et al., 2019), the testate amoeba family Rhogostomidae (Thecofilosea), dominated.Rhogostomidae was also a major family found in our study, but decreased with increasing Cu 2+ and with time, albeit with a noticeable spike at 2,500 mg Cu 2+ kg − 1 dw soil (Supplementary Table 2).
The majority of the genera in the samples (Supplementary Table 3) were bacterivorous.However, besides the already mentioned members of the family Plasmophorididae, we also noticed the presence of Orciraptor, Viridiraptor and Vampyrella, which prey on (mainly green) algae.None of these genera showed any particular response to the Cu 2+ concentrations.Conversely, Proleptomonas, a monotypic genus encompassing only P. faecicola, which occurred relatively frequently in all samples, showed a lower frequency at the two highest Cu 2+ concentrations.As a non-phagotrophic coprophilic flagellate, Proleptomonas is unique among the cercozoans, and was originally described as a free-fiving Fig. 6.An NMDS ordination of the cercozoan community comparing both day 21 and day 43 in all Cu 2+ concentration treatments.In addition to NMDS ordination in three dimensions using metaMDS (based on Bray-Curtis dissimilarities calculated on a Hellinger transformed ASV table) and both a pairwise factorfit and a permanova was conducted.

Cercozoan alpha diversity
Cu 2+ negatively affected Chao1 taxa richness, though only for the highest concentrations and in particular on day 43 (Fig. 5a).This Cu 2+ effect agrees with previous studies of mixed protozoan communities exposed to heavy metals in soil (Foissner, 1991;Ekelund and Rønn, 1994;Domonell et al., 2013) and the effect on testate amoebae in moss (Nguyen-Viet et al., 2006).In line with our findings, Ekelund et al. (2003) found a more pronounced effect of Cu 2+ on alpha-diversity than on abundance, probably because some groups of Protozoa are more susceptible than others to toxic compounds (Ekelund et al., 2000;Ekelund, 2003).
Shannon diversity responded in a more irregular way compared to Chao1.After a steady decrease, Shannon diversity suddenly increased at 2,500 mg Cu 2+ kg − 1 dw soil (Fig. 5b), without an accompanying large quantitative shift in richness, even though some compositional changes in cercozoan groups did occur.Carley et al. (2020) have suggested that turnover of taxa is a driving factor for increasing diversity.Hence, we suggest that the shift in the community, reflected in the increase in Shannon diversity (Fig. 5b), was caused by a turnover of taxa, potentially in the direction of taxa tolerant to extreme conditions.This is also reflected in the very noticeable shift at 2,500 mg Cu 2+ kg − 1 dw soil on the ordination plot (Fig. 6), where the points move away and become more scattered.
Soil Protozoa are generally negatively affected when exposed to heavy metals and other xenobiotics, but some opportunistic genera such as Colpoda, Spumella, Bodo and Ceremonas may increase (Domonell et al., 2013;Ekelund and Rønn, 1994;Foissner, 1991;Mondragón-Carmarillo et al., 2020).These genera are some of the most abundant and common microbial predators in soils and are important in the transformation of nutrients from bacterial biomass through the trophic network (Domonell et al., 2013;Ekelund and Rønn, 1994;Foissner, 1991).Hence the overall Cu 2+ effect on soil function may be less pronounced.Mondragón-Carmarillo et al. (2020) suggested that generalist Protozoa might have a more effective survival strategy compared to specialists in systems exposed to xenobiotics, as the depletion of specialist genera after the exposure allows generalists to increase due to the reduction of competition.Although lower levels of Cu 2+ have been found to inhibit microbial growth impacting the protozoan population (Ekelund et al., 2003), the higher Cu 2+ concentrations had greater impacts on the protozoan community structure, thus allowing tolerant microorganisms to increase.As such, more generalist genera would be able to survive and increase in abundance compared to those with specialist strategies.Mondragón-Carmarillo et al. (2020) also suggested that the decrease in abundance and species richness of Protozoa could have been related to the death of sensitive individuals and the reduction of food sources (see also Rykiel, 1985), thereby simplifying the protozoan community.Therefore, depending on the responses of a cercozoan group to Cu 2+ , the life strategies and the overall impact on the soil food web, different cercozoan groups would be more favoured compared to others.In our study, we found that bacterivorous families such as Allapsidae, Sandonidae and Cercomonadidae increased in abundance compared to non-bacterivore families such as Rhogostomidae and Plasmophorididae, (Supplementary Tables).As Cercozoa are considered to be mainly bacterial feeders (Oliverio et al., 2020), in this study, we could consider bacterivore cercozoan groups as the generalists and nonbacterial feeders as the specialists.

Conclusion
It has long been accepted that trace metals may affect microbial populations, but with the implementation of molecular technologies, we can study this impact in much greater detail.Our study illustrates that an increase in heavy metals, such as Cu 2+ , has an overall negative impact on cercozoan richness and diversity.However, the changes in diversity, and shift in the cercozoan community, that we saw, suggest that there was a taxa turnover at higher Cu 2+ concentrations towards taxa that were more tolerant.To understand the impact of Cu 2+ and other heavy metals on soil Cercozoa in more detail, identified strains or cultures from specific families such as Allapsidae, Cercomonadidae, Rhogostomidae and Spongomonadidae (Bass et al., 2009;Howe et al., 2009;Cavalier-Smith and Chao, 2003;Dumack et al., 2017) could be used in further toxicological studies.To our knowledge, this is the first in-depth molecular-based study that investigates the impact of heavy metals on a cercozoan community, with the potential use of aiding in further explorations of how such changes in the community could impact agricultural soils.

Soil microcosms
We used coarse sandy agricultural soil from Jyndevad Field Station (Aarhus University), with an organic content of 4.2 %, pH of 5.9, and a Cu 2+ content of 2.2 mg kg − 1 (Henriksen, 1959).Samples from the top layer were air-dried, homogenised and sieved (<2mm; Ekelund et al., 2003).We mixed 1 g of dry soil with 1.125 mg of finely ground barley straw.We prepared microcosms in 116 ml serum bottles, each with 8.0 g of the soil/straw mixture.We then made a series of aqueous solutions of copper(II) chloride (CuCl 2 ) starting with 14.909 mg of CuCl 2 in 40 ml of water to make the highest concentration, followed by a fourfold dilution, except for the three lowest concentrations, which were prepared separately.When added to the microcosms, we obtained final concentrations of 0, 0.04, 0.2, 0.6, 2.4, 9.8, 39.1, 156.3, 625, 2500 and 10 000 mg kg − 1 dry weight (dw) soil of added Cu 2+ to our microcosms.We referred to other studies for this logarithmic scale, which is commonly done in ecotoxicology (Ekelund et al., 2003;Holtze et al., 2003;Johansen et al., 2018).CuCl 2 is very soluble and spreads uniformly; as such the Cu 2+ solutions were evenly spread across the soil surface using a pipette, and each dilution had six replicates.The microcosms were then sealed and incubated in the dark at 15 • C.
We measured soil respiration, as carbon dioxide (CO 2 ) content, in the microcosms 21 and 43 days after incubation start.Each sampling consisted of triplicate microcosms from each of the eleven Cu 2+ treatments.CO 2 was measured by injecting 5 ml of gas from the microcosms' headspace into a gas chromatograph equipped with a thermal conductivity detector (Mikrolab Aarhus, Denmark).

Protozoan measurements
After CO 2 measurement, we opened the microcosms and promptly transferred 1.0 g of soil taken from different parts of the microcosm to an Eppendorf tube and placed it at − 20 • C until further processing.For each sample, we conducted a threefold dilution series culturing technique using 96-well microtiter plates.The microtiter plates were incubated in the dark at 15 • C and inspected for the presence or absence of Protozoa using an inverted microscope after 1 and 3 weeks of incubation (Darbyshire et al., 1974), as modified by Rønn et al. (1995).This allowed us to estimate the concentration of viable Protozoa by the 'most probable number' (MPN) method a using a computer program (Ekelund unpublished, based on Fisher, 1922) to convert the microtiter plate patterns to most probable protozoan numbers.

Soil DNA extraction and sequencing
Based on distinct differences in accumulated CO 2 , we selected the microcosms with the Cu 2+ concentrations 0.0, 0.04, 39.1, 156.3, 625, 2,400 and 10,000 mg Cu 2+ kg − 1 dw soil for molecular analysis.We used DNeasy PowerSoil Pro (Qiagen, Germany) to extract DNA according to the manufacturer's instructions and amplified the V4 region of the 18S gene using the previously designed specific cercozoan primers, Cerc479F and Cerc750R (Harder et al., 2016).Polymerase chain reaction (PCR) was conducted in triplicate for each DNA extract.The PCR setup and thermocycling followed Harder et al. (2016), and the PCR products were visualised on a 2 % agarose gel and stored at − 20 • C until further processing.Both forward and reverse primers were constructed with 96 unique tags (MID/barcodes) of 6 bp at the 5′ end using a restrictive dual indexing approach, where no primer tag (forward or reverse) was used more than once in any sequencing library, and no specific combinations of forward and reverse tag were reused.We scored the PCR bands as weak, medium, or strong, and used 7.5 μl, 5 μl or 2.5 μl respectively, in the pooling of PCR products before library building.We built sequencing libraries using the Illumina TruSeq PCR-Free kit (Illumina, Denmark) according to the manufacturer's instructions.We performed sequencing at the Danish National High Throughput DNA Sequencing Centre, on the Illumina Miseq v.3 platform (Illumina) with samples divided into 5 libraries, 300 bp paired-end runs.

Sequence processing
Bioinformatic steps followed the general procedures of Frøslev et al., (2017Frøslev et al., ( , 2019)), with minor modifications.Demultiplexing of samples was done with a custom script that keeps R1 and R2 separate for divisive amplicon denoising algorithm 2 (DADA2) processing and is based on Cutadapt (Martin, 2011) searching for a sequence pattern matching the full-length combined tag and primer allowing for no errors, and removing possible remnants of the other primer at the 3′ end.We used DADA2 (v 1.8; Callahan et al., 2016) to identify amplicon sequence variants (ASVs) and remove chimeras (bismeras).We used ASVs as the molecular taxonomic unit to allow for better reproducibility of downward analyses, and not to impose any arbitrary ideas about molecular species limits (Callahan et al., 2017).Sequences were filtered and matched between R1 and R2 reads with DADA2.Taxonomic annotation of the ASVs was done with the tool implemented in DADA2 using the protist ribosomal 2 (PR2) reference database (Guillou et al., 2013).Sequences were then resampled down to the second quartile, triplicates per sample were combined, and an outlier was removed due to too few ASVs.

Statistical analysis
We conducted statistical and taxonomic analysis in R (v 3.6; R Development Core Team, 2014), excluding non-cercozoan ASVs prior to completing the richness and diversity plots and analyses.Plots were done with the vegan package v2.5-7 (Oksanen et al., 2020) or in Sig-maStat® (Systat Software Inc.).Taxonomy abundance plots were created using the phyloseq package v 1. 3.2 (McMurdie and Holmes, 2013).Furthermore, a nonmetric multidimensional scaling (NMDS) ordination plot was created using metaMDS on a Hellinger transformed ASV table using Bray-Curtis dissimilarities.Ordination was done in three dimensions but plotted on the first two.One 10 000 mg Cu 2+ kg − 1 dw soil sample was lost and therefore could not be added to the analysis.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. Protozoan abundance in microcosms amended with six different Cu 2+ concentrations and monitored on day 21 (red dots) and day 43 (blue dots) after experimental start.We fitted the data to a simple 3-parameter sigmoid equation, y = a/(1 + exp(-(x-x0)/b)), often used to model eco-toxicological data (solid lines).A two-way ANOVA showed that Cu 2+ affected protozoan abundance significantly negatively (p < 0.0001) and that there were significantly fewer Protozoa on day 43 than on day 21 (p = 0.011).We found no significant change in Cu 2+ effect with time (p = 0.26).We assessed homoscedasticity and normality of the model by residual and quantile plots (Supplementary Fig.2B).Please notice that due to the shortcoming of the MPN method to produce exact results for small protozoan densities, as an artefact; all numbers for the two highest concentrations were similar.However, for reasons of clarity, we placed results for day 21 slightly above those for day 43.The x-axis shows how much Cu 2+ we added to the soil, the amount already present in the soil was 2.2 mg Kg − 1 , and is negligible compared to the addition.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. Relative read abundance (%) of the most abundant cercozoan families across all Cu 2+ concentrations after 21 and 43 days in soil microcosms.Families that comprised of less than 10 % of the total relative abundance are denoted by "Other (<10 %)".

Fig. 5 .
Fig. 5. Chao1 richness and Shannon diversity comparing all Cu 2+concentration treatments between both day 21 and 43 of incubation.A two-way ANOVA illustrated that the diversity of both Chao1 and Shannon was impacted by concentration and day of sampling.We assessed homoscedasticity and normality of the model by residual and quantile plots (Supplementary Fig.3).We used Chao1 and Shannon diversity to visualize the community richness and diversity(Hughes et al., 2001;Haegeman et al., 2012).