Molecular Ecology (2014) 23, 3341–3355

doi: 10.1111/mec.12819

Effects of long-term differential fertilization on eukaryotic microbial communities in an arable soil: a multiple barcoding approach GUILLAUME LENTENDU,*† ‡ TESFAYE WUBET,‡ § ANTONIS CHATZINOTAS,¶ CHRISTIAN WILHELM,*§ FRANC ß O I S B U S C O T ‡ § and M A R T I N S C H L E G E L † § *Plant Physiology, Institute of Biology, University of Leipzig, Johannisallee 21-23, Leipzig 04103, Germany, †Molecular Evolution and Animal Systematics, Institute of Biology, University of Leipzig, Talstraße 33, Leipzig 04103, Germany, ‡Department of Soil Ecology, UFZ – Helmholtz Centre for Environmental Research, Theodor-Lieser-Str. 4, Halle/Saale 06120, Germany, §German Centre for Integrative Biodiversity Research (iDiv), Halle-Jena-Leipzig, Deutscher Platz 5e, Leipzig 04103, Germany, ¶Department of Environmental Microbiology, UFZ – Helmholtz Centre for Environmental Research, Permoserstraße 15, Leipzig 04318, Germany

Abstract To understand the fine-scale effects of changes in nutrient availability on eukaryotic soil microorganisms communities, a multiple barcoding approach was used to analyse soil samples from four different treatments in a long-term fertilization experiment. We performed PCR amplification on soil DNA with primer pairs specifically targeting the 18S rRNA genes of all eukaryotes and three protist groups (Cercozoa, ChrysophyceaeSynurophyceae and Kinetoplastida) as well as the ITS gene of fungi and the 23S plastid rRNA gene of photoautotrophic microorganisms. Amplicons were pyrosequenced, and a total of 88 706 quality filtered reads were clustered into 1232 operational taxonomic units (OTU) across the six data sets. Comparisons of the taxonomic coverage achieved based on overlapping assignment of OTUs revealed that half of the eukaryotic taxa identified were missed by the universal eukaryotic barcoding marker. There were only little differences in OTU richness observed between organic- (farmyard manure), mineral- and nonfertilized soils. However, the community compositions appeared to be strongly structured by organic fertilization in all data sets other than that generated using the universal eukaryotic 18S rRNA gene primers, whereas mineral fertilization had only a minor effect. In addition, a co-occurrence based network analysis revealed complex potential interaction patterns between OTUs from different trophic levels, for example between fungivorous flagellates and fungi. Our results demonstrate that changes in pH, moisture and organic nutrients availability caused shifts in the composition of eukaryotic microbial communities at multiple trophic levels. Keywords: community ecology, eukaryotic microorganisms, fertilization, metabarcoding, protists, soil biodiversity Received 28 February 2014; revision received 27 May 2014; accepted 28 May 2014

Introduction Eukaryotic microorganisms represent an important part of the soil microbial community (Adl & Gupta 2006; Coleman & Whitman, 2005). Traditionally, soil eukaryCorrespondence: Guillaume Lentendu, Fax: 49 (0) 345 558 5449; E-mail: [email protected] © 2014 John Wiley & Sons Ltd

otic microbial communities have been investigated using cultivation- or microscopy-based approaches (Esteban et al. 2006; Zancan et al. 2006; Yurkov et al. 2012). Moreover, studies on eukaryotic microorganism diversity and ecology have focused on a limited number of functional groups (Ekelund et al. 2001; Domonell et al. 2013). Recent investigations, however, indicate that eukaryotic microorganisms exhibit high levels of

3342 G . L E N T E N D U E T A L . molecular diversity, functional diversification and habitat specificity (Boenigk et al. 2005; Von der Heyden & Cavalier-Smith 2005). In contrast to bacterial communities, cultivation-independent rRNA gene–based techniques have not yet been widely used to analyse eukaryotic microbial communities in soils (Lawley et al. 2004; Lara et al. 2007b; Glaser et al. 2014). Despite the high taxonomic resolution of isolation-based approaches and the wealth of new sequences revealed by comparative analysis of 18S rRNA gene libraries, none of these approaches can cover a sufficiently high number of individuals or sequences to provide a comprehensive overview of eukaryotic microbial communities. Although the usefulness of eukaryotic microorganisms as bioindicators, particularly in agro-ecosystems, has been acknowledged (Foissner 1999; B
erard et al. 2005), a fast and reliable analytical method for their analysis within hundreds of samples remains to be developed. More comprehensive approaches based on analysing the 18S rRNA gene barcoding markers of whole eukaryotic communities using next-generation sequencing (NGS) and bioinformatics (Stoeck et al. 2009; Cheung et al. 2010; Nolte et al. 2010; Sugiyama et al. 2010) enable the analysis of large sequence data from a wide spectrum of eukaryotic microbial groups to uncover the richness and diversity of operational taxonomic units (OTUs) in hundreds of samples (i.e. alpha-diversity) and relate patterns of community composition to environmental gradients and changes (i.e. beta-diversity; Bik et al. 2012). Recently, the introduction of network and co-occurrence analyses in environmental microbial ecology shed light on potential interactions in complex microbial assemblages (Horner-Devine et al. 2007; Fuhrman 2009; Steele et al. 2011). However, the high proportion of plant, animal and fungal sequences within the soil 18S rRNA gene pool (62.9% in Urich et al. 2008; 55;.3% in Baldwin et al. 2013) presents challenges in the study of soil eukaryotic microorganisms using universal eukaryotic barcoding markers that are not encountered when dealing with aquatic systems (Stoeck et al. 2009) or soil crusts (Meadow & Zabinski 2012). A different approach is clearly required to cover the full taxonomic diversity of eukaryotic microorganisms in samples from more complex environments such as arable soils. Soils have been treated with organic and inorganic fertilizers for over a century to increase agricultural yields (Edmeades 2003; Merbach & Schulz 2012). Although the positive effects of fertilization on nutrient availability and plant growth are obvious, little is known about the impact of long-term fertilization on belowground microorganisms. Previous studies have identified changes in microbial biomass and community structure related to both organic and inorganic fertilization (Marschner 2003; B€ ohme et al. 2005; Lejon et al.

2007), but these studies suffer from lack of resolution or specificity for eukaryotic microorganisms. Eukaryotic microorganisms fulfil diverse functions in soil food webs and nutrient cycles, including primary production (e.g. microalgae), predation (e.g. protists) and decomposition (fungi; Bonkowski 2004; Botha 2011). Because fertilization induces changes in the relative abundance and availability of nutrients as well as the soil microbial biomass (Houot & Chaussod 1995; B€ ohme et al. 2005), eukaryotic microbial communities and their interactions in intensely fertilized soils should be expected to differ significantly from those in unfertilized soils. Previous experiments have demonstrated that soil microorganisms can take a long time to respond to changes in the soil’s management regime, indicating that long-term experiments are required to properly describe the influence of management on soil communities (Hedlund et al. 2003; Eisenhauer et al. 2010). In this work, we describe the use of a multiple barcoding approach to characterize the diversity and community composition of eukaryotic microorganisms in a way that avoids masking effects of multicellular organisms. We compared universal and group-specific primers with respect to their reliability in characterizing eukaryotic soil microorganisms. We then evaluated this approach to unravel the effects of 112 years of differential fertilization on community structure using soil samples from the long-term static fertilization experiment at Bad Lauchst€ adt (Germany). Three hypotheses were tested: (i) group-specific barcoding markers enable a more precise detection of eukaryotic microorganisms compared with a universal barcoding marker; (ii) as carbon is the most important resource for most microorganisms, organic fertilization causes more pronounced shifts in community composition than mineral fertilization; and (iii) given the variety of functions fulfilled by eukaryotic microorganisms and their interdependency in the soil food web, the interactions between soil microorganisms as revealed by co-occurrence or exclusion patterns differ among and between different trophic levels.

