The legacy effect of cover crops on soil fungal populations in a cereal rotation

The use of rotations and minimum tillage in agriculture can permit more sustainable production through increasing soil organic matter and nutrients, and breaking of pathogen lifecycles. Soil fungal populations make an important physical and chemical contribution to soil. For example, mycorrhizal species are important in plant nutrition but are often overlooked when considering management practices for ef ﬁ cient soil function. We undertook DNA metabarcoding (Ion Torrent) using novel PCR primers and high-throughput sequencing of the D1 region of the large ribosomal subunit of the rRNA locus, to assess the effect of different forages and cereal tillage methods on the soil fungal community. The study comprised ﬁ ve forage treatments, perennial ryegrass ( Lolium perenne ) with either low or high N, chicory ( Cichorium intybus ), red clover ( Trifolium pratense ) or white clover ( Trifolium repens ) grown over 3 harvest years (2010 – 2012). Cultivation of chicory, red clover or white clover led to signi ﬁ cantly divergent soil fungal communities, with a notably lower diversity of fungal populations under clover, suggesting a link to soil N dynamics. Consistent with this, was a negative correlation of soil nitrate-N levels with populations of arbuscular mycorrhizal fungi (AMF) and other root-associated fungal groupings (dark septate endophytes, ‘ CHEG ’ , Sebacinales and Ceratobasidiaceae). In contrast, abundance of Fungi belonging to the genera Mortierella and Cryptococcus were positively correlated with soil nitrate-N, with Mortierella also being negatively correlated with soil P. Spring wheat was sown on the same plots (April 2013) followed by winter barley (October 2013). Half of each plot was sown either after ploughing or by direct drilling. A legacy effect of the preceding forage crop on the fungal community was detected after both cereal crops, with plots previously cultivated with ryegrass being most divergent. No overall effect of establishment method on fungal communities was detected but AMF and CHEG fungi were more abundant on direct-drilled plots and pathogenic fungi were more abundant on ploughed plots after the sowing of winter barley. (http://creativecommons.org/licenses/by/4.0/).


Introduction
Given the need for the sustainable land use in agriculture, there is an urgent need to understand the complexities of plant-soilmicrobe interactions, the role of microbes in plant nutrition and maintenance of soil function and how different agricultural management practices may alter the soil ecosystem and environment. The use of forage and cover crops in rotations and a reduction in tillage frequency, can play an important role in mixed farms by increasing the organic matter content of soils, thereby contributing to carbon sequestration and soil productivity (Six et al., 1998;West and Post, 2002). Such practices also improve the hydrological properties of soil, for instance water infiltration, moisture retention and erosion rates (Augé et al., 2001).
The benefit of leguminous forage crops, such as red (Trifolium pratense) and white clover (Trifolium repens) in providing nitrogen for following crops has long been appreciated (Ebelhar et al., 1984), but these and other forage crops can provide additional benefits within agricultural rotations. Forage crops can be an important element in the control of certain pathogens through breaking the lifecycle. For example, Gaeumannomyces graminis var. tritici causes root rot in wheat and other susceptible crops such as barley and rye but can be controlled via rotation with non-susceptible forages such as clover and chicory (Cook, 2003). Chicory has also been shown to reduce the levels of helminth parasite infection of grazing sheep by interfering with the free-living stages of the life cycle of the parasites within the pasture environment and hence reducing their ability to infect grazing livestock (Marley et al., 2006).
Kingdom Fungi comprises a morphologically diverse group of organisms ranging from single celled yeast to macrofungi that can form networks in soil over many metres. In the context of arable agriculture, Fungi attract attention as major crop pathogens, for instance wheat rust (Aime, 2006) and ear blight but they also play a central role in nutrient cycling through the catabolism of dead organic matter and as mycorrhizal mutualists (Smith and Read, 2007). Historically the arbuscular mycorrhizal fungi (AMF) (Glomeromycota) have been viewed as the main mycorrhizal symbionts in arable and grassland soils (Smith and Read, 2007;Joergensen and Wichern, 2008). However, recent advances in plant-soil interactions have revealed that a wider range of Fungi than previously suspected may be involved in mutualistic interactions with higher plants in grasslands and thus play an important role in plant nutrition (Smith and Read, 2007;Heijden et al., 2015). Notable among these are the dark septate endophytes (DSE) a diverse group of filamentous ascomycetes, many belonging to the order Helotiales (Wilberforce et al., 2003;Mandyam and Jumpponen, 2014), as well as many members of order Sebacinales  and family Ceratobasidiaceae (Veldre et al., 2013). These groups comprise taxa with highly divergent mycorrhizal morphologies, which, depending on the host plant, can range from typical sheathing ectomycorrhizas to intracellular hyphal coils or undifferentiated intercellular hyphae. These groups also include a few taxa are known to cause plant disease and for the majority there is only sporadic evidence that the interactions are mutualistic (Jones and Smith, 2004). Thus the outcomes of these symbioses may vary with particular plant/fungus combinations and specific environmental conditions (the "mutualism-parasitism-continuum"; (Mandyam and Jumpponen, 2014)), with the interplay between direct and mycorrhizal pathways of nutrient uptake (Smith and Read, 2007). Therefore research to elucidate the role of these fungi in effective soil function is needed to assist in the development of our understanding of soil-plant interactions and, thereby reduce our reliance on inorganic inputs within agricultural systems.
The study of soil fungi, historically focused culture-based or microscopic approaches, has hampered attempts to reliably identify and quantify the abundance of different taxa. However, recent developments in DNA sequencing technology and the expansion of large databases of sequence and taxonomic data, such as NCBI and curated database of fungal sequences such as RDP (Cole et al., 2014) and UNITE (Abarenkov et al., 2010), mean that sequencing and identification from environmental sample such as soil is now possible. As a result we are now beginning to understand the distribution and the changes in fungal communities from diverse habitats (Geml et al., 2014;Jumpponen and Jones, 2014;Tedersoo et al., 2014;Franke-Whittle et al., 2015).
Here we present an analysis of fungal communities under a crop rotation of forage crops and cereal crops using a novel primer combination that amplifies sequences from a diverse range of fungi from basal groups such as Chytridiomycota to the more advanced dikarya (Basidiomycota and Ascomycota). The aims of this paper are to determine the effect of five different forage crop regimes on soil fungal communities and then to monitor the effects of either direct drill or ploughing on establishment of cereal crops on these soil communities. It was hypothesised that the type of forage crop will lead to the development of divergent soil fungal communities and that direct drilling to establish follow-on cereal crops will cause less disruption to these communities than ploughing.

Forages
The study comprised of 5 forage treatments, perennial ryegrass (Lolium perenne) with either low (PRG0N) or high N (PRG200N), chicory (Cichorium intybus)(CY), red clover (Trifolium pratense)(RC) or white clover (Trifolium repens)(WC) grown over 3 harvest years (2010)(2011)(2012). The experimental plots were established at the Institute of Biological, Environmental and Rural Sciences (IBERS), Aberystwyth University (52.4331, -4.0261) on an area of stony, well-drained loam of the Rheidol series. Full details of the experimental site, plot establishment and crop management are described in (Crotty et al., 2015) and a summary schematic of plot layout and crop/sampling timelines is presented in Suppdata 1. In brief, after initial application of herbicide (4 l ha À1 ) to the initial ryegrass ley in April 2009, the area was ploughed, dolomitic limestone (5 t ha À1 ) and inorganic fertiliser (60 kg P 2 O 5 and 60 kg K 2 O ha À1 ) was applied. The area was power-harrowed and rolled before the experiment was established as a randomised block design with four replicate plots (7.5 Â 12 m), on the 29th June 2009. Five treatments (20 plots in total) were established, consisting of perennial ryegrass with either low N (80 kg N ha À1 , applied once in March 2011) or high N (200 kg N ha À1 yr À1 ; as four monthly applications: April-July), chicory (200 kg N ha À1 yr À1 ; as for high N ryegrass), white clover and red clover (no N fertiliser addition). Plots were seeded at a rate of 33 (ryegrass), 16 (red clover) and 6 (chicory, white clover) kg seed ha À1 . Plots were cut and herbage removed five times each season (monthly: May-September; from 2010 to 2012). All plots received P and K fertiliser (as required to maintain a soil P index of 2), and N fertiliser was applied to the relevant plots as ammonium nitrate.

Cereals
In February 2013, all plots were treated with herbicide (360 g l À1 glyphosate at 4 l ha À1 ; Gallup 360) and each plot was split into sub-plots (3.75 Â 12 m) and allocated at random to two cultivation treatments, ploughing (P) or direct drill (DD) giving a total of 40 plots. The P sub-plots were ploughed to a depth of 175 mm (20th March) and power-harrowed (4th April), whilst DD sub plots were undisturbed prior to sowing. Spring wheat (Triticum aestivum) was sown (253 kg seed ha À1 ; 5th April) and fertiliser was applied with the seed (49 kg N, 9 kg P 2 O 5 , 28 kg K 2 O and 16 kg SO 3 ha À1 ). Prilled lime was top dressed at 370 kg ha À1 post-sowing to avoid any tillage effects on soil pH. A second fertiliser application (21st May) supplied 127 kg N, 22 kg P 2 O 5 , 72 kg K 2 O and 42 kg SO 3 ha À1 to achieve sufficiency status according to RB209 recommendations (soil N index 1 to target the recommended spring wheat N application of 180 kg N ha À1 ). Wheat was harvested on 29th August.
On 10 October 2013 all plots were treated with herbicide (4 l ha À1 ) and the same sub-plot establishment treatments remained for barley. The P sub-plots were ploughed on the 14th October 2013 and power-harrowed on 15th October. Winter barley (Hordeum vulgare) was then sown (196 kg ha À1 ) on the 15 October 2013. No fertiliser was applied during establishment. Fertiliser was applied on the 13 March 2014 (52 kg N, 52 kg P 2 O 5 and 52 kg K 2 O ha À1 ) and prilled lime was top dressed at 370 kg ha À1 . On the 16th April 2014 a second fertiliser application supplied 80 kg N, 14 kg P 2 O 5 , 45 kg K 2 O and 24 kg SO 3 ha À1 and the crop was harvested on 15th July.

Measurements
Soil samples were taken from each plot with a 25 mm auger to a depth of 15 cm (fourteen cores per plot) on 14th September 2012 (prior to removal of forage crops 20 plots), 16th September 2013 (2 weeks after harvesting of spring wheat 40 plots) and 13th August 2014 (4 weeks after harvesting of winter barley 40 plots). The cores were bulked into a single sample per plot (fresh wt approx. 400 g in total) giving a total of 100 samples.

Soil chemical analyses
The soil chemical analyses were performed as follows: Bulked samples from each plot were sieved through a 12 mm sieve, homogenised by hand, sub sampled for multiple analytical procedures. They were stored at 4 C in the dark prior to analysis.
Soil mineral-N was determined as nitrate (NO 3 À ÀN) and ammonium-N (NH 4 À ÀN) after extracting fresh soil samples with 2 M potassium chloride solution using an extractant/soil ratio of 5:1. Concentrations of NO 3 (including any nitrite present) in the extracts were determined colorimetrically using the method of Henriksen and Selmer-Olsen (1970). NH 4 À ÀN concentrations in extracts were determined using the method of Krom (1980).
The following procedures described below were carried out on air-dried and milled samples. Olsen P was determined by shaking 5 g of soil with 100 ml sodium bicarbonate (pH8.5) for 30 min. The extracts were filtered and P determined colorimetrically at 880 nm. Mg, K, Ca and Na were extracted by shaking 10 g of soil with 50 ml of 1 M ammonium nitrate for 30 min and then filtered. Mg, Ca and Na were measured by inductively coupled plasma atomic emission spectroscopy (ICP-AES) (Varian Liberty AX, Agilent, UK) and K using a flame photometer (Corning 410, Sherwood Scientific Ltd., Cambridge, UK). S, B, Fe, Mn, Cu, Zn and Al were measured using Mehlich III extractions at Brookside Laboratories, Knoxville, Ohio. Soil pH was determined after mixing 10 g of soil with 50 ml deionised water. The solution was allowed to settle for 30 min before the pH was measured using a SCHOTT CG 843 P pH meter (SI Analytics GmbH, Germany).

DNA amplification and sequencing
The soil samples from the three annual timepoints detailed in section 2.1 were frozen at À80 C and freeze-dried. Each sample was ground through a 500 mm sieve with care taken to thoroughly mix the soil samples (ca. 300 g dw each) before and after grinding. DNA was extracted from 200 mg of the sieved soil sample using a PowerSoil DNA extraction kit (MO BIO Laboratories, Inc. Carlsbad USA).
The D1 region of the large sub unit (LSU) of ribosomal DNA (rDNA) was amplified using fungal specific novel primers; D1F2 (CYYAGTARCTGCGAGTGAAG) and the reverse NLC2AF (GAGCTG-CATTCCCAAACAA). The reverse primer is similar to the NLC2 primer (Martin and Rygiewicz, 2005) but missing two bases at the 3 0 end. Amplification (one per DNA sample) was performed in a 50 ml Polymerase Chain Reaction (PCR) using Promega GoTaq G2 DNA polymerase (Promega, Madison USA). Each reaction contained 250 nM of each primer, 2.5 mM MgCl 2 , 1 mg ml À1 BSA, 200 mM dNTPs and 0.5 U of DNA polymerase in the supplied buffer.
Primer D1F2 was linked at the 5 0 end to the IonTorrent A-adapter sequence (CCATCTCATCCCTGCGTGTCTCCGAC), the TCAG key and an IonXpress Barcode. Primer NLC2AF was linked at the 5 0 end to Ion Torrent B adapter sequence (CCTCTCTATGGGCAGTCGGTGAT).
The PCR conditions were 94 C for 5 min (initial denaturation) followed by 30 cycles at 94 C, 30 s (denaturation); 52 C, 30 s (annealing); 72 C, 30 s (extension) and a final extension at 72 C for 5 min. The PCR reactions were cleaned using spin columns (NBS Biological, Huntingdon UK) and the amplified DNA quantified using NanoDrop (NanoDrop Products, Wilmington USA). The 100 samples were pooled in equimolar concentrations into 2 aliquots of 50 samples each and were further purified using E-gel, then analysed and quantified using a Bioanalyser 2100 (Agilent Technologies, Santa Clara, USA). The pooled samples were adjusted to a final concentration of 12 pM. Emulsion PCR was carried out using Ion-Torrent 200 bp template kit on the Ion Torrent One Touch 2 system according to the manufacturer's instructions. Amplified sequence particles were then enriched using the One Touch ES to remove non-template particles and sequenced on "316" (100 Mbp) microchips using the Ion Torrent Personal Genome Machine (Life Technologies, Waltham USA).

Sequence data processing
Sequence data were quality checked, trimmed to 200 bp and demuliplexed using MOTHUR (v. 1.31.2; (Schloss et al., 2009)). Sequences with mismatching barcode and primer sequences (length less than 100 bp and with an average Phred score less than 20) were discarded. No evidence was seen of any decline in Phred scores over the length of the sequence after quality checking (Suppdata 2). Sequences were then checked for putative chimeric sequences using the UCHIME algorithm (Edgar et al., 2011) against a reference database of curated fungal LSU sequences downloaded from the Ribosomal Database Project (RDP) website (Cole et al., 2014). Sequences for each sample were rarefied to the lowest sequencing depth using the subsample function of MOTHUR to randomly select sequences from each file. Rarefied files were dereplicated and singletons discarded as recommended in Tedersoo et al. (2015), and operational taxonomic units (OTUs) assigned using USEARCH/UPARSE (v7 (Edgar, 2013)) at 97% clustering ; clusters containing fewer than 10 sequences were discarded. The delineation of OTUs at 97% was chosen empirically to determine the level at which genera assignments determined by the RDP classifier were affected. Clustering at higher levels of similarity (98-99%) did not change genera assignments or relative abundances whereas clustering at lower levels (95-96%) did reduce the numbers of genera identified.
A taxonomy was assigned to each OTU (operational taxonomic unit) using the Naïve Bayesian Classifier (Wang et al., 2007), with a recommended cut off for 200 bp sequences of 50%, against the RDP LSU database version 11. This version of the database contains ca 12000 fungal sequences in 3200 genera and in contrast to previous versions contains ca. 2000 non-fungal sequences. Additionally, we have added ca. 200 sequences to our in-house database, in particular to provide reference sequences for phylum Cercozoa (kingdom Rhizaria) which comprised the majority of non-fungal sequences discovered in this study. This ensured that all nonfungal sequences were identified at least phylum level. An example of the output spreadsheet is shown in Suppdata 3.
Where genus was not assigned by the classifier (but only assigned to family or order, due to confidence being lower that the threshold), an OTU identifier was assigned to that cluster. Data were then rendered in Excel and standardized by dividing the number of reads in each taxonomic unit by the total number of fungal reads in each sample to give relative abundances of the assigned taxa for each quadrat; any non-fungal taxa were reported separately. Shannon (À X S i¼1 P i lnP i where P i = relative proportion of the ith taxa) and Simpson (1= X S i¼1 P 2 i Þdiversity indices were calculated for each quadrat. Broad ecological function of the fungi was assigned to each taxon (where identified) at genus or family level through searches of academic literature (Rinaldi et al., 2008;Öpik et al., 2010;Tedersoo et al., 2010;Mandyam and Jumpponen, 2014), following an approach similar to Tedersoo et al. (2014). If different ecological functions were identified within a taxon, a function was only assigned when >75% of known species within that taxon could be assigned to a single function. Otherwise, the function remained undetermined. Five groupings, mostly associated with plant roots were identified: SAP (saprotrophic fungi; fungi not colonising plant roots which decompose organic matter), PATH (pathogens; defined having a known negative effect on plant growth), AMF (arbuscular mycorrhizas; members of phylum Glomeromycota), DSE (members of the ascomycete families Helotiaceae, Dermateaceae, Herpotrichiellaceae, Hyaloscyphaceae: ascomycete endophytes with dark melanised septa) and CHEG fungi (members of the families Clavariaceae, Hygrophoraceae, Entolomataceae and Geoglossaceae). This last group comprises taxa mainly associated with grassland habitats (Griffith et al., 2002). Ecological, microscopic, isotopic and phylogenetic evidence is accumulating (Griffith et al., 2002;Tedersoo et al., 2010;Birkebak et al., 2013;Halbwachs et al., 2013) to suggest that these fungi are biotrophic, rather than saprotrophic as previously thought (Rinaldi et al., 2008) but the precise nature of their interaction with plant hosts is unknown. Finally, the assigned ecological functions were then validated using the FUNGUILD database (Nguyen et al., 2015) for consistency. The RDP phylogenetic structure was input to a modified FUNGUILD Python script (https://github.com/UMNFuN/ FUNGuild) to return the classification for the RDP database from FUNGUILD (http://www.stbates.org/funguild_db.php) and compared against the classification used here.
Sequence data has been submitted to the European Nucleotide Archive with reference number PRJEB12560.

Data analysis
The R package was used for visualisations of relative abundance matrices using non-metric multidimensional scaling (NMDS) ordination to identify patterns in the data. Analyses of variance (ANOVA) tests were performed on the ordination axes to determine if data varied on the primary or secondary axis. Permutation multivariate analysis of variance (PERMANOVA) was used to determine overall significant differences in community data and was performed in PRIMER 6 & PERMANOVA+ (versions 6.1.12 and 1.0.2 respectively; Primer-E, Ivybridge, UK). Abundance percentage data were subjected to square root transformation and Bray-Curtis distance matrices calculated. PERMANOVA was carried out using default settings with 9999 unrestricted permutations and the Monte Carlo P value was calculated. PERMANOVA was executed with a nested design for the 2013 and 2014 sample data to account for the split plot.
Analysis of Similarity (ANOSIM) was carried out in PRIMER 6 & PERMANOVA+ using the Bray-Curtis distance matrix calculated above. ANOSIM was executed with a nested design for the 2013 and 2014 sample data to account for the split plot. This analysis was used to provide a metric of the degree of divergence between communities as given by the R statistic.
To calculate the contribution of environmental data on fungal communities, distance based linear modelling (distLM) was used to calculate which environmental variables had a significant correlation with the community data. Significant variables were used in distance based redundancy analysis (dbRDA) (Legendre and Anderson, 1999) as implemented in PRIMER 6 & PERMANOVA +, with predictor variables normalised before analysis. All other statistical analyses were performed in Excel and GenStat 1 . The effect of treatments on relative abundances after the establishment of the wheat crop were analysed according to the split-plot design, with effects of previous forage treatment estimated in the whole plot stratum and effects of establishment method and previous forage-establishment interactions estimated at the subplot level.

Effect of different forage crops on soil fungal populations
Following cultivation of the forage crops after 3 harvest years (September 2012 sampling), a significant treatment effect on fungal diversity (P = 0.042) was found, with a lower Simpson diversity index lower (post hoc Tukey's test) in soils under RC or WC (mean 30.1 and 31.6 respectively) than under PRG200N (mean 39.3). The NMDS ordination (Fig. 2a) illustrates the separation of the fungal communities data under different forage crops. The PERMANOVA main test also showed significant treatment effects (P = 0.001; Pseudo F = 1.413), with pair-wise tests showing that the clover plots were divergent (Table 2) and this was confirmed by ANOSIM (P < 0.05), with the greatest difference between clover forages and ryegrass. Two negative but non-significant ANOSIM R values close to zero were obtained for the PRG0N treatment, which may indicate a larger degree of replicate variability for this as opposed to other treatments (Chapman and Underwood, 1999).
The relative proportions of the particular functional/ecological grouping (e.g. arbuscular mycorrhizal [AM] etc., following classifications consistent with Nguyen et al. (2015)) under the different forage crops were also tested. Significant differences were found for CHEG, DSE and AMF (Table 3), with CHEG most abundant under the low N rygrass plots, AMF under both ryegrass treatments and DSE most abundant under red clover. Analysis of the changes in the abundance of specific fungal taxa that were associated with the  PRG200N and WC being the most divergent forage treatments (Fig. 2a). DistLM analysis showed significant correlation with three environmental variables, nitrate-N (P = 0.001), calcium (P = 0.022) and manganese (P = 0.045), which together explain 21% of the variation in the fungal community data. Similarly, dbRDA analysis (Fig. 2b) shows the correlation of nitrate-N with the fungal community. Primary axes accounted for 16% of the variation and white clover plots are associated with high scores on the first axis. These findings are consistent with the significant effect (ANOVA P < 0.001) of forage crop treatment on soil Nitrate-N (Table 3).
In order to identify which particular fungi and which functional groupings of Fungi were most affected by these differences in nitrate-N, the relative proportions of the top 30 taxa were correlated against the values for nitrate-N for each plot. The abundance of eight taxa was found to correlate significantly with nitrate-N levels (Table 4). Of these, only Mortierella and Cryptococcus, the two most abundant genera across all samples (see above) showed a positive correlation. Of the six taxa showing a negative correlation, two are known to be mycorrhizal, namely Rhizophagus (AMF) and Sebacina; two others are root-associated and putatively mycorrhizal (Clitopilus [Entolomataceae] and Microscypha [Helotiaceae] (Tedersoo et al., 2010)).
Analysis of the correlation of the abundance of functional groups with environmental variables, revealed a significant negative correlation for CHEG (Table 4). In contrast, the abundance of pathogenic fungi (of which Didymella and Pyrenochaeta were the most abundant) was positively correlated (P = 0.009; r = 0.569) with nitrate-N.
As noted above correlations of community data with environmental variables were strongest for nitrate-N. However, Mortierella, the most abundant taxon, was also significantly negatively correlated with soil P concentrations (P = 0.02 r = À0.514)

Effect of the first cereal (spring wheat) crop
At the second sampling time point in September 2013, the spring wheat crop sown (either by direct drilling or ploughing) had  been harvested from plots previously sown to different forages. Multivariate analysis of fungal community data (Table 2) showed a significant legacy effect for preceding forage type (PERMANOVA main test: P = 0.001; Pseudo F = 1.653) but not establishment method (P = 0.744; Pseudo F = 0.8283). Similarly, ANOSIM tests (Table 2) showed that plots previously growing ryegrass differed from those with chicory, white clover or red clover. The degree of separation was smaller than that found in 2012, consistent with a poorer separation of plots in NMDS ordination (Fig. 3a,b). Fungal diversity (as measured by Simpson/Shannon indices) was also unaffected by either the preceding forage crop or wheat establishment method.
DistLM analysis showed nine environmental variables had a significant effect on the fungal community (sulphur, phosphorus, potassium, boron, iron, copper, zinc, macropores% volume and fine earth bulk density). Together these explained 35% of the variation in the fungal community data. The dbRDA primary axes displayed 19% of the variation but no separation was apparent by forage crop (Fig. 4a) or by establishment method (Fig. 4b).
The significant environmental variables were correlated against the top 30 taxa and in total 53 significant correlations were found. Of these, only six had a correlation coefficient (r) > 0.5 or < À0.5 (i.e. explaining more than 25% of the variation in the data; Table 5). Only one functional group, (AMF), showed significant correlation with an environmental variable, being negatively correlated with potassium (r = À0.526, p <0.001). The effect of previous forage crop or establishment method on functional groups was tested via split plot ANOVA. Only one functional grouping (ascomycete yeasts) showed a significant response to forage crop but there was a significant interaction with establishment method.

Effect of the second cereal (winter barley) crop
The final soil sampling was undertaken in August 2014, 10 months after the sowing of winter barley and 4 weeks after the crop was harvested (Suppdata 6). Multivariate analysis of fungal communities showed a continued legacy effect of the original forage crop (PERMANOVA P = 0.024; Pseudo F = 1.3692) but establishment method was not significantly different (P = 0.085; Pseudo F = 1.443) and ANOSIM analysis (Table 2) show that fungal populations differed significantly between CY and PRG0N and PRG200N, and also between WC and PRG0N plots (the largest separation being between CY and PRG0N [r statistic = 0.198]). In addition, ANOSIM shows a small (r=0.11) but significant (P = 0.01) separation between establishment methods.
Greater fungal diversity was present in direct drilled than ploughed plots (Simpson index 22.79 and 18.88 [F = 11.879 P = 0.004]; Shannon index 3.94 and 3.77 [F = 16.051 P = 0.001] respectively). This difference was also reflected in NMDS ordination (Fig. 5). A split-plot ANOVA test on the scores for the first and second axes showed a significant treatment effect (Axis 1 P = 0.039; F = 5.14. Axis 2 P = 0.010; F = 8.65).
Several environmental variables showed significant effects on the fungal community (elements Na, S, P, K, Fe, Zn, Ca and fine earth bulk density), together explaining 35% of the variation in the data. The dbRDA primary axes displayed 19% of the variation but no separation was apparent by forage crop (Fig. 6a) or by establishment method (Fig. 6b). The significant environmental variables were correlated against the top 30 taxa and in total 46 significant correlations were found. Of these only three had a correlation coefficient (r) > 0.5 or < À0.5 (i.e. explaining more than 25% of the variation in the data) (Table 5). No functional groups showed a positive correlation with any environmental variable.
The effect of previous forage crop or establishment method on functional groups was tested via split plot ANOVA. There was a significant difference of establishment method on both CHEG (P < 0.001) and AMF (P < 0.001) but no interaction with forage crop. All forages had a lower relative abundance of AMF and CHEG on the ploughed plots (1.05% and 2.10% respectively) than the direct-drilled plots (1.49% and 2.8% respectively). In contrast, those fungi defined as pathogenic had a significantly higher abundance on ploughed plots (6.92% vs. 5.46%; F = 7.329 P = 0.016). There was no significant effect on DSE by either forage crop or establishment method.

Comparison between years
Comparison of fungal communities between all the samples taken during the course of this study showed a significant effect of sampling time (Fig. 7; PERMANOVA P = 0.001; Pseudo-F = 11.524). ANOSIM shows all years are significantly different (P = 0.001), with the largest difference between 2012 and 2014 data (R statistic 0.564) and the smallest between 2013 and 2014 (R statistic 0.248).

Discussion
This study set out to determine whether the cultivation of different forage crops within agricultural systems affected soil fungal populations and how the degree of soil disturbance associated with the subsequent establishment of cereal crops would affect these communities. Across the whole dataset, the fungal communities remained broadly stable in response to the treatments imposed ( Fig. 1; Suppdata 4,5), in contrast to the more substantial differences in community composition reported when arable soils managed by organic and conventional methods are compared (Hartmann et al., 2014;Lentendu et al., 2014). Our work was done from large field-scale plots and care was taken to sample at multiple points (14) for each 90 m 2 plot, following by extensive mixing to ensure sample homogeneity prior to subsampling. As noted by Penton et al. (2014), this is important to avoid microscale community variability.

Legacy effect of preceding forage crops and effect of soil tillage
Three years after the establishment of various forage crops, the fungal communities had diverged under the different treatments. Given the wealth of published evidence of the effect of the addition of synthetic N fertilisers on soil microbial communities (Lentendu et al., 2014;Paungfoo-Lonhienne et al., 2015;Zhou et al., 2015), it was expected that the high and PRG0N regimes would be clearly differentiated, due to the differences in synthetic N additions. However, this was not the case and may have been due to the fact that soil nitrate-N levels were actually highest under clover forages. ANOSIM analysis showed that the fungal community structure was more highly separated under the clover plots from the other forage treatments. These findings suggest that the release of bacterially-fixed N from legume roots might be an important determinant of soil fungal communities, although a more complex interaction involving root metabolites, such as flavonoids, may also be involved (Cesco et al., 2012). Fungal-feeding nematodes were also found to be significantly more abundant under WC/RC and CY at this experimental site (Crotty et al., 2015) and it may be the case that total fungal biomass (not quantified here) is greater under clover than grass, as found previously (Six et al., 2006). It is widely found that pH is an important determinant of community structure for soil Bacteria and to a lesser extent Fungi (Rousk et al., 2010;Tedersoo et al., 2014) but here pH played only a minor role, probably because pH was maintained (range 5.8-6.4) by regular lime applications, with no significant effects reported for any treatments on soil pH.
Sampling in 2013 and 2014 aimed to elucidate whether any legacy effects of previous forage crops on subsequent fungal populations could be detected and how these were affected by ploughing versus direct drilling when establishing subsequent cereal crops. Significant legacy effects were detected in both years, most pronounced between plots previously planted with PRG200N regime and those sown with WC/RC/CY ( Table 2). The mechanism whereby this legacy effect was mediated is not known here but it is unlikely to be solely Table 4 Those of the top 30 taxa from the 2012 sample that are significantly correlated with the level of nitrate-N. SAP = Saprotroph, AMF = Arbuscular mycorrhizal fungi, CHEG = Clavariaceae, Hygrophoraceae, Entolomataceae and Geoglossaceae, DSE = Dark septate endophyte, PATH = Pathogen.

Taxon
Order ( Group n/a n/a 2.9% (2.2-3.4%) AMF À0.315 0.177 Group n/a n/a 7.1% (6.4-7.4%) DSE 0.15 0.528 Group n/a n/a 0.2% (0.0-1.6%) CHEG À0.495 0.026 Group n/a n/a 6.0% (4.9-8.4%) PATH 0.528 0.018 due to synthetic N addition since the chicory cover crop also received 200 kg N ha À1 yr À1 in 2009-12, and all plots received synthetic N in 2013-14. Legacy effects in soil fungal populations could occur for a number of reasons, most obviously the particular qualities of different plant residues and consequently the fungal populations best adapted to decompose these (Allison et al., 2013). Additionally, differential colonisation of roots by mycorrhizal or other endophytic fungi, can in turn alter palatability of litter to soil invertebrates (Kostenko et al., 2012) and the residues of fungi themselves (e.g. chitin, melanins) have long residence times in soil and may have a  significant affect of the composition of subsequent microbial populations (Gleixner et al., 2002). Soil tillage has the effect of increase rates of decomposition and thereby release of nutrients to newly sown crops. Such radical disruption can affect soil fungal populations in several ways: i) The degree of disturbance, lower in direct drilling, tends to favour fungal growth and more K-selected species, as it allows for extension of larger mycelial networks (Wardle,1995); ii) Soil moisture under nontillage tends to be higher resulting in increased fungal biomass (Frey et al.,1999); iii) reduced incorporation of dead organic matter in nontillage systems and more surface decomposition of plant remains (Six et al., 2006). In previous comparisons of these two establishment methods, it has been reported that reduced tillage leads to increased fungal biomass and fungal:bacterial ratios (Six et al., 2006) and changes in bacterial populations (Carbonetto et al., 2014). Rates of root colonisation by AMF (Jansa et al., 2003) and diversity of AMF, as found here, are reduced by tillage (Säle et al., 2015). It might be expected that fungal populations would be greatly altered by tillage, especially for root-associated fungi and taxa forming extensive mycelial networks (eg. Agaricales). We found that CHEG fungi (which are long-lived and form large mycelial networks) were negatively affected by tillage, however other abundant taxa in Agaricales (e.g. Conocybe and Vascellum) were not significantly affected.

Metabarcoding of the LSU region
The metabarcoding method deployed here is novel in that the primers used are new and unusually amongst studies of fungal biodiversity (Lindahl et al., 2013;Tedersoo et al., 2014), we have examined the D1 variable region of the LSU gene rather than the adjacent ITS1/2 spacer regions, although comparison of these different approaches suggests that congruent data are obtained (Brown et al., 2014). For many but not all fungal groups (Schoch et al., 2014), more reference sequences are available via GenBank and UNITE databases for the ITS region Schoch et al., 2014) than for other LSU and other loci. However, for several fungal groups, e.g. phyla Glomeromycota (Öpik et al., 2010) and Neocallimastigomycota (Callaghan et al., 2015), orders Pucciniales (Aime, 2006) and Hypocreales, and families such as Clavariaceae (Birkebak et al., 2013), other loci are preferred and thus are better populated with reference sequences. Since the level of sequence conservation within the LSU region is greater and occurrence of insertion/deletions less frequent, the 'binning' of sequences to genera and higher taxonomic levels is more efficient. As a result of this, 79 and 70% respectively of sequences were attributed to family or genus level (with <4% of sequences not assigned to phylum), in contrast to the many datasets focused on the ITS region, where in some cases >30% of all sequences remain unidentified even to phylum level and thus any interpretation of data dependent on the remaining minority of sequences (e.g. Paungfoo-Lonhienne et al., 2015). The consequent disadvantage of use of the LSU locus is the reduced resolution at species level, although for many fungal groupings (e.g., Clavariaceae; (Birkebak et al., 2013)) reliable discrimination at species level is possible.
There is also much debate about the relative merits of the two ITS spacer regions and the particular primers used to amplify these (Tedersoo et al., 2014), and potential biases of the commonly used primers against Glomeromycota and Archaeorhizomycetes and Tullasnellaceae have been highlighted (Schadt and Rosling, 2015). To test the possibility that our methodology might bias against Glomeromycota, we undertook an analysis of LSU sequences and found the ten base positions nearer the 3 0 end of both primers to be well-matched (identical in >96% of sequences from GenBank). For Archaeorhizomycetes there were slightly more mismatches (Suppdata 7) but genera affected by these mismatches appear in this analysis. The fact that generic fungal-specific primers also find Glomeromycota to comprise 1-6% of total soil fungal DNA in habitats as diverse as sand dunes 5.29% (Geml et al., 2014), tallgrass prairie 3% (Jumpponen and Jones, 2014), pea fields <1% (Xu et al., 2012b), low input meadow 3.6% (Schmidt et al., 2013), pasture, grassland and woodland 1.5% (Lauber et al., 2008), grasslands and shrublands 1.4% (Tedersoo et al., 2014), suggests that our primers are not significantly biased against Glomeromycota.

Specific changes in the fungal communities
The composition of the fungal communities at phylum level was similar to that reported from other studies of soil fungi, dominated by Ascomycota (mean 42%; 26-57% of sequences), followed by Basidiomycota (mean 34%; 22-67%) and Zygomycota (mean 13%; 5-34%) (Hartmann et al., 2014;Lentendu et al., 2014;Liu et al., 2015). Although more than 800 fungal OTUs were identified across the 100 soil samples analysed here, 58-82% (mean 71%) of sequences are attributed to 30 of these taxa. In particular, three taxa (Mortierella,Cryptococcus and Veronaea spp.), one from each of the dominant phyla and each comprising on average 10-12% of total population, accounted for 27-46% of sequences across all samples. Members of the genera Mortierella and Cryptococcus have been found to be abundantly represented in other NGS investigations of soil fungi, for example featuring amongst the top ten taxa in arable (Xu et al., 2012b;Hartmann et al., 2014;Nallanchakravarthula et al., 2014;Liu et al., 2015) and woodland soils (Franke-Whittle et al., 2015). Furthermore, both taxa are readily isolated in soil plating experiments (Yurkov et al., 2015). In the current study, the abundances of the two taxa were positively correlated and Mortierella showed a positive correlation with soil N levels (Table 4) and negative correlation with soil P  (Table 5). Thus Mortierella may play a role in P mobilisation from inorganic sources via secretion of organic acids (Zhang et al., 2014). Mortierella has been found to increase plant P uptake in pot trials and to act synergistically with AMF (Osorio and Habte, 2013). The possibility that either Mortierella or Cryptococcus are endophytic/mycorrhizal is unlikely given that in a recent study of strawberry (Nallanchakravarthula et al., 2014), these taxa were amongst the most abundant in rhizosphere soil but absent from root tissues.
In attempting to interpret the complex lists of taxa and their relative abundance, ecological interpretation can be derived by assigning taxa to particular functional groups or guilds (Nguyen et al., 2015). Since such guilds often comprise phylogenetically related groups wherein particular traits are evolutionarily conserved, it is then possible to infer ecological function of poorly studied taxa based on their relatedness to well-studied species (Talbot et al., 2015). Although the cryptic nature of fungal growth makes this more difficult than for higher plants and animals, this approach has proven successful for some groups of fungi, for  example members of phylum Glomeromycota and the various groups of ectomycorrhizal fungi which dominate temperate and boreal woodlands (Smith and Read, 2007;Heijden et al., 2015). Several recent metabarcoding studies have defined functional or morphological grouping in the interpretation of their data (Clemmensen et al., 2013;Tedersoo et al., 2014) and there are attempts to broaden and standardise this approach (Nguyen et al., 2015). Here we have defined several functional/ecological groupings, one well-established (AMF), one morphological grouping (DSE-defined here as members of the family Herpotrichiellaceae [order Chaetothyriales] and several families within order Helotiales [Dermateaceae, Helotiaceae, Hyaloscyphaceae]) and one group of macrofungi (CHEG). This last grouping comprises both basidiomycetes (Clavariaceae, Hygrophoraceae, Entolomataceae) and ascomycetes (family Geoglossaceae) which co-occur in undisturbed nutrient-poor grasslands, especially in Europe (Griffith et al., 2002), and for which there is evidence of mycorrhizal trophic mode (Halbwachs et al., 2013).
For all three groups (AMF, DSE and CHEG), a negative correlation with soil nitrate-N was found, consistent with the pattern widely observed for mycorrhizal fungi from diverse ecosystems (Treseder, 2008), though eutrophication also leads to shifts in the relative abundances within these groupings depending on differential tolerance of elevated soil N (Lilleskov et al., 2002). Of the putatively mycorrhizal fungi found in the present study, the most abundant were DSE, of which Veronaea, Cadophora, Tetracladium and Microscypha spp. were the most abundant, together accounting for between 16% (2012) and 18.6% (2013) of all fungal sequences. Abundance of Sebacina, known to form diverse mycorrhizal associations  showed a similar negative correlation with soil nitrate-N and was more abundant in soils under ryegrass in 2012 (P < 0.001). More detailed analysis of these sequences (data not shown) found that all were identical to sequences from the roots of haymeadow plants . Unexpectedly (due to the absence of any woody hosts), sequences of several ectomycorrhizal basidiomycetes were detected, the most abundant being an unidentified OTU in family Thelephoraceae (mean 0.73%; present in 98/100 samples) which was identical to a cloned sequence from arable soil in Michigan (Lynch and Thorn, 2006). This suggests that some Thelephoraceae may form mycorrhizal or endophytic associations with non-woody hosts.
One key benefit of crop rotation is in the reduction of fungal root disease. However, reduced tillage, whilst having benefits for soil structure and erosion and soil tillage may exacerbate disease problems (Miller and Lodge, 2007). No obvious disease symptoms were recorded during the course of the current experiment and the yields reported were comparable to expected cereal yields for the region (data not shown). The taxa assigned to the ecological grouping 'pathogens' ('PATH', of which the genera Didymella, Pyrenochaeta, Cochiobolus and Thanatephorus were most abundantmean of 0.38-2.74%) showed a positive correlation with soil nitrate-N under the forage crops in 2012 (the opposite of the response shown by putatively mycorrhizal taxa), however no treatment effects were found.
The inverse response of pathogens and the putative mycorrhizal groups is consistent with the role of soil microbes in disease suppression in arable systems, for instance in the case of Fusarium vascular wilts and take-all decline (Garbeva et al., 2004;Penton et al., 2014). High levels of AMF colonisation of roots are associated with reduced levels of disease (Gianinazzi et al., 2010) and more recent NGS studies of fungi in arable soils have also shown inverse correlations in the abundance of soilborne pathogens and particular soil fungi. For instance Xu et al. (2012b) found abundance of Veronaea (DSE taxon, 3rd most abundant in the present study) to be inversely correlated with disease symptoms in pea crops.
Furthermore, isolate MSP6 identical in sequence to Cadophora (DSE taxon, 4th most abundant here) and also the most commonly isolated fungus from the roots of grassland plants has been shown to inhibit Fusarium colonisation of grass roots (Wilberforce et al., 2003). Cadophora is close to (and probably synonymous with) Leptodontidium spp., which is linked to suppression of Verticillium dahliae infection of strawberry (Nallanchakravarthula et al., 2014). Our data related to bulk soil and it is likely analysis of fungi inside roots would yield stronger correlations since it is at the surface and within internal tissues of roots that antagonism between pathogens and mycorrhizal fungi would be focused (Xu et al., 2012a;Penton et al., 2014).

Conclusions
This study has shown that forage crops significantly affect the soil fungal community with lower levels of root endophytic fungi on the white clover treatments and soil nitrate-N being a key determinant of the abundance of several ecological groupings. The contrasting legacy effects of these forage crops on soil fungal populations persisted under subsequent cereal crops. Reduced tillage had a lesser effect but did lead to increases in populations of AMF and CHEG fungi. Dias et al. (2015) have noted the importance of crop rotation for sustainable agriculture and the need to better elucidate the role of the soil biota in soil nutrient cycles and plant nutrition. The interactions between Fungi, Bacteria and other soil biota merit particular investigation, a challenge now becoming more tractable via metagenomic analyses (Carbonetto et al., 2014). The findings of this study highlight the potential importance of crop rotation for sustainable agriculture and the need to better elucidate the role of the soil biota in soil nutrient cycles and plant nutrition. Equally important is the negative role of some pathogenic fungi and the need to better understand the biotic component of disease suppression and the consequences of reduced tillage. Apparent from our study and many others relating to arable soils is the abundance of groups of soil fungi (Mortierella, Cryptococcus, Verronaea, Cadophora) about which very little is known. Since Mortierella is implicated in the mobilisation of phosphorus and Cadophora in the suppression of root infecting fungi, there is a need for more detailed pot-based investigation of the role of these common but neglected fungi in the efficient function of agricultural soils.