S. Emilia Hannula
Johannes A. van Veen- 1Department of Microbial Ecology, Netherlands Institute of Ecology (NIOO-KNAW), Wageningen, Netherlands
- 2Insititute of Biology, Leiden University, Leiden, Netherlands
Phylogenetic diversity of soil microbes is a hot topic at the moment. However, the molecular tools for the assessment of functional diversity in the fungal community are less developed than tools based on genes encoding the ribosomal operon. Here 20 sets of primers targeting genes involved mainly in carbon cycling were designed and/or validated and the functioning of soil fungal communities along a chronosequence of land abandonment from agriculture was evaluated using them. We hypothesized that changes in fungal community structure during secondary succession would lead to difference in the types of genes present in soils and that these changes would be directional. We expected an increase in genes involved in degradation of recalcitrant organic matter in time since agriculture. Out of the investigated genes, the richness of the genes related to carbon cycling was significantly higher in fields abandoned for longer time. The composition of six of the genes analyzed revealed significant differences between fields abandoned for shorter and longer time. However, all genes revealed significant variance over the fields studied, and this could be related to other parameters than the time since agriculture such as pH, organic matter, and the amount of available nitrogen. Contrary to our initial hypothesis, the genes significantly different between fields were not related to the decomposition of more recalcitrant matter but rather involved in degradation of cellulose and hemicellulose.
Introduction
Fungi are ubiquitous and diverse; the global number of fungal species has been estimated to be more than 1.5 million species (Hawksworth and Rossman, 1997; Hawksworth, 2001). At the local scale, fungal diversity is thought to have important consequences for ecosystem functioning particularly through their contribution to nutrient cycling and transport, moisture retention, and plant growth (van der Heijden et al., 1998; Setälä and McLean, 2004). Despite the presumptive great importance of the relationship between fungal diversity and ecosystem function (Nielsen et al., 2011), little is known about extent of the diversity even at the local scale and about the factors generating and maintaining it mainly because the majority of the species are not yet been described and are potentially unculturable. Next-generation amplicon sequencing technique targeting for example intergenic transcribed spacer (ITS)-region have greatly improved our understanding of the community diversity (Meiser et al., 2014; Tedersoo et al., 2014) but this technique offers little insight into function.
Traditionally, fungi have been divided into discrete ecological guilds, such as leaf litter-decomposers, humus saprobes, white-, and brown-rot wood decaying fungi, parasites, and mycorrhizal symbionts (Carlile et al., 2001). However, the functionality of individual species, and the synergistic effects among them are often obscure. For instance, from the well-studied wood-decaying systems it is known that fungal species differ markedly in the way they decompose wood (Worrall et al., 1997), with the decomposition rates depending on the fungal community structure (Boddy et al., 1989). Furthermore, every ecosystem has many fungal types with overlapping functions (Liers et al., 2011).
There is a clear succession in the communities that are involved in the decomposition of litter and organic matter. First, when fresh material is released to the system, r-strategists (both bacteria and fungi) will appear in high numbers, and after the easily degraded material is gone, K-strategists will start the more complicated degradation processes (de Boer et al., 2005). For these steps different sets of genes are required to produce enzymes needed to breakdown the organic material. The first colonizers include yeasts, zygomycetes, and several ascomycetes that produce cellulases and use easily available soluble compounds before being replaced by asco- and basidiomycetes that produce ligninolytic enzymes and glycosyl hydrolases (GHs) to attack the insoluble and more recalcitrant substrates such as lignin. Cellulases play a key role in organic carbon turnover in soil ecosystems (Rabinovich et al., 2002). The ability to degrade cellulose is distributed across the entire fungal kingdom (Lynd et al., 2002) while the ligninolytic enzymes (Wong, 2009) such as laccases, lignin peroxidases, and manganese peroxidases (Martinez et al., 2009; Edwards et al., 2011) are produced only by certain specialized saprotrophic fungi and especially basidiomycetes (Paul, 1989; Liers et al., 2011). Most of the actions are however synergistic: easily degradable cellulose from plant cell walls, an important source of energy for saprotrophic micro-organisms, is bundled into crystalline fibrils, which are linked by hemicellulose to form larger, lignin-encrusted microfibrils (Koshijima and Watanabe, 2003). Thus, a community of fungi and bacteria with different enzymatic capabilities rather than individual species is needed to break down all the cellulose, hemicellulose, and lignin.
Recent studies have confirmed that the content and composition of the gene families encoding enzymes involved in decomposition processes such as fungal oxidative lignin enzymes (FOLy) and carbohydrate-active enzymes (CAZy; Lombard et al., 2014) could be used as functional indicators for fungi (Eastwood et al., 2011; Floudas et al., 2012). Glycosyl hydrolases (EC 3.2.1-) are genes primary involved in cellulose and hemi-cellulose degradation. The division into GH families is based on functionality of the enzymes produced and different clades may be identified through sequences of the functional genes encoding the enzymes and subsequently their prevalence in the environment assessed. Furthermore, although bacteria and archaea are assumed to be the most important microbes participating in nitrogen cycling, fungi are known to contribute to the cycling of nitrogen as well. Fungal contribution to nitrogen cycling in this system was evaluated using nitrate reductase (NIAD) gene composition and fragment richness (Gorfer et al., 2011).
The sequencing of genomes of various species of fungi (http://1000.fungalgenomes.org/) has enabled us to investigate the presence and types of functional genes and subsequently design probes to target as many fungal groups as possible. Despite the increased capacity and decreased price of sequencing, due to the low number of transcripts of fungal genes involved in plant cell wall degradation, making up often <1% of the metatranscriptomic library (Damon et al., 2012), targeted approaches are still needed in ecological studies in order to compare statistically meaningful sample sets with each other. Therefore, the main objectives of this study were to characterize major fungal gene fragments associated with degradation of organic matter in grassland soils and to study if conversion of fields from agricultural used soils to species rich grasslands affects the diversity of selected genes.
To meet these objectives, we developed new primer sets and validated existing ones targeting functional genes and applied these to study the fungal functionality in chronosequence of abandoned agricultural land. The genes included in this study were selected based on existing metagenomics and metatranscriptomic information of soils and are thought to be the types most commonly found in soils (Kellner et al., 2010; Damon et al., 2012). This system has earlier been characterized to have a stable fungal biomass after initial increase after 2 years since abandonment from agriculture (Van der Wal et al., 2006). The fungal community structure has been shown to be different between the short term and long term abandoned fields in an earlier study (Thomson et al., 2015) but the functionality of fungi is not yet studied. The fields included in this study were abandoned for either <10 years or over 25 years. We hypothesized that the fields abandoned for longer time would have more genes coding for enzymes involved in the degradation of recalcitrant organic matter and the recently abandoned fields would have more fungi capable of producing enzymes degrading simple sugars. Furthermore, we hypothesize that the diversity and richness of genes would increase along the secondary succession gradient.
Materials and Methods
Field Sites and Sampling
We selected six fields of which three were abandoned from agricultural use for <10 years [fields Reijerskamp (RE) and Oud Reemst (OR) for 6 years and field Telefoonweg (TW) for 10 years] and three sites for over 25 years ago [Mosselseveld (MV) for 26 years and Boersbos (BB) and Dennenkamp (DE) for 29 years, respectively]. All fields, located in the Veluwe region in the Netherlands, were managed by low-intensive grazing of natural and introduced herbivores (roe deer, fallow deer, red deer, horses, Scottish Highland cattle) to prevent forest development (Van der Wal et al., 2006). Details on the sampling locations, soils, and the soil properties are presented in Table S1.
Three replicate cores of 5 cm diameter and 20 cm deep were collected from the six fields and kept separate. The total number of samples was thus 18. The soil was sieved (4 mm mesh size) and stored at −20°C for molecular analyses. Fresh soil was used to culture fungi in order to obtain pure fungal cultures from each soil type.
DNA Extractions from Soil and Fungal Cultures
DNA from soils was extracted from ~0.5 g of soil (wet weight) with a Power Soil DNA isolation kit (MOBIO Laboratories, Inc.) using a bead beating system. Yields of genomic DNA were inspected on 1% agarose gel and visualized under UV after ethidium bromide staining.
Pure fungal cultures were obtained from sieved fresh soil collected from the same soils by culturing them in three different media. The media used were water agar (WA), 1/10 potato dextrose agar (PDA), and basidiomycete specific media (Blanchette et al., 2010) with added TRITON-X100 (Sigma-Aldrich) to prevent colony expansion. Representative colonies were collected and grown on respective media without added TRITON-X100. The DNA from these fungi was extracted using Zymo Research Fungal/Bacterial DNA kit according to manufacturer's instructions with added β-mercaptoethanol.
Primer Design and Validation
The target genes tested are presented in Table 1. These genes were selected based on their presence in soils and availability of primers. If primers were not found in literature but the gene was found to be common in grassland soils in previous metagenomics studies, an attempt to design new primers was made. In addition to the genes included in the study, two other genes (namely GH1 and GH61) were evaluated for possible primer sites, but as such were not found (due to similarity with bacterial sequences or lack of conservation) these genes were not investigated further. The primers designed in this study are based on 21–49 sequences (depending on the available sequences per gene) of desired gene from fungal origin in CaZY database in November 2012. CaZY database was selected as it is a curated database for carbohydrate genes. In order to reach the best possible coverage for the primers, sequences from all fungal orders represented in CaZY were used in analyses and bacterial and plant genes were used as outgroups. All sequences from desired gene origin were aligned using MEGA4 software and conserved regions were detected visually. Primers were designed to amplify coding regions by using known protein sequences as references. The primers were made as degenerate possible to cover most fungal taxa possible. Gradient PCRs were performed for DNA from pure cultures and soil DNA with the newly developed primers to determine annealing temperatures and if multiple options of primers were designed, best one was selected (Table S2). The criteria for selection of best primer combination were: amplification of as many fungal pure cultures as possible and a clear PCR product from the soil samples. All the primers found from literature were further tested and validated using eight pure cultures of fungi. The optimized PCR conditions and primer sequences are presented in Table S2 and detailed information on primer design in Supplementary File 1. All PCRs were performed with positive (DNA isolated from fungal pure culture) and negative controls.
Table 1. Name used for the gene, enzyme it codes and reference for the primer set used.
PCR products from all 18 soil samples (6 soils and three replicates from each soil) were purified using Qiaqen PCR purification kit and cloned into Escherichia coli JM109 using the pGem-T Easy System II cloning kit (Promega, UK) with a vector: insert ratio of 3:1. Approximately 40 successful transformants per sample were selected for amplification and further 24 clones per soil were sequenced. The selection of clones to be sequenced was based on covering as many types as possible by selecting amplicons of different lengths. If multiple fragments of same size were present, the clones were selected randomly. In total, 144 clones per functional gene were sequenced per sample resulting to total of around 2500 clones sequenced.
Sequence Analyses
Nucleotide sequences in a same sample showing 97% similarity to each other were classified into operational taxonomic units (OTUs) by aligning the sequences with each other (with ClustalW) and calculating their pairwise similarity in MEGA4 (Tamura et al., 2007). This was done to calculate the richness of the unique OTUs per sample. OTUs covering >10% of the library per treatment were defined as dominant. Minimum of eight reference sequences from each gene family were retrieved from CAZY database and trimmed in MEGA4 to correspond the length of sequences obtained in this study (reference strains presented in Table S3). The sequences were aligned with each other using ClustalW (Larkin et al., 2007) to determine their similarities across samples and to the reference sequences. Taxonomic assignment of sequences was performed 75, 85, 90, and 95% for order, class, family and genus designations, respectively (Chen et al., 2013). Unique sequences (aligned across samples in order to determine one representative sequence for each type) were deposited in NCBI GenBank under accession numbers KP748529-KP749176.
In order to evaluate the composition of protein sequences, Augustus software was used to predict the partial protein (and mRNA) sequences from the DNA sequences using fungal genomes as a reference to enable correct annotation of the coding region (Keller et al., 2011). The sequences that did not code a protein when compared with multiple fungal strains were discarded from the mRNA data-set. The similarity of the obtained proteins was compared to known protein sequences obtained from the CaZy database in a similar way than was done for nucleotide sequences.
Statistical Analyses
UniFrac-weighted PCA (weighted using number of OTUs in a sample) was applied to examine differences in the functional communities between soil samples (Lozupone and Knight, 2005) using the data on unique OTUs in a sample and the information on how many times this sequence was present in the sample. Trees were built using ClustalW aligned data and exported to UniFrac.
The dissimilarities in gene composition between samples were calculated based on their UniFrac distance and correlated with each other and with dissimilarities between samples in soil physio-chemical factors measured, using linear correlation. Dissimilarities in plant communities which was a multivariate data-set were calculated using Principle Component Analyses (PCA) and calculating the dissimilarities between treatment for the two main axis. The two first axis explaining together 79% of the total variation in plant cover were used to calculate the dissimilarities in plant community between samples. The differences in richness of the genes and genes assigned into function were evaluated using the one-way analyses of variance testing. The assumption of normality was tested with Shapiro–Wilk statistics and homogeneity of variances was assessed with Levene's test. Differences between short-term and long-term abandoned fields were tested with Tukey's HSD test.
Results
The primers in this study were designed to catch as broad spectrum of fungi as possible (see Supplementary File 1 for details). For example, the cellobiose dehydrogenase primers amplified in silico 25 of the 29 selected reference strains from 15 (14) orders. Only Stachybotrys bisbyi (Hypocreales), Grifola frondosa (Agaricales), and two Verticillium species (Glomerellales) had mismatches in their cellobiose dehydrogenase gene which could prevent primer attachment. In vitro, these primers amplified five out of the eight fungal strains tested. The selected glucose oxidase primers amplified in silico 19 out of the selected 21 strains and in vitro seven out of eight fungi. The details on the performance of all the primers designed can be found in the Supplementary File 1. All existing primers were tested with the same eight pure cultures of soil fungi for their specificity (Table S3). Some primers amplified all the selected cultured fungi while some were specific to groups or phyla. Details on the specificity can be found in Table S3.
After optimization of methods, the PCRs for the soils gave positive signal for all except one gene tested. Unfortunately we could not obtain enough PCR products from three of the soils (MV, TW, and BB) for successful cloning of the endoglucanase gene GH45 and this gene was excluded from further analyses. The specificity of the other primers varied and the percentage of sequences that could be correctly assigned to the desired gene family was between 62.5% for oxalate decarboxylase to 100% for seven gene fragment types (GH5, GH18, GH31, GH51, GH74, GH76, and glucose oxidase). The percentage of sequences predicted by Augustus to code a protein sequence was between 52.5% for GH5 and 88.3% for nitrate reductase genes. The sequences that did not code a (partial) protein were removed from the dataset and not included further in the analysis.
The number of unique sequence types obtained per sample ranged from a single dominant type of GH51 fragment in field RE to 18 types of GH18 fragments in field DE, respectively (Table 2). There were differences in the number of unique types between soils but the richness of none of the genes studied was significantly affected by the time since abandonment (Table 2). However, when looking at groups of genes together clustered by function, we observed more genes involved in degradation of lignin, cellulose and hemicellulose in fields abandoned for longer than 25 years [Figure 1; F(1, 4) = 9.3, p = 0.038].
Table 2. The number of unique fragments of each gene studied detected in each field and averaged over treatments.
Figure 1. The number of unique gene fragments determined as number of unique sequences in a sample in long-term (blue) and short-term (red) abandoned fields. The bars represent the averages of 3 fields per treatment and error bars depict standard deviation. The list of how genes are divided into categories is presented in Table 2. Statistically significant differences (one-way ANOVA) between long-term and short-term abandoned fields are marked with asterisks.
As measured by using the calculated UniFrac significant differences, the short-term abandoned fields differed significantly from the long-term abandoned fields in the composition of GH3, GH6, GH10, GH31, heme-thiolate peroxidases, and nitrate reductases (Table 3). In most cases species composition in an individual fields was unique to the field and phylogenetically similar gene-types were not shared within abandonment category (Table 3, Figure 3). The relationship of the samples from the short term and long-term abandoned fields is reflected in the PCA plots (Figure 3).
Table 3. UniFrac significances between treatments based on both nucleotide and protein distances.
Due to limited numbers of sequenced clones in most of the gene families and the lack of typical grassland species in sequenced genomes, it was impossible to assign sequences from environmental samples to the genus level. In many cases, however, it was possible to assign the sequences to a certain order (with similarities of >75% to known sequenced species). The mean BLAST N (based on DNA) and BLAST P scores here varied between 55–100% and 38–95% identities, respectively. After aligning the sequences obtained in this study with selected sequences from the databases representing as diverse group of fungi as possible (Table S4), most similar reference sequence was resolved for each sequence indicating roughly to which phyla and order the type belongs to. The highest similarities with known sequences of the most common sequence type in each field are presented in Table 4. The explained differences in GH3 gene composition could be explained with the absence of sequences related to known GH3 β-glucosidases from Agaricales in short term abandoned fields where sequences related to Sordariomycetes were dominant (Figure 3). The separation of one field observed for GH5 could be explained by dominance of a type related to exo-β -1,3-glucanase from Laccaria bicolor (86–96% similarity) in all except that field. The dominant GH6 fragment type in all the long-term abandoned fields was related to a GH6 gene from Eurotiales while the most common type in short-term abandoned fields was related to type identified earlier from species belonging to Sordariales. The significant differences between short-term and long-term abandoned fields were explained by absence of few cellobiohydrolases from basidiomycete origin in short term abandoned fields. The GH10 types dominating the long-term abandoned fields were most related to basidiomycete fungus and in short-term abandoned fields dominant types were related to ascomycete xylanase explaining observed differences between the treatments (Figure 3). The GH31 types responsible for significant differences observed between short-term and long-term abandoned fields (Figure 3) were attributed to ascomycetes Myceliophthora thermophile and Pyrenophora tritici-repentis, and basidiomycete Agaricus bisporus, all present only in long-term abandoned fields. The GH51 types explaining the pattern seen in Figure 2 were types related to Magnaporthales sp. and to Penicillium sp. present only in short-term abandoned field OR. The dominant GH51 types in long-term abandoned fields were related to gene for α-L-arabinofuranosidase from Pleurotus sp. and in short-term abandoned fields α -L-arabinofuranosidase from Trichoderma.
Table 4. The dominant species (and similarity %) in the fields.
Figure 2. The network of connected edaphic and biological variables based on correlation of dissimilarity matrix. The cellulose and hemicellulose related genes are marked with red balls, nitrogen uptake genes with green, lignolytic genes with blue, glucose oxidase in purple, edaphic factors with yellow, and turquois and plant cover (calculated from PCA) in dark green. Only significant correlations are presented and the r value of the linear correlation is shown next to the connecting lines.
The dominant types of laccases were from basidiomycete origin but also fragments similar to laccases from ascomycete species were detected. Field RE was significantly different from other fields mostly due to absence of Neurospora, Aspergillus, and Rhizoctonia type fragments which led to dominance of type similar to basidiomycetes (Polyporales). Mn-peroxidases in all fields were dominated with fragments from basidiomycete origin. Field MV separated due to dominance of a type similar to uncultured Cortinarius clone from a previous study and presence of a type categorized to Russulales (83% similarity). The dominant types of heme-thiolate peroxidases were categorized as Capnodiales in most of the fields. Both fields BB and RE had many unique types that were similar to sequences from a previous study but not closely related to known sequences from any organism (Kellner et al., 2010). The dominant types of nitrate reductases were similar to nitrate reductases from P. chrysosporium and F. oxysporum and the dominant types were not significantly different between fields. The observed differences in community structure could be explained with presence and absence of rarer gene types such as uncultured Pezizomycotina nitrate reductase (HQ243641) that was primary present in fields abandoned for short time.
There were significant differences in functional gene composition between the fields which could not be explained by time since agriculture (Figure 3). Thus, the dissimilarity of the fields to each other were calculated using UniFrac distances and correlated with dissimilarities between measured environmental parameters. Significantly correlated dissimilarities are connected with lines or overlapping balls between genes and environmental parameters in the network Figure 2. Furthermore, to learn about connections between the genes, the dissimilarities were correlated with each other and significant correlations were marked with lines connecting the genes (Figure 2). There were total of 19 significant dissimilarity correlations between gene fragments (Figure 2). Eight of these correlations were negative and 11 positive. All other genes except GH18 family chitinases were correlated with at least one other gene or environmental factor. Cellobiohydrolase II gene (GH6) fragment dissimilarity was negatively correlated with cellobiohydrolase I (GH7) fragment dissimilarity. Endo-β-1,4-xylanase dissimilarities of families GH10 and GH11 were also negatively correlated with each other.
Figure 3. PCA ordination of the functional community structure in long-term (blue) and short-term (red) abandoned fields based on weighted UniFrac distances. Significance is based on UniFrac significance between treatments, is marked in the figure and significantly different treatments/samples are circled.
Mn-peroxidase and laccase fragment composition dissimilarity between fields was positively correlated with the dissimilarity in pH and associated (metal) ion contents in soils (Figure 2). pH differences in soils further affected the basidiomycete community responsible for nitrate uptake. Nitrogen content dissimilarity in soils (calculated dissimilarities between soils in their nitrogen content) was correlated with dissimilarity of the fungal genes responsible for nitrate uptake. The nitrate reductases from basidiomycetes were positively correlated with dissimilarity of nitrogen contents in soils while ascomycete nitrate reductases were negatively correlated. GH5 and GH3 were negatively and GH76 positively associated with soil phosphorous dissimilarity in soils (Figure 2). Dissimilarity in plant cover between fields calculated from principal component (PC) axis explaining 79% of the variation in plant community (44 and 35% for PCA1 and PCA2, respectively). The plant cover expressed using the PCA axis was related to GH10, Heme-thiolate peroxidase, and Mn-peroxidase composition dissimilarities in the soils.
Discussion
In contrast to the commonly used analyses of soil fungi using phylogenetic markers such as ITS genes (Lindahl et al., 2013), we investigated the presence and composition of 20 functional genes encoding mostly carbon cycling related and one nitrogen related genes. To achieve a goal to evaluate functional shifts in soil fungi along a chronosequence of soils abandoned from agriculture, new primers were developed and existing primers tested. An average of 138 unique sequences derived from 19 functional genes was obtained per soil sample which makes it higher coverage than obtained in metatransctiptomics or metagenomics studies. For example in a study by Damon et al. (2012) maximum of six types of genes from any given GH family were found in a sample. The genes studied here were all pre-selected for their presence and importance in soils and are involved in crucial soil processes such as cellulose, hemicellulose, and lignin degradation and nitrate reduction in which fungi are key agents. By selecting a broader range of genes as targets, we could inspect the community richness and function as a whole instead of focusing on only one or few genes of interest (Edwards et al., 2008; Kellner et al., 2009).
Unfortunately, the majority of the types of genes detected in this study could not be linked to a known species. This is an expected result as public databases have limited number, often <50, of correctly annotated reference genes and further reinforces the idea that identities of most of the genes from soil fungi are yet unknown: sequences produced in this study with highest similarities to known sequences were GH7 fragments and this was only due to the large amount of environmental sequences of cellobiohydrolases in databases. When using an approach of aligning gene fragments (mRNA) from this study to annotated genes (mRNA) from sequenced species, the similarities were in average around 85%. These results are in the range of other studies finding that average BLAST P scores obtained were around 80% for GH7 cellobiohydrolase (e.g., Edwards et al., 2008). Here, using cloning and subsequent Sanger-sequencing it was not needed to select for certain sized fragments allowing us to capture gene fragments from as diverse origin as possible. This is an advantage of this technique compared to for example some next generation sequencing (NGS) techniques such as Illumina HiSeq or MiSeq, and pyrosequencing. These techniques will generate more sequences but only fairly short and uniformly sized amplicons can be sequenced (Goodwin et al., 2016). Here, we detected large variation in intron positions and resulting amplicon sizes making the techniques sensitive to this variation less suitable for this purpose.
In some cases, the primer sets are selective to certain group of organisms. Furthermore, there is a need for primer sets that amplify fungi but exclude bacterial and plant derived sequences. Some of the primers used in this study are known to preferentially amplify only ascomycetes (GH3 and GH76) or basidiomycetes (i.e., Mn-peroxidases; Kellner et al., 2010) and our data concurs with these view (Table S3). This is primarily due to the natural lack of conservation between orders and phyla and technical design issues such as avoidance of degeneracy and limitations in the amount of reference sequences (Kellner et al., 2010). Yet, we feel that the primer sets designed in this study will enable further characterization of unknown fungal genes in soils and the long sequences generated here enable further primer development for the further capture of larger part of the soil functional biodiversity.
The depth of analyses for individual genes was rather low in this study especially when considering the diversity of soil fungi detected using analysis of ITS region (Thomson et al., 2015) and it is likely that some more quiescent groups of fungi might have not been detected. Nevertheless, the strength of this study is the multigene approach proposed and we are confident that we have captured the dominant fungal functional genes in these soils and increased, in some cases even doubled, the amount of sequences derived from soil and deposited in databases in the case of all the genes. Using this integrated multi-gene approach, we showed that the total number of fragments identified was affected by the time since abandonment from agriculture (Figure 1). The long-term abandoned fields had higher total number of unique fragments when all genes involved in cellulose, hemicellulose, and lignin degradation studied were calculated together. This could indicate that these soils abandoned for over 25 years have a larger pool of potential functions and that this functional diversity would increase in time since agriculture. Earlier study using ITS-region as marker found no significant differences between recent and long-term abandoned fields in their alpha-diversity of fungi (Thomson et al., 2015). This is in accordance with results presented here when single marker genes were used (Table 2) but contrary to the results obtained using a multi-gene approach (Figure 1).
By using ITS region as a marker, difference in community structure between short- and long-term abandoned fields was detected (Thomson et al., 2015). Here we could show that the community composition of GH3, GH6, GH10, GH31, heme-thiolate peroxidase, and nitrate reductase genes was different between short- and long-term abandoned fields making them good candidates for functional screening of soils showing differences corresponding to ITS-region. The reason why we detected changes in these genes but not in the community structure of GH5 and GH11 for example, is unclear. One possible explanation is that different fungal species possess different sets of CaZymes (Veneault-Fourrey et al., 2014) and that the diversity of the community producing certain enzymes is very small (i.e., only basidiomycetes), making it likely that differences in the community structure do not exist or are not large enough to be detected. For example, laccases can be produced by asco- and basidiomycetes and some bacteria (i.e., actinomycetes) while Mn-peroxidases are produced only by basidiomycetal fungi (Baldrian, 2006; Bödeker et al., 2009). Surprisingly, in our soils the cellobiose dehydrogenases were mostly ascomycetal types and involved in cellulose decomposition even though earlier studies claiming them to be predominantly basidiomycetal and involved also in lignin degradation (Harreither et al., 2011). This finding is related to the use of new primers which target broader range of fungi and especially ascomycetes and highlights the importance of development of primers to target broader groups of fungi in order to gain better understanding of the function of the soil fungal communities.
For GH3 and GH10, we saw a clear shift in communities in phylum level: in short-term abandoned soils the dominant types of this gene were related to ascomycetes while long-term abandoned fields were dominated with basidiomycetal types. In case of GH6 the corresponding shift explaining the observed differences between treatments was within Ascomycota (i.e., from Eurotiales to Sordariales). These are major shifts and partly corroborated with the changes observed in the taxonomic study of the same soils (Thomson et al., 2015). Probably these are a relic of the supply of substrates, but as this is not measured here, we can only speculate. Earlier studies have found that GH7 composition or richness in soil or litter is not significantly affected by tree species (Edwards et al., 2008), genetic modification of potato leaves (Hannula et al., 2013), elevated CO2 (Hassett et al., 2009), or N fertilization (Weber et al., 2012). Our results corroborate the lack of differentiation power of the GH7 gene to environmental disturbances. Based on evidence gathered here, we recommend that instead of focusing on one gene, studies should include multiple genes in the analyses to determine which ones are the most responsive to the disturbance in question.
The difference observed in heme-thiolate peroxidase composition can be explained by the observed differences in plant cover: fresh plant litter is rich in phenolics with high solubility and this varies greatly between plant species (Meier et al., 2008) and heme-thiolate peroxidases are involved in breakdown of these aromatic compounds while the other lignolytic enzymes are used to degrade lignin. Here we observed no differences in the composition of fungal genes involved in lignin degradation between the short-term and long-term abandoned fields. Earlier studies have found large differences in the presence and activity of genes involved in lignin decomposition between sampling times (Kellner et al., 2009), along soil depth gradients (Snajdr et al., 2008a) and between sampling sites (Snajdr et al., 2008b) in forest soils. The difference between our study sites and the earlier studies is that we were investigating grassland soils where the importance of lignocellulolytic genes might be overshadowed by the cellulolytic activities and the so-called r-strategist fungi.
By adopting an approach to look at the dissimilarity correlations between the genes and environmental factors, we could tease out the connections between the genes studied and edaphic factors (Figure 2). It is known that soil characteristics affect laccase gene composition (Lauber et al., 2009; Chen et al., 2013) and especially pH is known to influence the activity of ligninolytic enzymes in soils (Fujii et al., 2013). Furthermore, in pure culture studies it has been shown that ligninolytic enzymes have different ranges of optimum pH (Wong, 2009). We acknowledge that the use of DNA based methods tells only about the presence of these genes and not about actual production/activity of the enzyme but we are confident based on earlier research that these factors correspond well (Chen et al., 2013). Here, we detected a strong dependence of both laccases and Mn-peroxidases on the pH of the soils (Figure 2) and the conclusion may, thus, be drawn that the soil pH is the most important factor governing the composition of ligninolytic enzymes in these soils. The composition of GH3, GH5, and GH76 genes seemed to respond to differences in soil phosphorous content while, for example, sodium dissimilarities were correlated with cellobiose dehydrogenases. Logically, the nitrogen uptake genes were related to soil nitrogen content and could thus be used to study the effects of e.g., nitrogen fertilization on soil communities (Gorfer et al., 2011). Another factor affecting basidiomycetal nitrogen is, again, soil pH (Figure 2). Indeed, we can conclude that processes dominated by basidiomycetes were strongly affected by the soil pH while other soil and plant related factors were shaping the predominantly ascomycetal processes. Interestingly, nitrogen uptake was one the genes/process affected by time since abandonment and was tightly linked with dissimilarities in nitrogen contents in the soils. No gradient of nitrogen content was present in the soils but it has been hypothesized that the N mineralization is more important in later successional stages (Holtkamp et al., 2011) which coincides with increased importance of fungal communities in this system.
Fungi use a multitude of enzymes in synergy to break down plant cell walls and thus connections among genes are to be expected. For example, oxidoreductive enzymes such as cellobiose dehydrogenase can act synergistically with canonical GHs, accelerating the breakdown of litter (Langston et al., 2011). On the other hand, we saw that also genes producing similar enzymes (such as cellobiohydrolases GH6 and GH7) group together (Figure 2). The gene family GH18 involved in production of fungal chitinases was the only family not related to any edaphic factor or other gene. This is in accordance with the assumption that this gene may be involved in degradation of fungal hyphae and not involved in the decomposition of plant derived material (Hartl et al., 2012). Surprisingly, oxalate decarboxylase composition was not related to lignin degradation pathways even though it is known to be involved in regulation of intra- and extracellular quantities of oxalic acid, which is one of the key components in biological decomposition lignin (Mäkelä et al., 2009). In this system the composition of oxalate decarboxylases seems to, however, be linked to GHs from families 10 and 74 and to nitrogen uptake by Ascomycota. The cellobiose dehydrogenenase was also grouping with cellulose degradation process instead of lignin decomposition despite its known involvement in both processes (Hyde and Wood, 1997).
Using this targeted approach, we obtained over 100 unique sequences involved in carbon and nitrogen cycling per soil type which we used to characterize both the richness and community composition of the soils. The main finding is that different processes in soils are driven by different edaphic factors (such as pH) and that conversion of former agricultural land into species rich grassland affects (via changes in plant cover) the cellulolytic fungal community rather than the community that is involved in the decomposition of recalcitrant materials such as lignin.
Author Contributions
SEH and JAvV designed the experiments and wrote the manuscript. SEH performed the experiments and did the statistical analysis.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
This study was carried out as part of the EcoFINDERS research project (EU- FP7–264465). We thank Elly Morrien for help in field sampling and Wim van der Putten for comments on the manuscript. This is NIOO-KNAW publication 6201.
Supplementary Material
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmicb.2016.01897/full#supplementary-material
References
Baldrian, P. (2006). Fungal laccases - occurrence and properties. FEMS Microbiol. Rev.30, 215–242. doi: 10.1111/j.1574-4976.2005.00010.x
Blanchette, R. A., Held, B. W., Arenz, B. E., Jurgens, J. A., Baltes, N. J., Duncan, S. M., et al. (2010). An Antarctic hot spot for fungi at Shackleton's historic hut on Cape Royds. Microb. Ecol. 60, 29–38. doi: 10.1007/s00248-010-9664-z
Boddy, L., Owens, E. M., and Chapela, I. H. (1989). Small-scale variation in decay-rate within logs one year after felling - effect of fungal community structure and moisture-content. FEMS Microbiol. Lett. 62, 173–184. doi: 10.1111/j.1574-6968.1989.tb03691.x
Bödeker, I. T., Nygren, C. M., Taylor, A. F., Olson, A., and Lindahl, B. D. (2009). ClassII peroxidase-encoding genes are present in a phylogenetically wide range of ectomycorrhizal fungi. ISME J. 3, 1387–1395. doi: 10.1038/ismej.2009.77
Carlile, M. J., Watkinson, S. C., and Gooday, G. W. (2001). The Fungi. San Diego, CA: Academic Press.
Chen, X. B., Su, Y. R., He, X. Y., Liang, Y. M., and Wu, J. S. (2013). Comparative analysis of basidiomycetous laccase genes in forest soils reveals differences at the cDNA and DNA levels. Plant Soil 366, 321–331. doi: 10.1007/s11104-012-1440-z
Damon, C., Lehembre, F., Oger-Desfeux, C., Luis, P., Ranger, J., Fraissinet-Tachet, L., et al. (2012). Metatranscriptomics reveals the diversity of genes expressed by eukaryotes in forest soils. PLoS ONE 7:e28967. doi: 10.1371/journal.pone.0028967
de Boer, W., Folman, L. B., Summerbell, R. C., and Boddy, L. (2005). Living in a fungal world: impact of fungi on soil bacterial niche development. FEMS Microbiol. Rev. 29, 795–811. doi: 10.1016/j.femsre.2004.11.005
Eastwood, D. C., Floudas, D., Binder, M., Majcherczyk, A., Schneider, P., Aerts, A., et al. (2011). The plant cell wall-decomposing machinery underlies the functional diversity of forest fungi. Science 333, 762–765. doi: 10.1126/science.1205411
Edwards, I. P., Upchurch, R. A., and Zak, D. R. (2008). Isolation of fungal cellobiohydrolase I genes from sporocarps and forest soils by PCR. Appl. Environ. Microbiol. 74, 3481–3489. doi: 10.1128/AEM.02893-07
Edwards, I. P., Zak, D. R., Kellner, H., Eisenlord, S. D., and Pregitzer, K. S. (2011). Simulated atmospheric N deposition alters fungal community composition and suppresses ligninolytic gene expression in a Northern Hardwood Forest. PLoS ONE 6:e20421. doi: 10.1371/journal.pone.0020421
Floudas, D., Binder, M., Riley, R., Barry, K., Blanchette, R. A., Henrissat, B., et al. (2012). The paleozoic origin of enzymatic lignin decomposition reconstructed from 31 fungal genomes. Science 336, 1715–1719. doi: 10.1126/science.1221748
Fujii, K., Uemura, M., Hayakawa, C., Funakawa, S., and Kosaki, T. (2013). Environmental control of lignin peroxidase, manganese peroxidase, and laccase activities in forest floor layers in humid Asia. Soil Biol. Biochem. 57, 109–115. doi: 10.1016/j.soilbio.2012.07.007
Goodwin, S., McPherson, J. D., and McCombie, W. R. (2016). Coming of age: ten years of next generation sequencing technologies. Nat. Rev. Genet. 17, 333–351. doi: 10.1038/nrg.2016.49
Gorfer, M., Blumhoff, M., Klaubauf, S., Urban, A., Inselsbacher, E., Bandian, D., et al. (2011). Community profiling and gene expression of fungal assimilatory nitrate reductases in agricultural soil. ISME J. 5, 1771–1783. doi: 10.1038/ismej.2011.53
Hannula, S. E., de Boer, W., Baldrian, P., and Van Veen, J. A. (2013). Effects of genetically modified amylopectin-accumulating potato in decomposer processes and fungal diversity in litter and soil. Soil Biol. Biochem. 58, 88–98. doi: 10.1016/j.soilbio.2012.11.008
Harreither, W., Sygmund, C., Augustin, M., Narciso, M., Rabinovich, M. L., Gorton, L., et al. (2011). Catalytic properties and classification of cellobiose dehydrogenases from Ascomycetes. Appl. Environ. Microbiol. 77, 1804–1815. doi: 10.1128/AEM.02052-10
Hartl, L., Zach, S., and Seidl-Seiboth, V. (2012). Fungal chitinases: diversity, mechanistic properties and biotechnological potential. Appl. Microbiol. Biotechnol. 93, 533–543. doi: 10.1007/s00253-011-3723-3
Hassett, J. E., Zak, D. R., Blackwood, C. B., and Pregitzer, K. S. (2009). Are basidiomycete laccase gene abundance and composition related to reduced lignolytic activity under elevated atmospheric NO3- Deposition in a Northern Hardwood Forest? Microb. Ecol. 57, 728–739. doi: 10.1007/s00248-008-9440-5
Hawksworth, D. L. (2001). The magnitude of fungal diversity: the 1.5 million species estimate revisited. Mycol. Res. 105, 1422–1432. doi: 10.1017/S0953756201004725
Hawksworth, D. L., and Rossman, A. Y. (1997). Where are all the undescribed fungi? Phytopathology 87, 888–891. doi: 10.1094/PHYTO.1997.87.9.888
Holtkamp, R., van der Wal, A., Kardol, P., van der Putten, W. H., de Ruiter, P. C., and Dekker, S. C. (2011). Modelling C and N mineralisation in soil food webs during secondary succession on ex-arable land. Soil Biol. Biochem. 43, 251–260. doi: 10.1016/j.soilbio.2010.10.004
Hyde, S. M., and Wood, P. M. (1997). A mechanism for production of hydroxyl radicals by the brown-rot fungus coniophora puteana: Fe(III) reduction by cellobiose dehydrogenase and Fe(II) oxidation at a distance from the hyphae. Microbiology 143, 259–266. doi: 10.1099/00221287-143-1-259
Keller, O., Kollmar, M., Stanke, M., and Waack, S. (2011). A novel hybrid gene prediction method employing protein multiple sequence alignments. Bioinformatics 27, 757–763. doi: 10.1093/bioinformatics/btr010
Kellner, H., Luis, P., Schlitt, B., and Buscot, F. (2009). Temporal changes in diversity and expression patterns of fungal laccase genes within the organic horizon of a brown forest soil. Soil Biol. Biochem. 41, 1380–1389. doi: 10.1016/j.soilbio.2009.03.012
Kellner, H., Zak, D. R., and Vandenbol, M. (2010). Fungi unearthed: transcripts encoding lignocellulolytic and chitinolytic enzymes in forest soil. PLoS ONE 5:e10971. doi: 10.1371/journal.pone.0010971
Koshijima, T., and Watanabe, T. (2003). Association Between Lignin and Carbohydrates in Wood and Other Plant Tissues. Berlin: Springer.
Langston, J. A., Shaghasi, T., Abbate, E., Xu, F., Vlasenko, E., and Sweeney, M. D. (2011). Oxidoreductive cellulose depolymerization by the enzymes cellobiose dehydrogenase and glycoside hydrolase 61. Appl. Environ. Microbiol. 77, 7007–7015. doi: 10.1128/AEM.05815-11
Larkin, M. A., Blackshields, G., Brown, N. P., Chenna, R., McGettigan, P. A., McWilliam, H., et al. (2007). Clustal W and Clustal X version 2.0. Bioinformatics 23, 2947–2948. doi: 10.1093/bioinformatics/btm404
Lauber, C. L., Sinsabaugh, R. L., and Zak, D. R. (2009). Laccase gene composition and relative abundance in oak forest soil is not affected by short-term nitrogen fertilization. Microb. Ecol. 57, 50–57. doi: 10.1007/s00248-008-9437-0
Liers, C., Arnstadt, T., Ullrich, R., and Hofrichter, M. (2011). Patterns of lignin degradation and oxidative enzyme secretion by different wood- and litter-colonizing basidiomycetes and ascomycetes grown on beech-wood. FEMS Microbiol. Ecol. 78, 91–102. doi: 10.1111/j.1574-6941.2011.01144.x
Lindahl, B. D., Nilsson, R. H., Tedersoo, L., Abarenkov, K., Carlsen, T., Kjøller, R., et al. (2013). Fungal community analysis by high-throughput sequencing of amplified markers – a user's guide. New Phytol. 199, 288–299. doi: 10.1111/nph.12243
Lombard, V., Golaconda Ramulu, H., Drula, E., Coutinho, P. M., and Henrissat, B. (2014). The carbohydrate-active enzymes database (CAZy) in 2013. Nucleic Acids Res. 42, 21. doi: 10.1093/nar/gkt1178
Lozupone, C., and Knight, R. (2005). UniFrac: a new phylogenetic method for comparing microbial communities. Appl. Enivron. Microbiol. 71, 8228–8235. doi: 10.1128/AEM.71.12.8228-8235.2005
Lynd, L. R., Weimer, P. J., van Zyl, W. H., and Pretorius, I. S. (2002). Microbial cellulose utilization: fundamentals and biotechnology. Microbiol. Mol. Biol. Rev. 66, 506. doi: 10.1128/MMBR.66.3.506-577.2002
Mäkelä, M. R., Hilden, K., Hatakka, A., and Lundell, T. K. (2009). Oxalate decarboxylase of the white-rot fungus Dichomitus squalens demonstrates a novel enzyme primary structure and non-induced expression on wood and in liquid cultures. Microbiology 155, 2726–2738. doi: 10.1099/mic.0.028860-0
Martinez, D., Challacombe, J., Morgenstern, I., Hibbett, D., Schmoll, M., Kubicek, C. P., et al. (2009). Genome, transcriptome, and secretome analysis of wood decay fungus Postia placenta supports unique mechanisms of lignocellulose conversion. Proc. Natl. Acad. Sci. U.S.A. 106, 1954–1959. doi: 10.1073/pnas.0809575106
Meier, C. L., Suding, K. N., and Bowman, W. D. (2008). Carbon flux from plants to soil: roots are a below-ground source of phenolic secondary compounds in an alpine ecosystem. J. Ecol. 96, 421–430. doi: 10.1111/j.1365-2745.2008.01356.x
Meiser, A., Bálint, M., and Schmitt, I. (2014). Meta-analysis of deep-sequenced fungal communities indicates limited taxon sharing between studies and the presence of biogeographic patterns. New Phytol. 201, 623–635. doi: 10.1111/nph.12532
Nielsen, U. N., Ayres, E., Wall, D. H., and Bardgett, R. D. (2011). Soil biodiversity and carbon cycling: a review and synthesis of studies examining diversity–function relationships. Eur. J. Soil Sci. 62, 105–116. doi: 10.1111/j.1365-2389.2010.01314.x
Rabinovich, M. L., Melnik, M. S., and Boloboba, A. V. (2002). Microbial cellulases (Review). Appl. Biochem. Microbiol. 38, 305–321. doi: 10.1023/A:1016264219885
Setälä, H., and McLean, M. A. (2004). Decomposition rate of organic substrates in relation to the species diversity of soil saprophytic fungi. Oecologia 139, 98–107. doi: 10.1007/s00442-003-1478-y
Snajdr, J., Valaskova, V., Merhautova, V., Cajthaml, T., and Baldrian, P. (2008a). Activity and spatial distribution of lignocellulose-degrading enzymes during forest soil colonization by saprotrophic basidiomycetes. Enzyme Microb. Biotechnol. 43, 186–192. doi: 10.1016/j.enzmictec.2007.11.008
Snajdr, J., Valásková, V., Merhautová, V., Herinková, J., Cajthaml, T., and Baldrian, P. (2008b). Spatial variability of enzyme activities and microbial biomass in the upper layers of Quercus petraea forest soil. Soil Biol. Biochem. 40, 2068–2075. doi: 10.1016/j.soilbio.2008.01.015
Tamura, K., Dudley, J., Nei, M., and Kumar, S. (2007). MEGA4: Molecular Evolutionary Genetics Analysis (MEGA) software version 4.0. Mol. Biol. Evol. 24, 1596–1599. doi: 10.1093/molbev/msm092
Tedersoo, L., Bahram, M., Põlme, S., Kõljalg, U., Yorou, N. S., Wijesundera, R., et al. (2014). Global diversity and geography of soil fungi. Science 346, 1078–1087. doi: 10.1126/science.1256688
Thomson, B. C., Tisserant, E., Plassart, P., Uroz, S., Griffiths, R. I., Hannula, S. E., et al. (2015). Soil conditions and land use intensification effects on soil microbial communities across a range of European field sites. Soil Biol. Biochem. 88, 403–413. doi: 10.1016/j.soilbio.2015.06.012
van der Heijden, M. G. A., Klironomos, J. N., Ursic, M., Moutoglis, P., Streitwolf-Engel, R., Boller, T., et al. (1998). Mycorrhizal fungal diversity determines plant biodiversity, ecosystem variability and productivity. Nature 396, 69–72. doi: 10.1038/23932
Van der Wal, A., Van Veen, J. A., Smant, W., Boschker, H. T. S., Bloem, J., Kardol, P., et al. (2006). Fungal biomass development in a chronosequence of land abandonment. Soil Biol. Biochem. 38, 51–60. doi: 10.1016/j.soilbio.2005.04.017
Veneault-Fourrey, C., Commun, C., Kohler, A., Morin, E., Balestrini, R., Plett, J., et al. (2014). Genomic and transcriptomic analysis of Laccaria bicolor CAZome reveals insights into polysaccharides remodelling during symbiosis establishment. Fungal Genet. Biol. 72, 168–181. doi: 10.1016/j.fgb.2014.08.007
Weber, C. F., Balasch, M. M., Gossage, Z., Porras-Alfaro, A., and Kuske, C. R. (2012). Soil fungal cellobiohydrolase i gene (cbhi) composition and expression in a loblolly pine plantation under conditions of elevated atmospheric CO2 and nitrogen fertilization. Appl. Environ. Microbiol. 78, 3950–3957. doi: 10.1128/AEM.08018-11
Wong, D. W. (2009). Structure and action mechanism of ligninolytic enzymes. Appl. Biochem. Biotechnol. 157, 174–209. doi: 10.1007/s12010-008-8279-z
Keywords: fungi, soil fungal community, functional genes, glycosyl hydrolases, soil microbiology, peroxidases, chronosequence
Citation: Hannula SE and van Veen JA (2016) Primer Sets Developed for Functional Genes Reveal Shifts in Functionality of Fungal Community in Soils. Front. Microbiol. 7:1897. doi: 10.3389/fmicb.2016.01897
Received: 07 April 2016; Accepted: 11 November 2016;
Published: 29 November 2016.
Edited by:
Jeanette M. Norton, Utah State University, USAReviewed by:
Vijai Kumar Gupta, NUI Galway, IrelandAsuncion De Los Ríos, Spanish National Research Council, Spain
David B. Wilson, Cornell University, USA
Copyright © 2016 Hannula and van Veen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: S. Emilia Hannula, e.hannula@nioo.knaw.nl