Belowground fungal community diversity, composition and ecological functionality associated with winter wheat in conventional and organic agricultural systems

Study sites

During the end of April 2015, winter wheat (Triticum aestivum) plants were sampled from 10 different sites located in the surroundings of Tübingen, Baden-Württemberg, Germany (see Fig. 1). We sampled agricultural fields treated under conventional farming management located in Pfrondorf (C1, C2 and C4) and Entringen (C3 and C5) as well as fields under organic farming management located in Unterjesingen (O1, O2, O3, O4 and O5). Within each agricultural field, plants were collected along a diagonal transect. Organic farming sites most commonly featured frequent and diverse co-occurring herbaceous plants. Site O2 was also surrounded by shrubs and site O3 was located next to a group of ectomycorrhizal (ECM) trees. For the conventional farming sites, sites C2 and C3 had more diverse co-occurring plant species. A complete list of co-inhabiting plants at each collection site is given as Supplemental Information (Table S1). The farmers Ulrich Bechtle (Unterjesingen) and Eckhart Wizemann (Entringen, Pfrondorf) allowed us to collect plant specimens from their crop fields.

Plant and soil sample collection

Plants from a sampling transect, defined as the longest diagonal between two of the corners of each site, were sampled. Along each transect, 10 plots (labelled A–J) were positioned at similar distance intervals, with the two outer plots (A and J) being placed directly in the corners, the distance between plots ranged from 8 m to 28 m (Table S1). Within each sampling plot, at least five different plants at the tillering stage were randomly chosen and carefully dug out of the ground with a shovel, retaining as much of the root system as possible and a portion of soil. Plant size ranged from approximately 15 cm to 30 cm; individuals growing in the corners of the sites were always smaller than those growing in the central area. Sampled plants were stored in plastic bags at 4 °C until further processing.

The root systems of each plant were soaked in tap water to loosen the soil and were then washed extensively under running tap water. Subsequently, the five plants per plot with the highest number of healthy roots were chosen for molecular analysis. The roots were washed carefully under running purified water, then three fine roots per plant were randomly selected, cut off with a sterile scalpel and forceps, and were placed into a Petri dish filled with double-distilled and sterilised water. Root surface sterilisation was performed via a modified version of the method of Sun et al. (2008), which consisted of rinsing the roots in ethanol (70%) for 5 min in a Petri dish, then in sodium hypochlorite (NaOCl 1%) for another 5 min and finally washing them with double-distilled water three times to remove any remaining chemicals. After that procedure, all 15 root portions from each plot (three root segments from five plants each) were pooled together as one sample, placed into 2-mL Eppendorf tubes and dried at 50 °C in a Dörrex device for several days. In total, the roots of 500 plants were randomly sampled (representing 100 pooled samples), 250 plants for each of the two different farming methods. Herbarised plants were deposited in the Herbarium Tubingense (TUB) as Nos. TUB 021591 to TUB 021628.

In addition, from each plot, portions of soil adhering to the roots of the sampled plants were removed and stored in paper bags for drying. Two soil samples per site were analysed independently: soil from the two edge plots was mixed and pooled as one sample, whereas the soil from the remaining middle plots was used as a second sample. Measurements of total carbon and nitrogen content as well as organic carbon content, pH and plant-available phosphorus content were conducted at the Soil Science and Geoecology Laboratory of the University of Tübingen. Additionally, the carbon-to-nitrogen ratio was calculated and the humus content was estimated by multiplying the organic carbon content by 1.72 as described by Ad-hoc Arbeitsgruppe Boden (2005).

DNA extraction, amplification, cloning and DNA sequencing