Materials and methods Fertilization experiment and soil sampling Soil was sampled in August 2011 from the long-term static fertilization experiment at Bad Lauchst€ adt established in 1902 on a Haplic Chernozem soil (FAO) [Molisol (USDA)] (51° 240 N, 11° 530 E, http://www.ufz.de/ index.php?de=15278; Blair et al. 2006). The experiment is based on an 8-year crop rotation in which a gradient of organic (farmyard manure) and mineral (NPK) fertilization is applied over 18 treatments. The amount © 2014 John Wiley & Sons Ltd

S O I L E U K A R Y O T I C M I C R O O R G A N I S M S M E T A B A R C O D I N G 3343 of fertilizer application in relation to crop type and rotation treatment is specified in Table S1 (Supporting information); additional details are given in Blair et al. (2006). Samples were taken at the end of an 8-year crop rotation during the second consecutive year of Alfalfa (Medicago sativa) cultivation. The last fertilization prior to sampling was applied before the Alfalfa was seeded. Four different treatments were sampled: no fertilization, mineral fertilization (NPK), organic fertilization (farmyard manure) and both mineral and organic fertilization. These treatments are henceforth referred to as ‘N’, ‘M’, ‘O’ and ‘OM’, respectively. Unfortunately, no independent replicates were realized while establishing the long-term fertilization experiment 112 years ago. Pseudoreplicate samples were therefore collected by dividing each treatment plot into four sections or subplots, as discussed previously (B€ ohme et al. 2005). Three of these subplots were then selected at random for sampling. A total of 10 soil cores from the topmost 10 cm soil layer of each subplot were collected, combined, sieved through a 2-mm mesh and homogenized directly at the sampling site. The 12 resulting composite samples were transported in a cool box to the laboratory and stored at 20 °C.

Microbial metagenomic DNA extraction and amplicon library preparation Total DNA was extracted from 0.5 g soil samples using the PowerSoilâ DNA isolation kit (Mo Bio Laboratories, Inc., Carlsbad, CA) following the manufacturer’s instructions. DNA concentrations were quantified using a Nanodrop spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). Ribosomal, plastid and intergenic transcribed spacer DNA was then amplified using six different primer pairs (Table S2, Supporting Information). These primer pairs, which are specific for Cercozoa (newly designed, see Appendix S1, Table S3, Supporting Information), Chrysophyceae/Synurophyceae, Eukaryota, Fungi, Kinetoplastida and both prokaryotic (i.e. Cyanobacteria) and eukaryotic photoautotrophs, will be referred to as ‘cer’, ‘chry’, ‘euk’, ‘fun’, ‘kin’ and ‘upa’, respectively. PCR amplifications were conducted in triplicate in a total reaction volume of 25 lL containing 12,5 lL of GoTaqâ green Master Mix (Promega GMBH, Mannheim, Germany), 0.5 lL of each primer (25 lM) and 1 lL of soil DNA as template (~10 ng/lL). A seminested PCR amplification strategy was adopted for both the chry- and kin-primer pairs, as described previously for the kin-primer pair (Glaser et al. 2014). PCR conditions for all six primer pairs are provided in the supplementary information (Appendix S2, Supporting Information). Each of the specific primers contained at their 50 end the Lib-L A or B Roche adapters for forward and reverse © 2014 John Wiley & Sons Ltd

sequencing, respectively, the forward adapter being followed by a four nucleotide key and a 10 nucleotide barcode sequence. Amplicons from the triplicate PCRs were pooled and purified using the Qiaquick Gel Extraction kit (Qiagen, Hilden, Germany). The DNA concentration of the purified samples was measured using the PicoGreen dsDNA assay (Invitrogen, Carlsbad, CA, USA). Samples were pooled at equimolar concentration and sequenced unidirectionally on a Roche GS FLX+ 454 automated pyrosequencer using Titanium chemistry at the department of Soil Ecology, UFZ.

Sequence processing, taxonomic assignment and phylogenetic analyses The recovered reads were quality trimmed using MOTHUR (Schloss et al. 2009). Sequences were only retained if they had at most one MID mismatch, at most four primer mismatches, no undetermined nucleotides, an average quality value of 25 or more, a maximum homopolymer length of eight nucleotides, and a read length within 2 SD of the average value for the relevant primer pair (the minimum read lengths were 300, 265, 300, 300, 290 and 200 nucleotides in the respective data sets). In addition, only reads containing both the forward and reverse primers were kept for the upa-data set. Chimeras were removed from each data set using UCHIME (Edgar et al. 2011) as implemented in MOTHUR. To avoid sampling size effects, the number of reads per sample was normalized for each data set by randomly subsampling to the lower number of reads among samples. Reads were then clustered in operational taxonomic units (OTU) for each data set using the MCL program (Van Dongen 2000), based on pairwise global alignment distances that were calculated with the Needleman–Wunsch algorithm as implemented in MOTHUR. Threshold distances for OTU clustering corresponding to the genetic distances that maximized the separation between intraspecific and interspecific intragenus distances were determined by in silico PCR (Appendix S3, Supporting Information). This approach has previously been recommended for the analysis of prokaryotic 16S rRNA gene sequences (Kim et al. 2011). The clustering thresholds approximately representing species level were thus fixed at genetic distances of 0.04, 0.02, 0.02, 0.03, 0.05 and 0.03, respectively (Fig. S1, Supporting Information). Only OTUs with more than two reads were retained for further analyses because singletons and doubletons are likely to be artefacts in pyrosequencing data sets (Reeder & Knight 2009; Huse et al. 2010; Behnke et al. 2011). Each dereplicated sequence was taxonomically assigned using the Wang algorithm as implemented in MOTHUR by keeping at each level the taxonomy shared by more than 60% of

3344 G . L E N T E N D U E T A L . the best hits after a bootstrap assignment of 100 iterations (i.e. 60% consensus threshold). Reference sequences originated either from the PR2 (cer-, chry-, euk- and kin-data sets; Guillou et al. 2012) or UNITE (January 2013 release, its-data set; Abarenkov et al. 2010) databases, or from a self-constructed ENA sequences database (release 116, Leinonen et al. 2010; upa-data set; Appendix S3, Supporting Information). To enable comparisons between data sets, all taxonomies were tuned to the 8-rank taxonomy of the PR2 database, which follows the most up-to-date higher classification of eukaryotes (Adl et al. 2012). Accordingly, each OTU was assigned to a consensus taxonomy shared by at least 60% of its sequences. To assess phylogenetic diversity within samples from the different fertilization treatments, representative sequences for each OTU (i.e. the most abundant sequence in the OTU) from each data set were subjected to multiple alignment using MUSCLE (Edgar 2004). This analysis was not conducted for the fun-data set because the ITS region is not suitable for fungal phylogenetic reconstruction (Kr€ uger et al. 2012). The five sequence alignments were then used to construct bootstrapped maximum-likelihood phylogenetic trees using RAxML (Stamatakis 2006).

Statistical analyses All statistical analyses were conducted on community matrices of trimmed, subsampled, clustered, singletons and doubletons removed and taxonomic target-filtered data sets. The alpha-diversity of the samples was evaluated based on OTU richness and diversity (Shannon and Simpson indexes), and differences between fertilization treatments were analysed using the Mann–Whitney U-test. Rarefaction and accumulation analyses were performed to evaluate sequencing depth at the sample level and sampling depth at the experiment level, respectively. Beta-diversity based analyses were further conducted on relative abundance transformed community matrices. Permutational multivariate analysis of variances (PERMANOVA) based on the Bray–Curtis community dissimilarity index was performed to analyse the effect of fertilization treatment on community composition in each data set, using 999 permutations for each test. This analysis was also applied to the most abundant taxonomical groups in the euk-data set after normalizing the groups to equalize the number of sequences per sample. PERMANOVA was applied to 100 independent normalized group matrices, with p-values being corrected for multiple comparisons by controlling the false discovery rate using the method of Benjamini & Hochberg (1995). Dissimilarities in the phylogenetic compositions for the different treatments were evalu-

ated with the UniFrac unweighted algorithm (Lozupone & Knight 2005) in conjunction with maximum-likelihood trees and then used to construct a dissimilarity matrix for a PERMANOVA test. In addition, to assess the effects of fertilization on the different taxonomic groups of each data set, the organic and mineral fertilization treatments were fitted to redundancy analysis (RDA) ordinations of hellinger-transformed community matrices (Legendre & Gallagher 2001). These RDA analyses were applied separately to OTUs with common taxonomies and at different taxonomic levels (for taxonomic groups that were represented at least once in each sample and by at least seven OTUs). Significance was tested using the ‘envfit’ function of the R vegan package, using 999 permutations. All analyses other than the unweighted UniFrac distance calculations were performed using the R version 3.0.2 software (R Core Team 2013).

