Peatlands worldwide, particularly northern (boreal and subarctic) parts are shown to be important participant in global carbon (C) cycle in the recent past [1]. Despite covering only 6-8% of the terrestrial ecosystems, northern peatlands store ~550 Pg C [2], which accounts for between one-quarter and one-third of the world’s soil carbon [3]. The sequestration of C arises of northern peatlands is a result of high productivity rates rather than low decomposition. In particular, the reason for a peatland functioning as a C sink is that its vegetation fixes more C than its C lost through outflow of dissolved organic C and emissions of CO₂ and CH₄ [4-6]. However, there is a growing concern that the on-going and future global warming will change the C cycling in these ecosystems, peatlands may return the previously captured C to the atmosphere via releasing CO₂and/or CH₄ , which would possibly accelerate the present warming [7-10].

Many studies have shown that microorganisms play a crucial role in the C cycling process [11-14]. The microbial communities in the peatlands of Europe, America, Canada and UK have been studied [15-17]. However, there are relatively few studies on the peatlands in more temperate parts. The Sanjiang Plain, located in the temperate climate region, is the largest area of freshwater marshlands in China [18]. The peatlands there could be more vulnerable to global warming because they are mainly located at the southern limit of northern peatlands [19]. Thus, it is imperative that we deepen our understanding of the microorganisms in this ecosystem.

Both fungi and bacteria have important functional roles in peat biogeochemical process. Bacteria are more competitive in anoxic soil environments, as they have the ability to utilize alternative electron acceptors (SO₄^2-, NO₃^-) beyond oxygen and simple organic molecules [20,21]. In contrast, fungi typically have relatively lower requirements of biomass N and other nutrients. They are capable of contributing to nutrient cycling within their plant host, as they produce the extracellular enzyme machinery required to breakdown complex plant polymers, including phenolic compounds [22]. However, there are few studies of the fungal community in the Sanjiang Plain, and they have suffered from the lack of a high-quality curated database for taxonomic assignment.

In the present study, high-throughput Illumina sequencing of ITS rRNA genes is used to study the fungal communities in the Sanjiang Plain, the southern edge of northern peatlands. Peat soils were collected from three fens which were started to develop during different periods in this area. The greater sequencing depth achieved by the high-throughput sequencing allows capture of the less abundant and uncultured taxa, thus will supply a more thorough characterization of peatland fungal diversity. The chronological characterization will further facilitate proposal of potential lineages between fungal communities and soil age as well as C accumulation rate.

Study area and sampling description

Peat cores were sampled from three different minerotrophic fens, Shenjiadian (S), Honghe (H), Qindelie (Q), in the Sanjiang Plain (129°11ʹ-135°05ʹ E, 43°49ʹ-48°27ʹ N), north-eastern China (Figure 1). The three fens were started to develop during different periods in Holocene [23]. Over 70% of this region is dominated by fresh water wetlands developing in ancient riverbeds and waterlogged depressions [24] and approximately 30% (or nearly 3.3 × 104 ha) of this region is covered by peatlands, which developed under certain topographic conditions during the Holocene or earlier [25].

Figure 1: Map of the study region showing locations of sampling sites in the Sanjiang plain, which was generated by Zhenqing Zhang using ArcGIS 10.0.

The study area has a temperate humid to sub-humid continental monsoon climate. The mean annual temperature ranges from 1.4 to 4.3°C, with an average maximum of 22°C in July and an average minimum of -18°C in January. The mean annual precipitation is 500-650 mm, and 80% of rainfall occurs between May and September [26]. The sampling fens receive water inputs from groundwater as well as precipitation, and are primarily covered with sedges (Carex lasiocarpa). In May 2012, triplicate cores were sampled from each fen using a Russia peat core. The cores were subsampled for chronological, microbiological and physico-chemical analyses. The soil sample handling for chronological analysis is described in more detail in the work of Zhang et al. [23]. For microbiological and physico-chemical analyses, surface peat soil (0 to 30 cm) was homogenized in sterile bags. The peat soils were stored in the dark and kept chilled before transferred to laboratory. Chronological and physico-chemical characterization was conducted right after the samples arrived. Samples for DNA extractions were immediately frozen at -80°C.

Chronological and physico-chemical characterization