Dried fine root pools were homogenised with a mixer mill (MM 300; Retsch GmbH, Haan, Germany) with two steel beads per Eppendorf tube at a frequency of 30 Hz for 10 min. Total gDNA was extracted with the InnuPREP Plant DNA Kit (Analytik Jena AG, Jena, Germany) following the manufacturer’s protocol, except that the amounts of lysis buffer, RNase, protein kinase and binding solution were doubled in order to dissolve the large volume of homogenised root material completely. To detect root-associated fungi, PCRs were performed to amplify the ITS region by following two different approaches: (i) the first approach was used to amplify the ITS region including the 5.8S domain with the fungal-specific primer ITS1F (Gardes & Bruns, 1993) and the universal primer ITS4 (White et al., 1990); (ii) in the second approach, the ITS region and the large subunit domain D1/D2 were amplified with the fungal-specific primers NSI1 (Martin & Rygiewicz, 2005) and NL4 (White et al., 1990). PCR products from the second approach (for samples producing negative or weak PCR products) were subsequently used for nested PCRs with the primer combination ITS1F–ITS4. In both cases, MangoTaq™ DNA Polymerase (Bioline GmbH, Lückenwalde, Germany) was used for amplification and the reactions consisted of a volume of 25 µL containing 5 µL coloured reaction buffer, 14.5 µL double-distilled H₂O, 0.75 µL MgCl₂ (50 mM), 1 µL dNTPs (5 mM), 0.5 µL of the forward primer (25 pmol/µL), 0.5 µL of the reverse primer (25 pmol/µL), 0.25 µL MangoTaq DNA Polymerase (2 unit/µL) and 2.5 µL of the DNA template. For direct PCRs with the NSI1–NL4 primers, PCR reactions consisted of half the volume. PCRs were run with the thermo programme TOUCH 60 under the following conditions: initial denaturation at 94 °C for 3 min; 10 cycles each of denaturation at 94 °C for 30 s, annealing at 60 °C for 45 s, with each cycle decreasing by 1 °C, and elongation at 72 °C for 1 min 15 s, followed by 26 cycles each of denaturation at 94 °C for 30 s, annealing at 50 °C for 45 s, elongation at 72 °C for 1 min 15 s and a final extension phase at 72 °C for 7 min. PCR products were verified by running an agarose gel electrophoresis with 0.7% agarose and 1× Tris-Borat-EDTA (TBE) buffer for 15 min at 143 V. The PCR products were subsequently stained with Midori Green (NIPPON Genetics EUROPE GmbH, Dueren, Germany) for 30 min and visualised under UV illumination.

Positive PCR products that could not be sequenced directly were cloned into competent Escherichia coli cells with the TOPO TA Cloning® Kit for Sequencing (Invitrogen, Life Technologies GmbH, Darmstadt, Germany). The cloning reaction for one sample contained 0.33 µL of a salt solution, 0.33 µL double distilled H₂O, 0.33 µL of pCR®4-TOPO vector (Invitrogen, Carlsbad, CA, USA) and 1 µL of the undiluted PCR product. Samples were incubated at room temperature for 30 min. Afterwards, 19 µL of One Shot® TOP10 chemically competent E. coli cells (Invitrogen, Carlsbad, CA, USA) was carefully added to the cloning reaction and the mixture subsequently incubated on ice for 30 min. Heat shock was conducted in a water bath at 42 °C for 30 s, after which the cells were kept on ice for 2 min. Finally, after adding 80 µL of Super Optimal broth Catabolite repression (S.O.C.) medium, the mixture was placed into an incubator at 37 °C and 200 rpm for 1, then plated onto Lysogeny Broth (LB)-medium including Kanamycin (1 mL/L). Plates were incubated at 37 °C overnight. Up to 12 bacterial colonies per sample were picked and directly used for PCRs with the forward and reverse primers M13F/M13R (Invitrogen, Carlsbad, CA, USA) and MangoTaq DNA Polymerase in 25-µL reactions. Amplification in the thermo programme M13 consisted of initial denaturation at 94 °C for 7 min, followed by 32 cycles each of denaturation at 94 °C for 20 s, annealing at 55 °C for 40 s, elongation at 72 °C for 2 min 40 s and a final extension at 72 °C for 7 min. For samples with a large number of positive PCR products, eight clones were digested enzymatically with Hinf1 (New England BioLabs GmbH, Frankfurt am Main, Germany) to detect restriction fragment length polymorphisms. To check for different restriction patterns, gel electrophoresis was run with 1% agarose and 1× Tris-EDTA (TE) buffer at 120 V for 40 min, followed by staining with Midori Green (Nippon Genetics EUROPE GmbH, Dueren, Germany) for 30 min and visualisation under UV illumination. All PCR products showing different restriction patterns were sequenced, as well as some with the same pattern.