Co-occurrence and network analysis To compare community assemblages between the OTUs of the different data sets, a co-occurrence analysis was conducted on a merged presence/absence matrix of OTUs from the six data sets. A general co-occurrence profile for the whole OTU community was constructed by calculating C-scores and comparing them to a null model based on 1000 randomized matrices subject to a fixed rows and columns constraint (i.e. with constant row and column sums), using the procedure implemented in MOTHUR. Nonrandom co-occurring OTU pairs were then identified with the Pairs program (Ulrich 2008), using the same null model to randomize the presence/absence matrix and C-score as a pairwise measure of co-occurrence. To overcome computational limitations, the presence/absence matrix was broken down into a set of submatrices representing 70 OTUs each, and all submatrix pairs (each of which represented 140 OTUs) were analysed using Pairs. This procedure was repeated 10 times, and nonrandom cooccurrence pairs found in at least nine replicates were retained for further analyses. The nonrandom co-occurrence network was visualized using Cytoscape (Smoot et al. 2011). The preferential nonrandom co-occurrence branching of the OTUs (i.e. the proportion of edges connecting group x OTUs to group y OTUs) was assessed at the phylum level. This observed network topology was compared to 100 simulated network topologies in which the degree value for each node (i.e. for each OTU, the number of first-neighbours OTUs directly connected by a nonrandom co-occurence edge) was conserved, using the degree.sequence.game tool of the R package ‘IGRAPH’ (Csardi & Nepusz 2006). The significance of the between phyla branching preferences for each OTU was tested by performing a Mann–Whitney © 2014 John Wiley & Sons Ltd

S O I L E U K A R Y O T I C M I C R O O R G A N I S M S M E T A B A R C O D I N G 3345 tion in the cer- and fun-data sets, respectively. Community composition was also affected by mineral fertilization in the chry-data set, and by the interaction of mineral and organic fertilization in the kin- and upadata sets. In addition, when the euk-data set was splitted into subsets of OTUs sharing the same taxonomic assignment, organic fertilization was found to have significant effects on community composition for Fungi, Cercozoa and Ciliophora (Table S5, Supporting Information). A RDA analysis (Fig. S4, Supporting Information) showed that the soil chemical parameters (pH, moisture, available carbon and nitrogen, Table S6, Supporting Information) were significantly and positively correlated with the organic fertilization treatment (999 permutations, P < 0.05). Consequently, only the treatment parameter was conserved when studying the responses of the different eukaryotic microbial communities to environmental parameters. At the domain and kingdom levels, the RDA analysis (Fig. S5, Supporting Information) revealed a similar significant effect of organic fertilization on community composition to that seen in the PERMANOVA analysis. At deeper taxonomic levels, the taxonomic groups from the euk-data set that overlapped with other data sets exhibited variable responses to organic fertilization. All OTUs representing the Ascomycota orders Dothideomycetes and Sordariomycetes as well as Agaricomycetes (Basidiomycota) and Mortierella (Zygomycota) orders from both the euk- and fun-data sets were significantly affected by organic fertilization. However, while OTUs affiliated with almost all of the higher groups of the Cercozoa (Edomyxa, Filosa-Imbricatea, Filosa-Sarcomonadea, Filosa-Thecofilosea), and the four classes they represented were significantly affected by the organic fertilization in the cer-data set, organic fertilization did not affect cercozoan class OTUs in the euk-data set. Similarly, while Chlorophyta OTUs of the

U-test to compare the observed and simulated networks for each pair of phyla. The full analysis scripts are deposited in Dryad (Lentendu et al. 2014).

Results The influence of fertilization on diversity, community composition and the core microbiome

150

In 1 treatment In 2 treatments In 3 treatments In 4 treatments Core microbiome

100

ab ab a b

e ef ef f c

c cd d

N

M O OM

0

50

Number of OTUs

200

The total number of analysed sequences after trimming, subsampling and removal of singleton/doubleton OTUs ranged from 6.598 (kin) to 29.766 (euk), with an average of 550 (kin) to 2.481 (euk) sequences per sample (Table S4, Supporting Information). The eukaryotic microbial alpha-diversity ranged from 16 (kin M1) to 190 (euk O3) OTUs (Fig. 1) and was relatively well saturated (Fig. S2, Supporting Information). The OTU richness data revealed small variations among the fertilization treatments, with significantly more OTUs in the samples from the O (cer-data set) or OM treatments (kin- and upa-data sets). Such low variations were also observed for Shannon and Simpson diversity indexes (Fig. S3, Supporting Information). Interestingly, the majority of the OTUs, ranging from 46.5% for the euk-data set to 77% for the chry-data set, were represented in at least one sample from each fertilization treatment. The set of OTUs that were present in all samples for all treatments was defined as the core microbiome (i.e. OTUs present in all samples, see Shade & Handelsman 2012; Huse et al. 2012) and comprised 16.1–37.8% of each sample OTUs in the respective data sets (Fig. 1). Eukaryotic microbial beta-diversity responded strongly to the fertilization treatments. Organic fertilization significantly modified the community composition in all data sets except euk- (Table 1) and explained 46% and 48% of the dissimilarities in community composi-

N M O OM

cer

N

M O OM

chry

© 2014 John Wiley & Sons Ltd

N

M O OM

euk

N

M O OM

fun

kin

N

M O OM

upa

Fig. 1 Operational taxonomic units (OTU) richness in the different fertilization treatments for all six data sets. Bars are split according to the occurrence patterns of the corresponding OTUs. The standard deviation of the mean OTU richness in each treatment is shown above the bars. Different letters alongside the bars within each data set denote significant differences in OTU richness between fertilization treatments (Mann-Whitney test, P value < 0.05).

3346 G . L E N T E N D U E T A L . Table 1 Permutational multivariate analysis of variance exploring the differences in community composition (from Bray–Curtis dissimilarity matrices) in each data set separately for the different fertilization treatments cer

chry

euk

fun

kin

upa

Explanatory parameters

R2

P value

R2

P value

R2

P value

R2

P value

R2

P value

R2

P value

Organic* Mineral Organic & Mineral

0.46 0.06 0.08

0.001 0.21 0.126

0.40 0.20 0.08

0.001 0.01 0.167

0.13 0.04 0.11

0.173 0.954 0.31

0.48 0.10 0.08

0.001 0.07 0.122

0.37 0.07 0.17

0.001 0.216 0.022

0.26 0.12 0.14

0.001 0.024 0.015

*All samples with organic fertilization (O and OM) were tested against all samples without organic fertilization (M and N), and so on for the two other tests.

upa-data set were affected by both mineral and organic fertilization, no such effect was observed for Chlorophyta OTUs from the euk-data set. Additionally, Oxytrichidae (Ciliophora, euk-data set), Streptophyta (eukand upa-data sets), stramenopiles (Clade-c from chrydata set, euk-data set, Bacillariophyceae from upa-data set) and Neobodo/Rhynchobodo (kin-data sets) OTUs were found to be structured by organic fertilization. The only group that reacted significantly to mineral fertilization only was the Tubulinea (Amoebozoa, eukdata set). Organic fertilization was also observed to structure phylogenetic diversity. Between 28 and 59% of the phylogenetic diversity in the cer-, euk-, kin- and upa-data sets was explained by organic fertilization, while mineral fertilization explained no more than 10% of these dissimilarities (Table S7, Supporting Information).

Overlaps in the taxonomic compositions The taxonomic compositions for the six different data sets (Fig. 2) partially mirrored the differences between the four fertilization treatments previously identified by the beta-diversity analyses. However, taxonomic composition appeared more stable when only OTUs incidence was taken into account (Fig. S6, Supporting Information). The euk-data set was dominated by Streptophyta sequences (46 OTUs, 51.8% of reads, from which 93% were assigned to a Medicago sp. OTU), as well as animals (49 OTUs, 15.56% of reads) and fungal sequences (199 OTUs, 24.4% of reads). The primer pairs targeting protists groups were completely specific for their targeted groups (cer-, chry and kin-data sets). Similarly, only fungal sequences, in majority Ascomycota (224 OTUs, 59.1% of reads), were recovered in the fundata set. The upa-data set mainly contained Cyanobacteria (16 OTUs, 24% of reads), Chlorophyta (65 OTUs, 16.5% of reads) and stramenopiles (13 OTUs, 13.9% of reads). A total of 67 OTUs (23.2% of reads) corresponded to nontarget Bacteria and were removed from this last data set.