Subsamples with a volume of 3 cm³ were used for loss-on-ignition (LOI) with sequential combustion at 500°C to estimate organic matter [27] The bulk density with a 1 cm interval of each peat core was calculated according to the dry weight and volume of each subsample. Ash-free (organic matter) bulk density was calculated from the measurements of the bulk density and organic matter contents. All subsamples for AMS dating were dated with an accelerator mass spectrometry system at the Institute of Earth Environment, CAS. The AMS 14 C dates were calibrated converted into calendar ages using the program Calib 7.02 based on the INTCAL 13 calibration dataset [28]. The apparent carbon accumulation rates were calculated using calibrated AMS 14C ages, ash-free bulk density measurements and the C contents of peat organic matter in peatlands (using 52% C in peat organic matter) [29]. Soil pH was measured in a 1:5 soil/water suspension [30]. The total nitrogen (N) in the soil was determined by dichromate oxidation using the Continuous Flow Analytical System (SAN++, SKALAR, Netherland).

DNA extraction, sequencing

Genomic DNA was extracted from 0.5 g of peat soil using a FASTDNATM SPIN Kit for soil (MPBio, Santa Ana, CA, USA) according to the manufacturer’s instructions. The DNA concentration was measured using a NanoDrop 2000 spectrophotometer (NanoDrop Technologies, Inc., Wilmington, DE, USA). The V3-V4 region of the bacterial 16S rRNA gene was amplified using the primer pair ITS3F (GCATCGATGAAGAACGCAGC), ITS4R (TCCTCCGCTTATTGATATGC) combined with Illumina adapter sequences, and barcodes [31]. PCR reactions were performed in a 30 μL mixture containing 3 μL of each primer (2 μM), 10 μL of template DNA (1 ng/μL), 15 uL of Phusion® High-Fidelity PCR Master Mix (BioLabs, Inc., New England, USA) and 2 μL of water. The following thermal program was used for amplification: 95°C for 1 min, followed by 30 cycles of 98°C for 10 s, 50°C for 30 s, 72°C for 30 s and a final extension step at 72°C for 5 min. Each sample was amplified in triplicate, and the PCR products were pooled and purified using the Qiagen Gel Extraction Kit (Qiagen, Hilden, Germany). Metagenomic sequencing libraries were generated using the TruSeq® DNA PCR Free Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions and pooled at an equimolar ratio. The 250 bp paired-end sequencing was performed on an Illumina HiSeq2000 platform at Novogene Bioinformatics Technology Ltd.

Sequence data pre-processing and statistical analysis