Amplified PCR products were cleaned up with diluted (1:20) ExoSAP-IT® (Affymetrix UK Ltd., High Wycombe, UK) following the manufacturer’s protocol. Cycle sequencing was performed on 4 µL of a 1:6 diluted ABI PRISM BigDye® Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Life Technologies GmbH, Darmstadt, Germany) and 1 µL of the primers ITS1F–ITS4 and, if necessary, 5.8SR (Vilgalys & Hester, 1990), ITS2 or ITS1 (White et al., 1990) (12 pmol/µL) for directly sequenced PCR products and the primers M13F–M13R (12 pmol/µL) for cloned PCR products. The thermo programme ABI31CYCL was used for amplification under following conditions: initial denaturation at 96 °C for 1 min, followed by 27 cycles each of denaturation at 96 °C for 10 s, annealing at 50 °C for 5 s and elongation at 60 °C for 4 min. Afterwards, the cycle sequencing products were precipitated in 75% isopropanol, washed with 80% ethanol and dried in a vacuum centrifuge for approximately 5 min. DNA pellets were dissolved in 15 µL Hi-Di™ formamide (Applied Biosystems, Life Technologies, Warrington, UK) for 10 min and sequencing was conducted on an ABI3130xl Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). A list of primers used in this study, including their sequences, is given in the Supplemental Information (Table S2).

Sequence editing, chimaera checking and BLAST search

Forward and reverse ITS sequences were assembled automatically and edited manually by Sequencher v.4.9 (Gene Codes Corporation, Ann Arbor, MI, USA).

All sequences were first checked for PCR chimaeras with VSEARCH v2.4.3 (Rognes et al., 2016), implementing the UCHIME algorithm (Edgar et al., 2011). All chimeric and ‘borderline’ chimeric sequences were excluded from further analyses. Non-chimeric ITS sequences with GenBank accession nos. KY430446 to KY430584 were used as queries in a search against the GenBank ITS RefSeq Database (BioProject PRJNA177353, accessed 12 December 2018) with the programme MALT (Herbig et al., 2016), using the default parameters. The ITS sequences generated in this study were then annotated with the most specific taxonomic label inferred.

Fungal diversity estimation

Fungal ITS sequences were clustered into OTUs based on a ≥97% similarity threshold with the UCLUST algorithm implemented in VSEARCH (Edgar, 2010; Rognes et al., 2016). We estimated the alpha diversity of those samples by calculating Hill numbers in the R package ‘vegan’ (Oksanen et al., 2016). We also computed the number of OTUs that were shared by the two farming systems and the number of OTUs that were unique to a specific system. A pairwise species distance matrix was computed by first calculating the pairwise sequence distances via the alignment-free approach of spaced-words frequencies (Leimeister et al., 2014). Pairwise sequence distances were subsequently averaged to pairwise species (OTU) distances.

Phylogenetic community structure

We calculated the standardised effect sizes of the mean phylogenetic distance (MPD, equivalent to—net relatedness index (NRI), Webb et al., 2002) and the mean nearest phylogenetic taxon distance (MNTD, equivalent to—nearest taxon index (NTI)) metrics as implemented in the R package ‘picante’ to assess the composition of fungal samples in conventional vs. organically managed farming sites (Kembel et al., 2010). We compared the MPD and MNTD values for our samples with random samples generated under a null model that shuffled species labels randomly across the complete phylogeny.