To compare the taxonomic coverage of the universal eukaryotic barcode marker to the five specific ones, the number of OTUs with a common taxonomic assignment among markers was calculated (Fig. 3). At the kingdom level (Fig. 3a), the euk-data set overlapped completely with the cer- (Rhizaria), chry- (stramenopiles) and fundata sets (Fungi/Ophistokonta) but exhibited only 86% overlap with the upa-data set (Archaeplastida, stramenopiles) because this primer pair also amplified Cyanobacteria. The euk-OTUs did not share any taxonomic information with the kin-OTUs because the eukdata set did not contain any Excavata sequences. The number of OTUs that could be taxonomically assigned decreased as the level of taxonomic identification went from phylum to species (Fig. 3, Table S8, Supporting Information), with a corresponding gradual decrease in the degree of taxonomic overlap between euk- and the other data sets. Thus, at the phylum level, 155 cerOTUs, 315 fun-OTUs and 29 upa-OTUs shared taxonomic assignments with 59, 160 and 30 euk-OTUs, respectively (Fig. 3b). At the genus level, 21 cer-OTUs, 53 fun-OTUs and 8 upa-OTUs shared taxonomic assignments with 19, 35 and 3 euk-OTUs, respectively (Fig. 3d). None of the chry-OTUs shared a taxonomic affiliation with a euk-OTU at the phylum level because none of the euk-OTUs were assigned to the Chrysophyceae/Synurophyceae. Overall, only 169 of the 471 eukOTUs shared a taxonomic assignment with 205 of the 760 OTUs from the specific primers data sets.

Co-occurrence patterns in the eukaryotic soil microbial community The co-occurrence C-score was significantly higher in the observed community matrix for the OTUs of all data sets than in randomized matrices (observed Cscore = 0.479, average randomized C-score = 0.476, SD = 3.5 E-4, P < 0.001), indicating the presence of nonrandom co-occurrence patterns among OTUs. Of the 759 528 possible OTU pairs, the PAIRS software identified 2010 pairs that co-occurred nonrandomly and 1627 that © 2014 John Wiley & Sons Ltd

S O I L E U K A R Y O T I C M I C R O O R G A N I S M S M E T A B A R C O D I N G 3347

O

OM

0.4 0.2

0.2 N

M

O

OM

(e)

M

O

O

OM

Eukaryota Archaeplastida Chlorophyta Chlorophyceae Trebouxiophyceae Ulvophyceae Streptophyta Embryophyceae Klebsormidiophyceae Stramenopiles PX_clade

0.6 0.4 0.2

Bacillariophyta Alveolata Dynophyta Dinophyceae Bacteria Cyanobacteria Chroococcales Nostocales Nostocaceae

0.0

0.2

OM

M

0.8

0.8 0.4

0.6

Kinetoplastida Eubodonid Neobodonid Dimastigella Neobodo Rhynchobodo Rhynchomonas Parabodonid Parabodo

0.0

0.0

0.2

0.4

0.6

0.8

Fungi Ascomycota Dothideomycetes Eurotiomycetes Leotiomycetes Sordariomycetes Basidiomycota Agaricomycetes Tremellomycetes Zygomycota Chytridiomycota Chytridiomycetes Glomeromycota Glomeromycetes

N

(f)

1.0

1.0

(d)

N

0.6

0.8

1.0 0.8 0.4

0.6

Stramenopiles Chrysophyceae Clade−C Ochromonas Pedospumella Poteriospumella Spumella

0.0

M

Eukaryota Stramenopiles Alveolata Ciliophora Amoebozoa Lobosa Archaeplastida Chlorophyta Streptophyta Fungi Ascomycota Basidiomycota Chytridiomycota Glomeromycota Mucoromycota Metazoa Annelida Arthropoda Nematoda Tardigrada Rhizaria Cercozoa

1.0

N

(c)

Leukarachnion

0.0

0.0

0.2

0.4

0.6

0.8

Cercozoa Filosa−Granofilosea Limnofilida Filosa−Imbricatea Euglyphida Spongomonadida Thaumatomonadida Filosa−Sarcomonadea Cercomonadida Glissomonadida Endomyxa Novel−clade−9 Vampyrellida Endomyxa−Phytomyxea Plasmodiophorida Filosa−Thecofilosea Cryomonadida Novel−clade−10−12 Tremulida

1.0

(b)

1.0

(a)

N

M

O

OM

N

M

O

OM

Fig. 2 The average taxonomic composition in each fertilization treatment was calculated based on the proportion of reads. OTUs gathering less than 0.7% of reads were merged to a higher taxonomic level. To enhance readability, OTUs assigned to a deeper taxonomic rank were merged into a common higher taxonomy as follows: Cercozoa OTUs are represented until the family level (a), Chrysophyceae OTUs until genus level (b), Eukaryota OTUs until class level (c), Fungi OTUs until order level (d), Kinetoplastida OTUs until genus level (e) and upa-data set OTUs until order level (f).

(a)

chry 56

(b)

upa 102

fun 359

cer 187

kin 46

euk 459

(c)

upa 86

cer 187 fun 345

euk 432

kin 46 chry 56

(d)

upa 27

cer 128 euk 162

fun 204

kin 43 chry 47

upa 7

cer 75 euk 152 fun 219

kin 25

chry 19

Fig. 3 Venn diagrams accounting for the overlap between the taxonomic assignment of the operational taxonomic units (OTUs) from the different data sets at the kingdom (a), phylum (b), family (c) and genus (d) levels, respectively. Taxonomic assignments are based on the PR2 (cer-, chry-, euk- and kin-data set), UNITE (fun-data set) and NCBI (upa-data set) taxonomies. The sizes of the circles and overlap regions correlate with the number of OTUs in the data set and the number of shared OTUs between two or more data sets, respectively. The numbers indicate the quantity of OTUs successfully identified at the displayed taxonomic level for each data set. Venn diagrams were generated with the venneuler R-package.

excluded themselves nonrandomly. Due to its ubiquity, no nonrandom co-occurrence or exclusion was observed for the core microbiome (Fig. S7, Supporting Information). Indeed, the null model used to randomize the presence/absence matrix conserved the row and column sum, making such ubiquitous OTUs insensitive to randomization. Nonrandom co-occurrence was only © 2014 John Wiley & Sons Ltd

observed for 583 OTUs occurring in three to 10 samples, and nonrandom exclusion was identified for 464 OTUs occurring in four to nine samples. The node degree values for both networks ranged from 1 to 28, with an average value of 4.5. OTUs with a degree of 10 or more were considered to exhibit high connectivity. Based on this criterion, three highly intra-connected

3348 G . L E N T E N D U E T A L . OTU groups were identified in the co-occurrence network (Fig. 4, Fig. S8, Supporting Information), while no comparable structure was found in the exclusion network (Fig. S9, Supporting Information). These groups contained 37, 32 and 11 OTUs, respectively, and occurred in 5.4 samples on average. Their occurrence was related to the applied fertilization regime: group (a) occurred almost exclusively in organic fertilized samples, group (b) was detected almost exclusively in nonorganic fertilized samples, and group (c) was detected only in nonfertilized samples (Table S9, Supporting Information). By comparing the observed network to simulated networks with identical topologies (i.e. in which the degree values for all nodes were conserved), we were able to identify some important branching preferences between OTUs of different phyla (Fig. 5). For example, Ascomycota OTUs co-occurred nonrandomly with Cercozoa OTUs more often than would be expected by chance alone (Mann–Whitney U-test, P value < 0.05), indicating that these groups have similar environmental preferences and that there may be direct or indirect interactions between the corresponding OTU pairs. Interestingly, such preferential branching was generally observed between phyla of organisms from different trophic levels; the only exception was the co-occurrence of OTUs from the autotrophic stramenopile PX-Clade and the Chlorophyta.

(a)

Discussion The multiple barcoding approach and taxonomic overlap Previously, characterization of soil eukaryotic communities by high-throughput pyrosequencing used universal eukaryotic barcode primers targeting the 18S rRNA gene (Urich et al. 2008; Meadow & Zabinski 2012; Baldwin et al. 2013; Bates et al. 2013). Here, we followed this strategy by using the euk-primer pair but paralleled it with a multiple barcoding approach to enhance the coverage of eukaryotic soil microorganisms. Likewise, 454 pyrosequencing was preferred over other NGS technics (e.g. Illumina), as it allowed the best compromise between a reliable description of microbial alpha- and beta-diversity due to its moderate sequencing depth (Lundin et al. 2012) and a better taxonomical resolution due to its long reads. The taxonomic overlap analysis confirmed our first hypothesis, as 71% of the taxa recovered by the specific primers were not detected by the universal eukaryotic barcode primer pair. This is consistent with the observations of Stoeck et al. (2006), who found that about half of the eukaryotic microbial diversity is missed when only one universal eukaryotic barcoding primer pair is used. The use of universal primers would have the advantage to more broadly compare samples and grasp the dominant microbial