Raw sequences were divided into sample libraries via samplespecific barcodes and truncated after cutting off the barcode and the primer sequence. Forward and reverse reads with at least 10-bp overlaps and less than 5% mismatches were merged using FLASH [32]. Qualities filtering on the raw tags were performed according to the QIIME (V1.7.0, http://qiime.org/index.html) quality control process [33] and all sequences shorter than 200 bp in length and a quality score lower than 25 in the raw reads were removed. The remaining sequences were subjected to chimera removal using UCHIME Algorithm. Uparse (Version 7.0.1001, http://drive5.com/uparse/) was used to classify the operational taxonomic units (OTUs) at the 97% similarity level [34]. The taxonomic identity was annotated using a Blast algorithm against sequences within the UniteDatabase (https:// unite.ut.ee/) using QIIME software [35]. To study the phylogenetic relationship of different OTUs and the differences between the dominant species in different samples (groups), multiple sequence alignments were conducted using the MUSCLE software (Version 3.8.31, http://www.drive5.com/muscle/) [34]. Alpha and beta diversity analyses were performed using harmonized data by random subsampling 25,728 sequences. The alpha diversity indices, including Chao1, Shannon, Simpson, ACE and equitability, were calculated. The beta diversity between microbial communities was evaluated using both weighted and unweighted unifrac distances. The linear discriminant analysis (LDA) effect size (LEfSe) algorithm was used for high-dimensional biomarker discovery [36] using the non-parametric factorial Kruskal-Wallis (KW) sum-rank test [37] to detect features with significant differential abundance with respect to the class of interest; biological consistency was subsequently investigated using a set of pairwise tests among subclasses using the (unpaired) Wilcoxon rank-sum test [38]. As a final step, LEfSe used LDA to estimate the effect size of each differentially abundant feature and to perform dimension reduction, when necessary. Pearson correlation analyses were used to correlate the relationships between the soil geochemical and microbial parameters [39]. Differences in soil properties across samples were determined using ANOVA followed by Least Significant Difference (LSD) test in IBM SPSS (version 19.0, Chicago, IL, USA) [40]. The beta-diversity community was compared by permutational MANOVA [41]. The Illumina sequencing data in the present study has been deposited into NCBI SRA database with the accession number as No. SRP082472.

Physic-chemical and chronological characterization of peat

9 peat cores were retrieved from S, in Sanjiang Plain, north-eastern China. The AMS dating results indicated that S H, and Q fens was developed during different periods, and had different C accumulation rates. A summary of soil physic-chemical characteristics was presented in (Table 1). Soil pH was all acidic and varied from 4.80 to 5.44. Soil total C and N ranged from 286.88 to 430.12 g kg^-1 and from 11.62 to 34.05 g kg^-1, respectively.

Table 1: Soil properties including the AMS dating result and carbon accumulation rate of samples from 9 peat cores in Sanjiang plain.

Fungal community diversity

A total of 533,323 sequences targeting the ITS gene were obtained from 9 surface (0-30 cm) soil samples using Illumina HiSeq sequencing, ranging from 53,490 to 63,861 reads per sample. After OTU clustering at 97% sequence identity, a total of 989 OTUs were subsequently generated after resampling with 43,303 sequences per sample. Shannon, Simpon, Chao1, ACE index and equitability were calculated to estimate microbial richness and evenness (Table 2). However, there were no significant differences between these three fens in all diversity indexes (p>0.05).

Table 2: Diversity indexes of fungal community in Shengjiadian (S), Honghe (H) and Qindelie (Q) fens.

Microbial community composition

The dominant fungal phyla across all soil samples were Ascomycota and Basidiomycota, with relative abundances ranging from 79.75% to 54.07% and 10.04% to 25.38%, respectively (Figure 2).

Figure 2: Fungal community structure variation in Shengjiadian (S), Honghe (H) and Qindelie (Q) fens. Relative abundance of bacterial at phylum level was shown. Each bar represents the relative abundance of each sample. Each color represents a particular phylum. The numbers in the sample names indicate the sampling depth. The mean of 3 samples taking from the same site is shown.

The relative abundances of minor phyla Chytridiomycota, Zygomycota, Glomeromycota and Rozellomycota were all lower than 1%. In addition, numerous sequences could not be classified to known fungi with relative abundances varying from 1.65% to 19%. The 10 most abundant fungal OTUs were affiliated with two different phyla (Ascomycota and Basidiomycota) and four different orders (Incertae sedis , Agaricales, Helotiales, Tremellales) (Figure 3).

Figure 3: Phylogenetic tree of the top 10 OTUs revealed in the fungal sequences. Each color represents a particular OTU. The second and outermost layer indicated the relative abundance and confidence, respectively.

A Venn diagram was used to compare the similarities and differences between the communities in the three fens (Figure 4a). The Shenjiadian, Honghe, and Qindelie fungal communities had 134 OTUs in common, and 149, 228, 240 unique OTUs, respectively. The unique OTUs accounted for 34%, 45% and 43% of the total detected OTUs in Shenjiadian, Honghe and Qindelie fens.

The fungal community structure was distinct in the three sites (r=0.375, p=0.035) and dominated with sequences belonging to Ascomycota (Figure 1). The highest relative abundance (80%) of Ascomycota sequences was observed in the youngest fen, Shenjiadian. In contrast, less than 54% of Zygomycota sequences observed in the most ancient fen, Qindelie. The three sites showed a similar relative abundance of Zygomycota sequences, while the Glomeromycota and Chytridiomycota sequences in Honghe were approximately five-fold and two-fold in relative abundance compared with the other two sites, respectively. The biomarkers explored using the LEfSe analysis of Shenjiadian, Honghe and Qindelie fens were affiliated with Tremellomycetes, Eurotlomycetes and Helotiaceae, respectively (Figure 4b).

Figure 4: LEfSe anaylsis (a and b) and Venn diagram (c) shown the unique and shared OTUs of fungal community between Shengjiadian (S), Honghe (H) and Qindelie (Q) fens. In the cladogram (a), the circles from the inside out radiation represent bacterial taxon from kingdom to family, and the diameter of the circle is in proportion to the relative abundance of this taxon. The taxon with significant difference is marked with same color as the sampling site in which the taxon appears the highest-ranked variety, and the branch area is shaded correspondingly. The taxon without significant difference is marked yellow. In the histogram (b) of LDA score, the bacterial groups with statistically significant among sites and depths were shown (P<0.05, LDA score>4.0), and the column length represents the LDA score.

Fungal distribution link to the soil properties

Pearson correlation analyses were used to correlate the relationships between soil properties and relative abundance of the most abundant (top 10) fungal taxa at different levels (Figure 5). At phylum level, only the relative abundance of Basidiomycota was significantly (r=0.708, p=0.033) correlated with soil C/N ratio positively. Another dominant phylum exhibited no correlation with measured soil properties. The correlation analysis of the dominant fungal classes and soil properties revealed that the abundance of both Lecanoromycetes and Eurotiomycetes increased with increasing TN and C accumulation rate. At order level, we found that the relative abundance of Sebacinales was positively correlated to TN (r=0.890, p=0.001) and C accumulation rate (r=0.847, p=0.004). Pleosporales was negatively correlated to soil bulk density (r=-0.727, p=0.026), and in the meantime positively correlated to soil age (r=0.873, p=0.002). Filobasidiales also exhibited a positive correlation with soil age (r=0.681, p=0.043). No significant correlation between soil properties and relative abundance of fungal taxa was observed at genus level. And, among all the measured soil properties, soil pH was the only one which was not found to be correlated to fungal abundance at any level (data not shown).

Figure 5: The relationship between the relative abundances of abundant fungal phyla (A), fungal classes (B), fungal orders (C), fungal genera (D) and soil properties.

The fungal taxa observed in the present study primarily belonged to 2 phyla: Ascomycota and Basidiomycota, consistent with previous studies of peat soils [42-44]. Previous studies reported that Ascomycota and Basidiomycota were capable of aerobic degrading dissolved organic matter (DOM) including cellulose and polyphenolic compounds [45]. The predominance of Ascomycota and Basidiomycota at the surface peat soil was consistent with this ability. The relative abundance of Zygomycota of different sites was similar, suggesting that their physiology is distinct from that of the Ascomycota and Basidiomycota. Zygomycota can survival over long periods of dormancy by producing thick-walled, resistant spores [46]. Additionally, members of Zygomycota were not capable of using cellulose and sucrose degradation products, but instead could use carbon substrates of animal and fungal origin, such as fungal hyphae [47]. For example, the important member of Zygomycete, Mortierella spp., was observed in all samples. They can degrade chitin, the essential component of fungal hyphae, as efficiently as chitinolytic actinomycetes [48,49]. We believe that Zygomycota may play an important role in the peatlands carbon cycle as well as Ascomycota and Basidiomycota.

In this study, fungal communities appeared to respond differently to soil physic-chemical and chronological characterizations. Statistical analysis indicated that pH is less important as an environmental force in shaping fungal community structure. We thought this finding could be attributed to the ability of fungi to tolerate a wider pH range for optimal growth as well as their optimal extracellular enzyme activity at low pH [50]. Furthermore, our results suggested that except for TC, TN, C/N ratio and bulk density, soil age and the C accumulation rate were also important in structuring fungal distribution. Several studies demonstrated that fungal community distribution pattern was mainly effected by peatland vegetation, DOC, DON [15,51,52]. To our knowledge, this was the first study to report the feature that fungal community were significantly correlated to soil age and C accumulation in peatland ecosystems. A mechanistic understanding of the role of important fungal taxa in peatland carbon cycling required additional field experiments and ecophysiological studies in the laboratory.

In conclusion, the relative abundance, distribution, and composition of fungal communities in three different minerotrophic fens distributed in the Sanjiang Plain, the southern edge of northern peatlands, was investigated by next-generation sequencing. We captured a rich fungal community and confirmed that the dominated taxa were also frequently detected in other northern peatland ecosystems. TC, TN, C/N ratio and bulk density are determined to be important environmental parameters shaping fungal community structure, however, pH was not. Additionally, for the first time, we found the distribution patterns of several abundant fungal taxa were closely related to the soil age and C accumulation rate. However, because of the current limited sampling sites in the experiment, we had to say it was a primary report which could be useful as a reference for researchers in this field. A further more detail detection based on large amounts of sampling sites would be necessary in the future.

This work was financially supported by the Science Foundation of Chinese Academy of Sciences (XDB15030103), the National Basic Foundation (2016YFC0501202), the National Natural Science Foundation of China (41401332, 41671094, 41571255) and the Natural Science Foundation of Jilin Province (20140101017JC).