Differences in environmental factors between farming methods

The differences in measured soil parameters between organic and conventional plots were visualised via non-metric multidimensional scaling (NMDS) using the Bray-Curtis index with functions implemented in the ‘vegan’ package (Oksanen et al., 2016) for R (R Development Core Team, 2011). Specifically, NMDS was performed to analyse the effect of, for example plant-available phosphorus, carbon-to-nitrogen ratio, pH value (pH CaCl₂) and humus content (Table S3).

We conducted a transformation-based Redundancy Analysis (tb-RDA) to assess the effect of the environmental parameters on the OTU distribution (Legendre & Legendre, 2012). Species abundances were transformed using Hellinger transformation (Legendre & Gallagher, 2001) and variables were selected using backward selection using functions from the ‘vegan’ package. Additionally, a PERMANOVA was conducted (using the Bray–Curtis index) to test whether the farming system has a significant influence on the OTU distribution (Anderson, 2001).

Functional guilds

We classified the identified taxa into ecological guilds with the FUNGuild tool and database (the database is available at www.stbates.org/funguild_db.php, accessed 23 January 2020) (Nguyen et al., 2016), to determine if specific functional groups of fungi (e.g. pathotrophs, saprotrophs, symbiotrophs, etc.) differed between sites and/or the farming methods applied. All OTUs that did not match any taxa in the database were categorised as ‘unassigned’.

PCR success and belowground fungal diversity

A total of 54% of the pooled root samples yielded PCR products; 39 samples resulted positive with the primer set ITS1F/ITS4 and 22 with NSI1/NL4, respectively (Table S2). Seven samples produced positive results with both sets of primers. Nested PCRs from the primer set NSI1–NL4 yielded strong bands but they were not directly sequenceable and, after cloning, the PCR produced only very few colonies. In total, 130 full-length non-chimeric ITS sequences clustered in 39 OTUs, 23 of which resulted in unique genotypes. In six cases, identical sequences for the same sample were obtained from both primer sets.

Belowground fungal communities were dominated by Ascomycota with 104 sequences (80%) and Basidiomycota with 25 sequences (19.23%), whereas the phylum Glomeromycota with one sequence (0.77%) was less frequent (Fig. 2). When we clustered these sequences into OTUs, 24 were classified as Ascomycota, 14 as Basidiomycota and one as Glomeromycota. Within the Ascomycota, the OTUs #20 (39 sequences), #7 (30 sequences) and #26 (five sequences) contained the most sequences and in the OTUs with the most sequences were the Basidiomycota #35 (four sequences), #5, #8 and #0 (three sequences each). The most frequent OTUs, #20 and #7, were detected in seven and eight different samples from conventionally and organically managed sites, respectively. The remaining OTUs were found either in a single sample or in up to four different samples. Up to 12 different fungal sequences (#20 on C2) from the same OTU were generated from the same sample. Within each farming system, some OTUs were found in samples along the entire transect. The OTUs show a nearly linear accumulation and the OTU sampling did not reach saturation (Table 1). Out of all OTUs detected, 10 occurred exclusively in conventional farms (three non-singletons) and 19 in organic farms (one non-singleton), whereas 10 OTUs were present in both farming systems. Based on a comparison with reference sequences, several OTUs have broadly disjunct distributions, either continental or locally restricted. The diversity indices of the belowground fungal communities are given in Fig. 3.

Belowground fungal community composition

Belowground fungal communities from both conventional and organic farming systems showed non-significant MPD values, whereas the MNTD values were significant for organic farming systems only (p < 0.05) (Table 1).

Environmental factors of different farming systems

The main physical and chemical properties of the soils are shown in Table S3. The C/N ratios were higher in soils under organic farming practices than those under conventional practices. On the contrary, the related parameters total N, total C, organic C, phosphorus and humus contents had some values that overlapped between both management systems. Nonmetric multidimensional scaling (NMDS) showed the clear separation of farming systems according to the soil parameters analysed (Fig. 4A).