(c)

(b)

Fig. 4 Subnetworks of operational taxonomic units (OTUs) (nodes) with high centrality (node degree > 9) in the initial nonrandom co-occurrence (edges) network representing the fertilization-selected microbiome. The OTUs were divided into three groups using the walktrap.community tool of the igraph R-package, one group occurred almost exclusively in nonorganic fertilized samples (a), one only in nonfertilized samples (b), and the third almost exclusively in organic fertilized samples (c). The node colours correspond to the different data sets (cer, orange; chry, dark yellow; euk, pink; fun, red; kin, green; upa, dark blue) and their sizes correspond to log-transformed OTU abundance (ranging from 3, p23_180 (b), to 1860 reads, euk_2(c)). Details on the taxonomic identification of OTUs are provided in the Table S9 (Supporting Information). © 2014 John Wiley & Sons Ltd

S O I L E U K A R Y O T I C M I C R O O R G A N I S M S M E T A B A R C O D I N G 3349

Fig. 5 Generalized co-occurrence and exclusion networks for the eukaryotic microbial phyla were generated by analysing the proportion and phyla of first-neighbour nodes for each operational taxonomic units (OTU) in the initial co-occurrence and exclusion networks (Fig. S8, Fig. S9, Supporting Information). They were compared to similar results generate from 100 randomly re-branched networks with conserved node degree values that were created using the degree.sequence.game tool of the igraph R package. Two phyla were connected if, in proportional terms, the OTUs of the first phylum co-occurred with (indicated by green, solid edges) or excluded (indicated by red, dashed edges) OTUs of the second phylum significantly more often than in the simulated network (Mann-Whitney test, P value < 0.05). The shapes and colours of the phyla nodes indicate their broad trophic status: primary producer (light green, diamond), predominantly predator (red, octagon), predominantly saprophyte (blue, triangle), and obligatory plant symbiont (purple, rounded rectangle).

community. However, such primers are limited to amplify rare taxa or taxa which have lower proportions of template DNA in a DNA extract. The approach we used in this article has the advantage to target specific groups regardless of their template concentration, thereby providing a better coverage of taxonomic groups as compared to universal primers. The presented method is thus suitable for analysing eukaryotic © 2014 John Wiley & Sons Ltd

microbial diversity in soil, although it could be further improved by adding another universal eukaryotic primer pair taking into account the variable primer amplification bias (Stoeck et al. 2006), as well as additional specific primers targeting different soil taxa, as recently introduced for ciliates (Lara et al. 2007b; Stoeck et al. 2014). We also want to point out that nontarget bacterial sequences can now be excluded from upa-amplicons by using a modified reverse primer (Hamsher et al. 2011). However, to improve the multibarcoding approach for protists, we need more primers specific for the eukaryotic supergroups and their major lineages. Future studies could use these primers for a PCR screening of samples prior to sequencing and may select the most promising target groups for the preparation and sequencing of libraries. Another way of increasing the proportion of microbial eukaryotic sequences identified when using a universal eukaryotic marker would be to use peptide nucleic acid-mediated PCR clamps (Von Wintzingerode et al. 2000) that specifically target plant, animal and fungal 18S rRNA gene to prevent their amplification during PCR. The new cercozoan primer pair that was designed to amplify a sequence length suitable for pyrosequencing proved to be specific for Cercozoa. According to the euk-data set, this was the second most abundant group of eukaryotic microorganisms in the studied soil after the Ascomycota, demonstrating their relative dominance in the soil food web. This again is consistent with previous results (Baldwin et al. 2013; Domonell et al. 2013). The taxonomic overlap analysis further demonstrated that the use of a universal eukaryotic barcode can cause a lack of precision when characterizing eukaryotic microorganisms for various reasons, such as primer amplification preferences, PCR competition, lack of taxonomic resolution (for further discussion see Stoeck et al. 2006). In our study, Kinetoplastida and Chrysophyceae/Synurophyceae were totally absent from the euk-data set, although they were successfully characterized based on the kin- and chry-data sets. Conversely, the euk-data set covered 21 fungal genera that were not contained in the fun-data set. This analysis also highlights the difficulty of comparing the taxonomic compositions of data sets generated using different molecular markers for which the corresponding databases have different levels of completeness. This was particularly obvious with the upa marker, for which the reference databases contained very few Xanthophyceae or Bacillariophyceae sequences. At the beta-diversity level, the apparent differences in community structure between organic and nonorganic fertilized soils identified using specific primers differed from those seen with the universal primers.

3350 G . L E N T E N D U E T A L . The data generated using the universal eukaryotic marker only revealed an effect of fertilization at the community level after the OTUs had been separated at the kingdom level (Table 1, Table S5, Supporting Information). Moreover, when fertilization effects were analysed at different taxonomic levels, significant changes in the OTU composition of groups such as Glissomonadida, Leotiomycetes and Chlorophyta were detected in the cer-, fun- and upa-data sets but not in the eukdata set, even though it also covered these groups (Fig. S4, Supporting Information). Such changes in beta-diversity patterns could be explained by the comparatively low coverage of each target group in the euk-data set: the coverage of individual target groups is inversely proportional to the number of groups targeted by the universal eukaryotic marker, thus decreasing the taxonomical overlap with the specific marker data sets at deeper taxonomical levels (Fig. 3) as previously discussed. The application of a rarefaction analysis to each group in the euk-data set prior to the PERMANOVA (Table S5, Supporting Information) further reduced the number of target sequences per subsample after the normalization of the initial reads (see section 2.3). Little or no control can thus be exerted over the sequencing depth of the different eukaryotic microbial groups targeted by the universal eukaryotic primer. However, such control is possible when using a primer pair that is specific for a single microbial group due to the equimolar pool of sample DNA in such cases. Analyses using the chry- and kin-specific primers revealed a high diversity of Chrysophyceae and Kinetoplastida in the soil samples, which had not previously been recognized. Only a few publications have discussed the presence of these organisms in soils (Ekelund et al. 2001; Lara et al. 2007a; Domonell et al. 2013), and there has only been one published cultureindependent diversity assessment of both groups in soils (Glaser et al. 2014). Moreover, these groups were not detected at all when using universal eukaryotic barcoding primers, indicating that either there is relatively little of their DNA in the soil metagenomic pool or the selected universal primer pair has low specificity for them. However, sequences from both groups were retrieved from all samples (other than sample 13-C for chry-data set), with similar levels of OTU richness, suggesting that they are both ubiquitous in soil communities. This clearly demonstrates the merits of using multiple specific primer pairs in parallel to characterize eukaryotic soil microorganisms because the strength of the signals generated by less abundant groups may be reduced by the presence of plant, animal and fungal sequences in the universal eukaryotic barcode data sets.

Fertilization drives changes in eukaryotic soil microbial communities The long-term fertilization treatments applied at the study site had relatively little impact on the richness of the detected eukaryotic microorganisms. We observed a slight increase in the OTU numbers for Cercozoa, Chrysophyceae/Synurophyceae and microalgae in the organic fertilized soils. The stable richness patterns observed in all six data sets were partially due to the relatively high proportion of OTUs belonging to the core microbiome (16.1–37.8% in each sample). For comparative purposes, previous analyses of pyrosequenced ITS fungal diversity in alpine meadows and eukaryotic 18S rRNA gene diversity in semi-arid floodplains provided no clear evidence of a core microbiome (less than 1% of OTUs; Lentendu et al. 2011; Baldwin et al. 2013). This was probably because the variety of soils, temporal variation in soil humidity and diversity of associated plants were all much greater in those investigations than in this work. The existence of a core microbiome in the studied environment can be explained by the relative similarity of the sampled soils due to the stability of the site’s historical soil management regime and the limited range of plants grown on the studied soil (Blair et al. 2006). Because plant roots exert strong selective effects on the microorganisms in their rhizosphere (Grayston et al. 1998), the core microbiome may have been associated with the plant that was being grown at the site, Medicago sativa, even though the rhizosphere and bulk soil were not separated during sampling. The core microbiome may also indicate that the diversity of eukaryotic microorganisms in agricultural soils is comparatively limited. This would be consistent with previous observations for algal populations in crop fields relative to those seen in the soils of other agro-ecosystems (Zancan et al. 2006), for grassland yeasts in intensely managed soils (Yurkov et al. 2012), and for soil protists in tilled agricultural soil compared with those in untilled soil (Adl et al. 2006). In general, based on the six data sets generated in this study, the eukaryotic microbial communities in the soils of the four fertilization treatments could be divided into three groups of OTUs: generalists (i.e. the core microbiome, Fig. 1), specialists (i.e. the fertilization-selected microbiome, Fig. 4, Table S9, Supporting Information) and OTUs with neutral occurrence patterns including rare OTUs (i.e. the “rare biosphere”). The microbiome associated with fertilization was responsible for the significant differences between the fertilization treatments with respect to both community compositions (Table 1, Table S5, Fig. S4, Supporting Information) and phylogenetic diversity (Table S7, Supporting Information). Organic fertilization with © 2014 John Wiley & Sons Ltd

S O I L E U K A R Y O T I C M I C R O O R G A N I S M S M E T A B A R C O D I N G 3351 farmyard manure seemed to be the strongest driver of changes in community composition; mineral fertilization had much less pronounced effects. This may reflect the greater increase in nutrient availability provided by manure treatments, which are known to primarily affect bacterial communities but have indirect effects on predatory protists (Frosteg ard et al. 1997). It was also suggested that both mineral (NPK) and organic (farmyard manure) fertilization affect the community composition of diatoms and testate amoebae (Heger et al. 2012), although no comparison to unfertilized soils was undertaken in this study. Notably, our results contradict those of previous studies on fungal communities in the surface layers of farmyard manure-amended soils. For example, (Lejon et al. 2007) did not detect significant differences within this taxonomic group when comparing manure-treated and untreated soils. This discrepancy may have been due to their use of F-ARISA, which provides a much lower resolution of fungal communities than amplicon sequencing. In addition, these authors investigated a short-term fertilization treatment (~10 years), whereas we studied an experiment that had been running for 112 years. It is important to note that each treatment with farmyard manure would add a new external community to the soil, and this may have affected the large changes in soil community composition and phylogenetic diversity that were associated with organic fertilization. However, the studied plots had not been treated with manure for 2 years prior to sampling, so any organisms in the manure that were not well adapted to the soil conditions should have died out before the soil samples were collected. Moreover, no significant changes in eukaryotic community structure were observed in a previous study of soils that had been fertilized with different organic carbon inputs, including manure, over a long period of time (Marschner 2003). This suggests that periodic inoculations of new microorganisms in this way may have little effect on the composition of the indigenous eukaryotic soil community. However, it would be of interest to test survival of microorganisms introduced to soil with farmyard manure in a future study.

Analysis of habitat preferences and species interactions based on co-occurrence and exclusion networks The structures of the nonrandom co-occurrence and exclusion networks provided further information on the fertilization-selected microbiome. The high degree nodes in the co-occurrence network were divided into three OTU groups corresponding to preferential occurrence in organic-, nonorganic or nonfertilized soils (Fig. 4, Table S9, Supporting Information). This demon© 2014 John Wiley & Sons Ltd

strated the existence of a connection between network structure and habitat preference. There was only one genus (Paracercomonas) common to multiple groups, reflecting the specialization of certain organisms towards the pH, moisture and nutrient availability conditions associated with the different fertilization treatments (Table S6, Supporting Information). The three groups contained members of 12, eight and 10 phyla, respectively. This emphasizes that most eukaryotic microbial groups demonstrate adaptation to specific abiotic and biotic conditions even within similar habitats. The analysis of the selected microbiome, as defined by the co-occurrence network, also provided information on the niche boundaries and life strategy types of eukaryotic microorganisms. For example, five of the six Kinetoplastida OTUs in the selected microbiome occurred almost exclusively in fertilized samples (Fig. 4, Table S9, Supporting Information). This suggests a preference for nutrient-rich soils that is commonly found in organisms adopting a K-type strategy. In contrary, the five arbuscular myccorrhiza Fungi in the selected microbiome only occurred in nonorganic fertilized soils. This could be due to a reduced need for the cultivated plants to form symbioses in the rhizosphere as they already were supplied with nutrients from the manure fertilization, as previously observed in another longterm field experiment (Oehl et al. 2004). The ecological interpretation of the co-occurrence and exclusion patterns of the OTUs was hampered by the lack of functional information associated with their taxonomical assignments, although this issue has recently been addressed for prokaryotes (Martiny et al. 2013). It should be noted that a unique nonrandom cooccurrence between two OTUs does not necessarily reflect a direct or indirect interaction. On the other hand, when OTUs are grouped at higher taxonomical levels, they can be associated to form a broad trophic status, even though eukaryotic microbial phyla may comprise various trophic groups (Adl & Gupta 2006). The assembly pattern of OTUs grouped at the phylum level reduces the complexity of the network (Fig. 5), thus enabling the extraction of some trends. The lack of exclusion or co-occurrence pattern between phyla at the same trophic levels indicates a balance between interphylum competition and habitat preference or niche overlap. On the other hand, the co-occurrence or exclusion of phyla between trophic levels probably indicates that their OTU members have similar or opposing habitat preferences, as well as direct relationships between their members. One potential direct relationship identified in the studied soils is the nonrandom co-occurrence of cercozoan Vampyrellida OTUs that are known to be fungivorous and algivorous (Hess et al. 2012) with seven fungal and six algal OTUs as well as four

3352 G . L E N T E N D U E T A L . Cercozoa OTUs. Although this analysis excludes OTUs that co-occur due to their similar habitat preferences, we cannot directly prove that Vampyrellida preys on those fungal and algal taxa. Moreover, certain pairs of phyla (e.g. Discoba and Chlorophyta) exhibited both significant exclusion and co-occurrence between their members, indicating that the biological processes governing OTU assembly, such as the boundaries of ecological niches, were also influential at taxonomic levels below the phylum, as recently deduced from a bacterial and fungal cultivation experiment (Lennon et al. 2012). Further studies on co-occurrence and exclusion patterns will thus be necessary to reveal general rules for eukaryotic microbial assembly in soil food webs.

Conclusion The level of soil eukaryotic microbial diversity in the Bad Lauchst€adt long-term fertilization experiment determined using a multiple barcoding approach was twice as high as that calculated using a single universal eukaryotic barcode marker. Although a significant proportion of the microbial community was identical in all treatments, the microbiomes associated with fertilization were distinct enough to identify significant differences in community composition and phylogenetic diversity between treatments. Organic fertilization (farmyard manure) appeared to be the strongest driver of community composition, while mineral fertilization did not induce many significant changes in community structure. In the analyses using the universal eukaryotic marker, these fertilization effects could only be detected after splitting the OTUs between the different taxonomic groups that were covered. The use of a connectivity parameter in the co-occurrence network also made it possible to precisely extract groups of OTUs sensitive to fertilization treatments, thus associating potential bio-marker OTUs with particular fertilization conditions. Conversely, the core OTUs that were identified may represent bio-markers of cropping soil planted with Medicago sativa. Further studies (e.g. on other crop soils) and meta-analyses are needed to divide this core microbiome into soil generalists and OTUs that are specifically associated with the cultivated plant. The results presented herein provide a framework for the characterization of soil eukaryotic microbial diversity as well as their habitat preferences, niche spaces and potential interactions within the soil food web.

Acknowledgements This study was financed through the HIGRADE PhD scholarship of Guillaume Lentendu, and by special funds of the rectorate of the University of Leipzig. We thank Ines Merbach

and the workers of the Bad Lauchst€ adt Long-Term Fertilization Experiment for their support during soil sampling, and Elke Schultz and Gabriele Henning for kindly providing the soil parameters. We are extremely grateful to Sigrid H€ artling and Beatrix Schnabel for their technical support during amplicon library preparations and sequencing. We finally thank three anonymous reviewers for their helpful comments.