Variable selection with tb-RDA identified three significant variables that account for 43.51% of the variation (Ntotal (p = 0.020), Stotal (p = 0.035), Pcal (p = 0.005)) of the OTU distribution (Fig. 4B). The PERMANOVA, however, did identify the farming system as only having a minor effect, while not being significant (p = 0.086, R² = 0.202).

Functional guilds

The fungal OTUs or species were grouped into 14 ecological guilds distributed across the following trophic modes: pathotrophs (two OTUs), saprotrophs (four OTUs), symbiotrophs (three OTUs), pathotroph-saprotrophs (six OTUs), pathotroph-symbiotrophs (four OTUs), pathotroph-saprotroph-symbiotrophs (four OTUs) and those with an unknown position in the food chain (16 OTUs) (Fig. 5; Table 1). In general, a higher diversity of trophic modes (involving some symbiotroph OTUs) were found on organically managed sites than in conventional farming systems (Table 1).

Taxonomic coverage and belowground fungal community diversity

The low number of PCR-positive samples obtained is possibly associated to the sample type rather than to the primers’ specificity and sensitivity. We combined two sets of primers since only seven samples were positive with both sets. The remaining samples were positive with only one primer set. We suspect that factors such as the inadequate disintegration of fine roots or the presence of inhibiting substances from the host plant could have affected the detection success of fungi (PCR amplification) (Schrader et al., 2012).

In this study, in roots of T. aestivum, up to eight different OTUs were identified per sample, which indicates that the roots of one plant were colonised by several fungal species simultaneously. Verbruggen et al. (2014) examined the fungal communities associated with wheat roots from organically and conventionally managed farming sites in Switzerland. They detected a remarkably greater number of OTUs (493) by sampling a smaller number of individual wheat plants from a larger number of sites. The higher diversity detected in their study may be partly caused by the different environmental conditions and sampling strategy, but was primarily caused by the molecular methodology used: they applied next-generation sequencing, in contrast to the Sanger sequencing we used here. Similar to this study, they identified most of the detected fungi as Ascomycota, with two of the most common OTUs clustered within the Sordariomycetes, and a few Basidiomycota (Fig. 2; Table 1). Similarly, Carter et al. (1999), who isolated fungi from the roots of field-grown wheat and used cultivation methods and subsequent fungal characterisation by rDNA sequencing, found the Ascomycota to be more frequent than the Basidiomycota.

Although the organic farming sites and two of the conventional farming sites were located within short geographical distance from each other, each site had an almost unique OTU composition and the overlap of OTUs was low (Fig. 3). Conversely, some OTUs were frequently detected across several sampling sites, independent of the farming method. The OTU #20, identified as Microdochium sp., was more commonly detected in conventional than in organic sites. Members of the genus Microdochium are well known as plant pathogens which can cause head blight of cereals (Xu et al., 2007). The fact that this fungus was isolated from roots of healthy wheat plants confirms the common definition of endophytes, which also includes latent pathogens that at some stage of their lifecycle inhabit plant tissues without causing apparent harm to their host (Petrini, 1991; Sun & Guo, 2012). The OTU #7 (without any assignment beyond Ascomycota in the ITS RefSeq database) was closely related to Lasiosphaeriaceae spp. (sequences from GenBank database) and it has also been isolated from roots of T. aestivum in Sweden (Grudzinska-Sterno et al., 2016). Lasiosphaeriaceae spp. apparently include endophytic species with a broad host range; some members of this fungal family have also been isolated as endophytes of various host plants from diverse habitats (Arnold et al., 2000; Su, Guo & Hyde, 2010). In addition, OTU #26, which was identified as Fusarium equiseti, was detected in all conventional farming sites. F. equiseti is a ubiquitous soilborne pathogen that causes vascular wilt on a wide spectrum of plants (Dean et al., 2012).

Impact of agricultural practices on belowground fungal communities