References Abarenkov K, Nilsson RH, Larsson K-H et al. (2010) The UNITE database for molecular identification of fungi – recent updates and future perspectives. New Phytologist, 186, 281– 285. Adl MS, Gupta VS (2006) Protists in soil ecology and forest nutrient cycling. Canadian Journal of Forest Research, 36, 1805– 1817. Adl SM, Coleman DC, Read F (2006) Slow recovery of soil biodiversity in sandy loam soils of Georgia after 25 years of notillage management. Agriculture, Ecosystems & Environment, 114, 323–334. Adl SM, Simpson AGB, Lane CE et al. (2012) The revised classification of eukaryotes. Journal of Eukaryotic Microbiology, 59, 429–514. Baldwin DS, Colloff MJ, Rees GN et al. (2013) Impacts of inundation and drought on eukaryote biodiversity in semi-arid floodplain soils. Molecular Ecology, 22, 1746–1758. Bates ST, Clemente JC, Flores GE et al. (2013) Global biogeography of highly diverse protistan communities in soil. The ISME journal, 7, 652–659. Behnke A, Engel M, Christen R et al. (2011) Depicting more accurate pictures of protistan community complexity using pyrosequencing of hypervariable SSU rRNA gene regions: protistan community complexity. Environmental Microbiology, 13, 340–349. Benjamini Y, Hochberg Y (1995) Controlling the false discovery rate – a practical and powerful approach. Journal of the Royal Statistical Society. Series B (Methodological), 57, 289–300. B
erard A, Dorigo U, Humbert JF, Martin-Laurent F (2005) Microalgae community structure analysis based on 18S rDNA amplification from DNA extracted directly from soil as a potential soil bioindicator. Agronomy for Sustainable Development, 25, 285–291. Bik HM, Porazinska DL, Creer S et al. (2012) Sequencing our way towards understanding global eukaryotic biodiversity. Trends in Ecology & Evolution, 27, 233–243. Blair N, Faulkner RD, Till AR, Korschens M, Schulz E (2006) Long-term management impacts on soil C, N and physical fertility: Part II: bad Lauchstadt static and extreme FYM experiments. Soil and Tillage Research, 91, 39–47. Boenigk J, Pfandl K, Stadler P, Chatzinotas A (2005) High diversity of the “Spumella-like” flagellates: an investigation based on the SSU rRNA gene sequences of isolates from habitats located in six different geographic regions. Environmental Microbiology, 7, 685–697. B€ ohme L, Langer U, B€ ohme F (2005) Microbial biomass, enzyme activities and microbial community structure in two European long-term field experiments. Agriculture, Ecosystems & Environment, 109, 141–152. Bonkowski M (2004) Protozoa and plant growth: the microbial loop in soil revisited. New Phytologist, 162, 617–631. © 2014 John Wiley & Sons Ltd

S O I L E U K A R Y O T I C M I C R O O R G A N I S M S M E T A B A R C O D I N G 3353 Botha A (2011) The importance and ecology of yeasts in soil. Soil Biology and Biochemistry, 43, 1–8. Cheung MK, Au CH, Chu KH, Kwan HS, Wong CK (2010) Composition and genetic diversity of picoeukaryotes in subtropical coastal waters as revealed by 454 pyrosequencing. The ISME journal, 4, 1053–1059. Coleman DC, Whitman WB (2005) Linking species richness, biodiversity and ecosystem function in soil systems. Pedobiologia, 49, 479–497. Core Team R (2013) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. Csardi G, Nepusz T (2006) The igraph software package for complex network research. InterJournal, Complex Systems, 1695. Domonell A, Brabender M, Nitsche F, Bonkowski M, Arndt H (2013) Community structure of cultivable protists in different grassland and forest soils of Thuringia. Pedobiologia, 56, 1–7. Edgar RC (2004) MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Research, 32, 1792–1797. Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R (2011) UCHIME improves sensitivity and speed of chimera detection. Bioinformatics, 27, 2194–2200. Edmeades DC (2003) The long-term effects of manures and fertilisers on soil productivity and quality: a review. Nutrient Cycling in Agroecosystems, 66, 165–180. Eisenhauer N, Beßler H, Engels C et al. (2010) Plant diversity effects on soil microorganisms support the singular hypothesis. Ecology, 91, 485–496. Ekelund F, Rønn R, Griffiths BS (2001) Quantitative estimation of flagellate community structure and diversity in soil samples. Protist, 152, 301–314. Esteban GF, Clarke KJ, Olmo JL, Finlay BJ (2006) Soil protozoa —an intensive study of population dynamics and community structure in an upland grassland. Applied Soil Ecology, 33, 137–151. Foissner W (1999) Soil protozoa as bioindicators: pros and cons, methods, diversity, representative examples. Agriculture, Ecosystems & Environment, 74, 95–112.  Petersen SO, B Frosteg ard A, a ath E, Nielsen TH (1997) Dynamics of a microbial community associated with manure hot spots as revealed by phospholipid fatty acid analyses. Applied and Environmental Microbiology, 63, 2224–2231. Fuhrman JA (2009) Microbial community structure and its functional implications. Nature, 459, 193–199. Glaser K, Kuppardt A, Krohn S et al. (2014) Primer pairs for the specific environmental detection and T-RFLP analysis of the ubiquitous flagellate taxa Chrysophyceae and Kinetoplastea. Journal of Microbiological Methods, 100, 8–16. Grayston SJ, Wang S, Campbell CD, Edwards AC (1998) Selective influence of plant species on microbial diversity in the rhizosphere. Soil Biology and Biochemistry, 30, 369–378. Guillou L, Bachar D, Audic S et al. (2012) The Protist ribosomal reference database (PR2): a catalog of unicellular eukaryote small sub-unit rRNA sequences with curated taxonomy. Nucleic Acids Research, 41, D597–D604. Hamsher SE, Evans KM, Mann DG, Poul
ıckov
a A, Saunders GW (2011) Barcoding diatoms: exploring alternatives to COI5P. Protist, 162, 405–422. Hedlund K, Santa Regina I, Van der Putten WH et al. (2003) Plant species diversity, plant biomass and responses of the

© 2014 John Wiley & Sons Ltd

soil community on abandoned land across Europe: idiosyncracy or above-belowground time lags. Oikos, 103, 45–58. Heger TJ, Straub F, Mitchell EAD (2012) Impact of farming practices on soil diatoms and testate amoebae: a pilot study in the DOK-trial at Therwil, Switzerland. European Journal of Soil Biology, 49, 31–36. Hess S, Sausen N, Melkonian M (2012) Shedding light on vampires: the phylogeny of vampyrellid amoebae revisited. PLoS ONE, 7, e31165. Horner-Devine MC, Silver JM, Leibold MA et al. (2007) A comparison of taxon co-occurrence patterns for macro- and microorganisms. Ecology, 88, 1345–1353. Houot S, Chaussod R (1995) Impact of agricultural practices on the size and activity of the microbial biomass in a long-term field experiment. Biology and Fertility of Soils, 19, 309–316. Huse S, Welch D, Morrison H, Sogin M (2010) Ironing out the wrinkles in the rare biosphere through improved OTU clustering. Environmental Microbiology, 12, 1889–1898. Huse SM, Ye Y, Zhou Y, Fodor AA (2012) A core human microbiome as viewed through 16S rRNA sequence clusters. PLoS ONE, 7, e34242. Kim M, Morrison M, Yu Z (2011) Evaluation of different partial 16S rRNA gene sequence regions for phylogenetic analysis of microbiomes. Journal of Microbiological Methods, 84, 81–87. Kr€ uger D, Kapturska D, Fischer C, Daniel R, Wubet T (2012) Diversity measures in environmental sequences are highly dependent on alignment quality—data from ITS and New LSU primers targeting basidiomycetes (JE Stajich, Ed,). PLoS ONE, 7, e32139. Lara E, Berney C, Ekelund F, Harms H, Chatzinotas A (2007a) Molecular comparison of cultivable protozoa from a pristine and a polycyclic aromatic hydrocarbon polluted site. Soil Biology and Biochemistry, 39, 139–148. Lara E, Berney C, Harms H, Chatzinotas A (2007b) Cultivationindependent analysis reveals a shift in ciliate 18S rRNA gene diversity in a polycyclic aromatic hydrocarbon-polluted soil. Fems Microbiology Ecology, 62, 365–373. Lawley B, Ripley S, Bridge P, Convey P (2004) Molecular analysis of geographic patterns of eukaryotic diversity in Antarctic soils. Applied and Environmental Microbiology, 70, 5963–5972. Legendre P, Gallagher ED (2001) Ecologically meaningful transformations for ordination of species data. Oecologia, 129, 271–280. Leinonen R, Akhtar R, Birney E et al. (2010) The European nucleotide archive. Nucleic Acids Research, 39, D28–D31. Lejon DPH, Sebastia J, Lamy I, Chaussod R, Ranjard L (2007) Relationships between soil organic status and microbial community density and genetic structure in two agricultural soils submitted to various types of organic management. Microbial Ecology, 53, 650–663. Lennon JT, Aanderud ZT, Lehmkuhl BK, Schoolmaster DR (2012) Mapping the niche space of soil microorganisms using taxonomy and traits. Ecology, 93, 1867–1879. Lentendu G, Zinger L, Manel S et al. (2011) Assessment of soil fungal diversity in different alpine tundra habitats by means of pyrosequencing. Fungal Diversity, 49, 113–123. Lentendu G, Wubet T, Chatzinotas A et al. (2014) Data from: effects of long term differential fertilization on eukaryotic microbial communities in an arable soil: a multiple barcoding approach. Dryad Digital Repository. doi: 10.5061/dryad.c4s8s.