The present study revealed that the fungal communities are phylogenetically clustered and have significant NRI values, indicating that each farming method contains a distinctive fungal community (Table 2). Here, the farming method does influence the structure and also the fungal OTU diversity. In particular, these results highlight that fungal community assembly could result from similarities in response to some environmental factors, such as abiotic habitat filters, plant-imposed habitat filters or interactions with other fungal species (Saunders, Glenn & Kohn, 2010). Conventional farming systems, in contrast to organic systems, are characterised by the application of agrochemicals, which can alter the soil’s chemical composition and influence fungal diversity directly (Ellouze et al., 2014). Verbruggen et al. (2014) found that different agricultural management practices did not have a significant effect on the fungal communities associated with roots of T. aestivum, but they observed a significant effect of some abiotic soil factors like pH and magnesium concentration. The results for this analysis were however obtained by pooling sequences from different sites under the same management type, and further studies with an increased replicate number would be needed to verify the results obtained here. In this study, total N, total S and Pcal explain around half of the variation in abundance and composition of the fungal communities associated with winter wheat. Studying the long-term effects of nitrogen and phosphorus fertilisation on soil microbial community structure and function under continuous wheat production, Li et al. (2019) observed that N fertilisation increased fungal abundance. Similarly, Jie et al. (2012) found that sulphur affects the composition of arbuscular mycorrhizal fungal communities of different soybean cultivars. Cai et al. (2019) working with a rice-fish system found that total N and available P regulated the abundance of dominant fungi. When assessing the influence of the farming method directly, however, only a weak (and non-significant) effect could be observed, which could indicate that the underlying soil parameters play a stronger role than the type of management (REF).

One of the most striking differences between sampling sites from the two farming systems was the abundance and diversity of the co-occurring vegetation, which was generally higher in organic farming sites. As the OTU diversity was similarly high in organic sites, we postulate that the co-occurring vegetation could serve as a reservoir source of some fungal taxa. Verbruggen et al. (2014) examined the diversity of co-occurring weed communities and their ability to significantly influence fungal diversity and suggested that this might be partly a result of the weed litter left in the soil, which could stimulate the growth of endophytes with saprobic activity. Furthermore, they suggest the influence of the co-occurring plant community if some endophytic fungi with a broad host range proliferate more extensively in neighbouring plant species than in wheat itself. Besides that, there is evidence that plants release a wide array of compounds via their roots that can act as nutrient sources for soil fungi and can stimulate their proliferation (Ellouze et al., 2014). Therefore, a wide array of co-occurring plant species that release equally diverse compounds may provide a basis for the establishment of a wide array of fungi. Riess (2009) studied endophytes associated with cereals, including T. aestivum, and the co-occurring herbaceous plant species and observed that in many cases, the neighbouring plants were colonised by the same fungal phylotype. This author suggested that the associated fungi could connect their host plants in belowground networks in a similar manner to ECM (Onguene & Kuyper, 2002) or arbuscular mycorrhizal fungi (Voets et al., 2009). This could also be supported by the fact that a great number of endophytic fungal species are primarily transmitted horizontally (Rodriguez et al., 2009).

The predicted functional guilds of the fungal sequences also varied between conventionally and organically grown winter wheat (Fig. 5; Table 1). Interestingly, an increase in the occurrence of potential plant pathogen fungi and a decrease in the symbiotrophs were predicted for organic farming. The fact that fungal orders containing mainly endophytic species were detected in the roots, as well as mycorrhizal, parasitic or saprobic species, is not surprising. For instance, members of the Glomeromycota, which comprise arbuscular mycorrhizal fungi and were detected twice in this study, can typically grow within plant roots (Parniske, 2008). Additionally, it is difficult to distinguish fungal species on the basis of their lifestyles, because some species change their ecological strategy during their lifecycle. For example there is evidence that some endophytes can also adopt saprobic lifestyles at senescence of their host plant (Promputtha et al., 2007; Porras-Alfaro & Bayman, 2011).