3354 G . L E N T E N D U E T A L . Lozupone C, Knight R (2005) UniFrac: a new phylogenetic method for comparing microbial communities. Applied and Environmental Microbiology, 71, 8228–8235. Lundin D, Severin I, Logue JB et al. (2012) Which sequencing depth is sufficient to describe patterns in bacterial a- and bdiversity? Environmental Microbiology Reports, 4, 367–372. Marschner P (2003) Structure and function of the soil microbial community in a long-term fertilizer experiment. Soil Biology and Biochemistry, 35, 453–461. Martiny AC, Treseder K, Pusch G (2013) Phylogenetic conservatism of functional traits in microorganisms. The ISME Journal, 7, 830–838. Meadow JF, Zabinski CA (2012) Spatial heterogeneity of eukaryotic microbial communities in an unstudied geothermal diatomaceous biological soil crust: Yellowstone National Park, WY, USA. FEMS Microbiology Ecology, 82, 182–191. Merbach I, Schulz E (2012) Long-term fertilization effects on crop yields, soil fertility and sustainability in the Static Fertilization Experiment Bad Lauchst€adt under climatic conditions 2001–2010. Archives of Agronomy and Soil Science, 59, 1041– 1057. Nolte V, Pandey RV, Jost S et al. (2010) Contrasting seasonal niche separation between rare and abundant taxa conceals the extent of protist diversity: high seasonal protist abundance turnover. Molecular Ecology, 19, 2908–2915. Oehl F, Sieverding E, M€ader P et al. (2004) Impact of long-term conventional and organic farming on the diversity of arbuscular mycorrhizal fungi. Oecologia, 138, 574–583. Reeder J, Knight R (2009) The “rare biosphere”: a reality check. Nature Methods, 6, 636–637. Schloss PD, Westcott SL, Ryabin T et al. (2009) Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Applied and Environmental Microbiology, 75, 7537–7541. Shade A, Handelsman J (2012) Beyond the Venn diagram: the hunt for a core microbiome. Environmental Microbiology, 14, 4–12. Smoot ME, Ono K, Ruscheinski J, Wang P-L, Ideker T (2011) Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics, 27, 431–432. Stamatakis A (2006) RAxML-VI-HPC: maximum likelihoodbased phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics, 22, 2688–2690. Steele JA, Countway PD, Xia L et al. (2011) Marine bacterial, archaeal and protistan association networks reveal ecological linkages. The ISME journal, 5, 1414–1425. Stoeck T, Hayward B, Taylor GT, Varela R, Epstein SS (2006) A multiple PCR-primer approach to access the microeukaryotic diversity in environmental samples. Protist, 157, 31–43. Stoeck T, Behnke A, Christen R et al. (2009) Massively parallel tag sequencing reveals the complexity of anaerobic marine protistan communities. BMC Biology, 7, 72. Stoeck T, Breiner H-W, Filker S et al. (2014) A morphogenetic survey on ciliate plankton from a mountain lake pinpoints the necessity of lineage-specific barcode markers in microbial ecology. Environmental Microbiology, 16, 430–444. Sugiyama A, Vivanco JM, Jayanty SS, Manter DK (2010) Pyrosequencing assessment of soil microbial communities in organic and conventional potato farms. Plant Disease, 94, 1329–1335.

Ulrich W (2008) Pairs – A FORTRAN Program for Studying Pairise Species Associations in Ecological Matrices.Available from www.keib.umk.pl/pairs/ Urich T, Lanz
en A, Qi J et al. (2008) Simultaneous assessment of soil microbial community structure and function through analysis of the meta-transcriptome (N Ward, Ed,). PLoS ONE, 3, e2527. Van Dongen S (2000) Graph Clustering by Flow Simulation. PhD thesis, University of Utrecht, The Netherlands. Von der Heyden S, Cavalier-Smith T (2005) Culturing and environmental DNA sequencing uncover hidden kinetoplastid biodiversity and a major marine clade within ancestrally freshwater Neobodo designis. International Journal of Systematic and Evolutionary Microbiology, 55, 2605–2621. Von Wintzingerode F, Landt O, Ehrlich A, Gobel UB (2000) Peptide nucleic acid-mediated PCR clamping as a useful supplement in the determination of microbial diversity. Applied and Environmental Microbiology, 66, 549–557. Yurkov AM, Kemler M, Begerow D (2012) Assessment of yeast diversity in soils under different management regimes. Fungal Ecology, 5, 24–35. Zancan S, Trevisan R, Paoletti MG (2006) Soil algae composition under different agro-ecosystems in North-Eastern Italy. Agriculture, Ecosystems & Environment, 112, 1–12.

G.L., C.W., M.S. and F.B. designed the experiment. G.L. performed soil sampling. G.L. designed the new cercozoan primer pair. G.L., A.C. and T.W. developed the experimental molecular methodologies. G.L. and T.W. performed pyrosequencing and conducted bioinformatic and statistical analysis. G.L. wrote the paper with advice from all co-authors.

Data accessibility Raw pyrosequencing reads and metadata are available at the European Nucleotide Archive: http://www.ebi. ac.uk/ena/data/view/PRJEB5170. Quality trimmed and subsampled reads: Dryad doi:10. 5061/dryad.c4s8s OTU community matrices and taxonomic assignment: Dryad doi:10.5061/dryad.c4s8s Bioinformatic, statistical and network analyses scripts: Dryad doi:10.5061/dryad.c4s8s

Supporting information Additional supporting information may be found in the online version of this article. Appendix S1 New cercozoan-specific primer pair suitable for 454 pyrosequencing. Appendix S2 PCR conditions. Appendix S3 Threshold determination.

© 2014 John Wiley & Sons Ltd

S O I L E U K A R Y O T I C M I C R O O R G A N I S M S M E T A B A R C O D I N G 3355

Table S1 Fertilization treatments in the Bad Lauchst€adt longterm fertilization experiment (modified in 1995). Table S2 Primer pairs. Table S3 In silico amplification performances of cercozoan-specific primer pairs. Table S4 Amount of sequences in the different data sets. Table S5 Permutational multivariate analysis of variance exploring the differences in community composition (from Bray–Curtis dissimilarity matrices) depending on different fertilization treatments in the main taxonomical groups of the euk-data set, subsampled at same amount of sequences per samples. Table S6 Mean value (and standard deviation) of physicochemical parameters in fertilization treatments.

Fig. S2 (a) Rarefaction curves were calculated at the minimal sequencing depth among samples belonging to the same data sets after the removal of singleton and doubleton OTUs. (b) Accumulation curves show the OTU discovery rate with samples added randomly one by one. Fig. S3 OTU diversity in the different fertilization treatments for all six data sets calculated with the Shannon (a) and Simpson (b) diversity indexes. The standard deviation of the mean OTU richness in each treatment is shown above the bars. Fig. S4 Redundancy analysis (RDA) ordinations were fitted to environmental variables for all six data sets ((a) cer, (b) chry, (c) euk, (d) fun, (e) kin, (f) upa). Fig. S5 Redundancy analyses (RDA) were applied to the main taxonomic groups in each data set, gathering at least seven OTUs and with at least one OTU occurring in each sample, using hellinger-transformed community matrices.

Table S7 Permutational multivariate analysis of variance exploring the differences in phylogenetic composition, based on Unifrac unweighted distances, depending on different fertilization treatments.

Fig. S6 The average taxonomic composition in each fertilization treatment was calculated in terms of the corresponding proportion of OTUs.

Table S8 Per cent of OTUs identified at the different taxonomic levels.

Fig. S7 OTU richness was plotted in terms of increasing incidence for each data set separately ((a) cer, (b) chry, (c) euk, (d) fun, (e) kin, (f) upa).

Table S9 Taxonomic identification and occurrences of OTUs from the selected microbiome (degree > 9) extracted from the nonrandom co-occurrence network. Fig. S1 Intraspecific vs. interspecific genetic distances were used as proxies for the in silico determination of clustering thresholds for the OTUs.

© 2014 John Wiley & Sons Ltd

Fig. S8 Network showing the OTUs (nodes) with nonrandom co-occurrence (edges) as defined by the Pairs software package. Fig. S9 Network showing the OTUs (nodes) with non-random exclusion (edges) as defined by the Pairs software package.