Higher fungal diversity is correlated with lower CO2 emissions from dead wood in a natural forest.

Wood decomposition releases almost as much CO2 to the atmosphere as does fossil-fuel combustion, so the factors regulating wood decomposition can affect global carbon cycling. We used metabarcoding to estimate the fungal species diversities of naturally colonized decomposing wood in subtropical China and, for the first time, compared them to concurrent measures of CO2 emissions. Wood hosting more diverse fungal communities emitted less CO2, with Shannon diversity explaining 26 to 44% of emissions variation. Community analysis supports a 'pure diversity' effect of fungi on decomposition rates and thus suggests that interference competition is an underlying mechanism. Our findings extend the results of published experiments using low-diversity, laboratory-inoculated wood to a high-diversity, natural system. We hypothesize that high levels of saprotrophic fungal biodiversity could be providing globally important ecosystem services by maintaining dead-wood habitats and by slowing the atmospheric contribution of CO2 from the world's stock of decomposing wood. However, large-scale surveys and controlled experimental tests in natural settings will be needed to test this hypothesis.

Global decomposition of wood releases CO2 (6 to 9.5 Pg C/year123) at similar rates to fossil-fuel combustion (9.5 Pg C/year in 20114). Decomposing wood also serves as essential habitat56. The factors controlling wood decomposition rates are therefore of broad importance to conservation and to carbon cycle-climate feedbacks.

However, temperature and moisture variables only explain minority portions of total variance in decomposition rates78. For instance, Bradford et al.9 reported that regional temperatures explain only 28% of local variance in mass loss.

The diversity of wood-decomposing fungi might explain much of the remaining unexplained variance. In laboratory-inoculation experiments using small numbers of culturable fungal species, wood pieces with higher final fungal diversity exhibited reduced decay rates101112. Inoculated wood placed in the field also showed a negative effect of final fungal species diversity on decay (R2 = 0.1513).

However, in contrast to laboratory experiments, natural wood decomposition involves much higher species diversity, more complex assembly histories, and selective faunal feeding on decomposers141516. Thus, it is important to examine the relationship between fungal diversity and decomposition rates in wood that is colonized and decomposing under natural conditions.

Natural fungal communities can be characterized using metabarcoding17, in which nuclear ribosomal internal transcribed spacer (ITS) regions are PCR-amplified and read using high-throughput sequencing18192021. ITS1 and ITS2 are each sufficiently variable to differentiate fungal species1819 and return similar estimates of OTU (Operational Taxonomic Units) richness and community structure1922.

Here we metabarcoded ITS2 to examine fungal communities in naturally colonized wood pieces sampled across a wide range of decay classes in the Ailao Mountain forest of Yunnan, China. These wood pieces were sampled from a larger experiment involving three tree species (LC: Lithocarpus chintungensis [Fagaceae], LX: L. xylocarpus, and SN: Schima noronhae [Theaceae]) from which naturally occurring dead-wood pieces were regularly measured for CO2 emission rates over three years8. We measured the extent to which variation in the species diversity and composition of fungal communities can explain variation in emission rates.

Results

Taxonomy results

Numbers of fungal OTUs ranged from 17 to 199 across wood pieces, tree species, and sampling dates, with means of 73.8 (LC, Sep 2012), 76.7 (LC, June 2013), 87.0 (LX, Sep 2012), 90.5 (LX, June 2013), 83.7 (SN, Sep 2012), and 88.4 (SN, June 2013).

41.1% of the 1,807 OTUs produced by uclust and 76.3% of the 1,565 OTUs produced by CROP were assigned to Fungi, and the proportions assigned to each fungal class were similar across assignment methods (Table 1). Because we removed non-Fungi reads from the dataset before taxonomic assignment, we attribute the taxonomically unassigned OTUs to the still highly incomplete UNITE and Genbank databases used for taxonomic assignment.

Table 1

Taxonomic assignments to Class level for the ITS2 Operational Taxonomic Units (OTUs).

CROP + Genbank (1,565 OTUs)
Kingdom Fungi	Ascomycota	Basidiomycota	Chytridiomycota	Zygomycota	Glomeromycota	Mortierellales	Phylumunidentified
1194 OTUs	632	319	1	3	1	3	235
76.3% of total	52.9% of Fungi	26.70%	0.10%	0.30%	0.10%	0.30%	19.70%
uclust + UNITE (1,806 OTUs)
Kingdom Fungi	Ascomycota	Basidiomycota	Chytridiomycota	Zygomycota	Phylum unidentified
742 OTUs	378	169	2	5	188
41.1% of total	50.9% of Fungi	22.80%	0.30%	0.70%	25.30%

Uclust and CROP refer to the two OTU-clustering methods used, and GenBank and UNITE refer to the fungal reference databases used (see Methods: Bioinformatic Analyses for details). The taxonomic assignments of the 200 largest OTUs by read number are listed in Supplementary Information S6.

Fungal diversity and CO2 emissions

In June 2013, the month with the highest CO2 emissions, emissions declined with fungal diversity in all wood species (R2 = 26% to 44%, Fig. 1C,F).

Figure 1

Linear regressions of CO2 emission rates on Shannon fungal diversities measured from individually metabarcoded wood pieces.

Top The solid black curve indicates the air temperature. Carets indicate times of CO2 measurements. Blue shading indicates the warm months when wood decomposition is >50% of maximum. (A–C) CO2 emissions decline with increased fungal species diversity in two of the species in June 2012 (LC and SN) and in all three species in June 2013. In September 2012, CO2 emissions are lower, and there is no relationship. The OTU-picking method is de novo clustering with CROP. (D–F). Same as (A–C) but the OTU-picking method is QIIME’s reference-based matching against the UNITE database, with de novo clustering of non-matched reads with uclust. Non-significant regressions are indicated by dashed lines. Shown here are the non-rarefied datasets. Rarefaction does not change the results (Supporting Information S1). LC = Lithocarpus chintungensis, LX = L. xylocarpus, SN = Schima noronhae.

In June 2012, CO2 emissions from LC and SN also declined with fungal diversity, even though we used a fungal diversity estimate taken three months later (September 2012) and even after conservatively omitting an influential datum from SN (high CO2, low diversity) (Fig. 1A,D). The third species (LX) did not return a significant regression, but its CO2-diversity relationship was visually nearly indistinguishable from its congener LC, suggesting that wood species partly governs the emissions-diversity relationship. Variances explained (26% to 28%, Fig. 1A,D) were lower than in June 2013.

Finally, in September 2012, CO2 emissions did not decline with higher fungal diversity (Fig. 1B,E), which is consistent with the generally lower CO2 emissions in September (Fig. 1).

The above results were robust to two OTU-picking methods (CROP and uclust, Fig. 1), rarefaction (non-rarefied shown in Fig. 1; rarefied in Supporting Information S1), and two diversity estimates (Shannon in Fig. 1, Simpson in S1). Regressions using Simpson diversities were generally statistically more significant (S1). We also analyzed after omitting single-read OTUs (which are more likely to be pipeline artefacts23) and achieved the same results, except that the previously non-significant SN regressions in Sept 2012 (Fig. 1A,D) became statistically significant (authors’ unpublished results). In short, the analyses presented in Fig. 1 are conservative.

Chemistry of decomposing wood

We analyzed the chemistry of a separate subset of 27 wood pieces from the larger experiment. Mean CO2 emissions from 2010–2012 showed no correlation with the densities of carbon, nitrogen, phosphorus, or lignin within any of the decay classes, with one exception, nitrogen density in decay class 1 (statistical details in Supporting Information S2).

Pure-diversity versus species-selection effects

Two general mechanisms could explain the observed diversity-function relationships. The first is a ‘pure diversity’ effect where species identity does not matter, only that increased species richness and evenness per se is somehow responsible for slower wood decomposition. The second is a ‘species-selection’ effect where more diverse fungal communities might be more likely to contain particular species that cause slow decomposition and somehow also govern the overall decomposition rate of the wood piece. To differentiate these two, we used a method devised by Sandau et al.24 to generate a parameter λ for each regression in Fig. 1 (statistical details in Supporting Information S3). λ ranges between 0 and 1, with 0 indicating that variation in species composition does not account for variation in CO2 emissions (i.e. a ‘pure-diversity’ effect). For two tree species, LC and LX, λ always took values near zero (Table 2). For the third tree species SN, λ was also nearly zero in June 2012 but took intermediate values in June 2013, suggesting that fungal composition in this tree species at this time had some explanatory power. The general failure to detect composition effects can be observed in the community ordinations (Fig. 2) by noting that the SN/June 2013 samples were the only ones to line up along the CO2 emissions gradient (except the two lowest diversity samples). Not surprisingly, conventional community-analysis tests returned the same conclusion: variation in community composition is not explained by CO2 emissions (statistical details in Supporting Information S4).

Table 2

Estimates of the contribution of fungal community composition to CO2 emissions, using the method of Sandau et al.24.

CO2 sample date	Species	OTU clustering method	λ
Jun-12	LC	CROP	0.000
		uclust	0.000
	LX	CROP	0.002
		uclust	0.000
	SN	CROP	0.001
		uclust	0.004
Jun-13	LC	CROP	0.000
		uclust	0.000
	LX	CROP	0.000
		uclust	0.001
	SN	CROP	0.366
		uclust	0.500

The generated parameter λ varies between 0 and 1, with 0 indicating that variation in composition does not explain variation in emissions. Composition only contributes to explaining variation in one tree species in one sampling date (SN, June 2013). Conventional community analysis (Supporting Information S4) also detected a contribution of composition in SN in June 2013. LC = Lithocarpus chintungensis, LX = L. xylocarpus, SN = Schima noronhae.

Figure 2

Correspondence analysis ordinations of fungal communities, by tree species and sampling date.

Point size is scaled to CO2 emissions, and the gradient represents fungal Shannon diversity. In all ordinations (A–F), CO2 emissions decrease with higher fungal diversity (point size decreases up the gradient, echoing Fig. 1). Also evident is that the lower diversity wood pieces are compositionally very dissimilar to each other and to the higher diversity wood pieces. Left-hand column (A,C,E). June 2012 CO2 vs. September 2012 fungal diversity. Right-hand column (B,D,F). June 2013 CO2 vs. June 2013 fungal diversity. (A,B) Lithocarpus chintungensis. Note that the label for point 14 at the top of A is obscured by the small point size. (C,D) L. xylocarpus. (E,F) Schima noronhae. Shown here are the non-rarefied datasets clustered using CROP (see Methods). Rarefaction or using uclust-clustering does not change the results (Supporting Information S4).

Discussion

We found that naturally colonised wood with more diverse fungal communities decomposes more slowly (Fig. 1), resulting in a negative relationship between fungal biodiversity and the ecosystem function of decomposition. This result suggests positive relationships between fungal biodiversity and the ecosystem services of carbon storage and the provision of decomposing-wood habitat in forests.

Our results are consistent with five published experiments using laboratory-inoculated wood, which have all found negative relationships between fungal diversity and decomposition rates1011121325. The one exception, Valentin et al.26, found a positive relationship, but in that study, field-collected microbial communities were serially diluted and re-inoculated into laboratory wood incubations. Serial dilution does not necessarily remove microbial species, but it does make all species less abundant27, which might have reduced decomposition rates.

In contrast, field studies to date have reported only ambiguous relationships between fungal diversity and decomposition rates. For instance, Hoppe et al.28 found non-significant correlations between fungal OTU richness and decay class (reflecting different numbers of years of decomposition), albeit negative relationships consistent with our results. Van Der Wal et al.29 measured tree-stump decomposition and reported weakly positive effects of fungal species richness (but not Simpson diversity) on sapwood decomposition, but only in late decay. Kurbartová et al.30 found no relationship between wood loss and fungal OTU diversity after 12 years of decomposition but also reported that the least-decayed logs had the highest community diversities, again consistent with our results. All three studies found differences in community composition for logs that differed in remaining undecomposed weights.

Importantly, none of those three field studies made concurrent measurements of fungal diversities and CO2 emission rates, as we did here (Fig. 1). We observe that relationships between CO2 emissions rate and fungal diversity varied from month to month and across tree species (Fig. 1, Table 2, S1), suggesting that fungal activity and composition are dynamic and environmentally responsive. Thus, the fungal community measured after years of decomposition might not reflect the communities that were active during decomposition, obscuring any relationship between mass loss and fungal diversity.

Our study helps to reconcile the differing results found in the published laboratory and field studies, by making concurrent measurements of emissions and fungal diversity in a field setting that is the most natural on the spectrum of possibilities: colonization of locally dominant tree species by the local fungal community, with uninterrupted and full exposure to local environmental variability and the local faunal community, including fungivores, and long-term succession of fungal and other microbial communities. Our collected wood pieces span a range of at least one to fifteen years of decomposition on the forest floor8. Our results suggest that laboratory experiments correctly reveal negative relationships between CO2 emissions and fungal diversity.

Mechanisms

In contrast to wood, microbial diversity is reported to accelerate the decomposition of soil organic matter3132, and this is thought to represent a general pattern3334 (but see Creed et al.35 for leaf litter). We hypothesise that because soil organic matter presents a much higher diversity of resources than does dead wood, niche complementarity amongst decomposer species drives positive relationships between diversity and soil organic matter decomposition.

By contrast, niche overlap provides a plausible biological mechanism for why wood decomposition should slow with fungal diversity. Interference competition has long been predicted to evolve when niche overlap is high and the disputed resource is valuable36. In forests, decomposing wood resources are available to many fungal species, and aggressive interactions are indeed observed among these fungi37. Elsewhere, it has been shown that interference competition reduces virulence (= host consumption rate) in endosymbioses383940 and productivity in bacterial communities4142. Consistent with those findings in other contexts, interference competition can also explain why fungal biomass has been found to explain variance in wood mass loss9. When a piece of wood is colonized by many fungal species, the hypothesized higher levels of interference competition would result in less wood converted into fungal biomass (or CO2). High niche overlap is also consistent with the observed ‘pure diversity’ effect of fungal diversity on emissions (Table 2, S3, S4) since any species should fight all others (note, however, that we could not test for community composition effects at higher taxonomic levels, see Maherali & Klironomos43). Thus, theory suggests that the arrow of causation can be drawn in the direction of fungal diversity driving decomposition rate. Finally, interference competition results in competitive exclusion, which will cause community composition to change over time. This means that measures of fungal diversity made after years of decomposition are unlikely to explain final variance in decomposition.

Sources of error and proposed future experiments

One source of error is that metabarcoding provides only approximate estimates of species frequencies, due to the many errors known to be introduced by metabarcoding, especially PCR-primer mismatches that lead to biased amplification2223. Also, DNA is environmentally persistent, so fungal species that are no longer represented by living colonies might still be detected by PCR. Nonetheless, we found that Shannon and Simpson indices, which both incorporate species frequency information by discounting rare species (here, rare = low-read OTUs), were able to explain variation in CO2 emissions. There are two likely and non-exclusive explanations. (1) Low-read OTUs were more likely to have been the remains of dead species and/or sequence artifacts from the metabarcoding pipeline and thus indeed should be discounted, and (2) in a Norway spruce forest, Ovaskainen et al.21 found that abundances of fungal fruiting bodies and OTU read numbers were positively correlated, suggesting that low-read OTUs indeed represent low-biomass species, which should have weaker influence on decomposition rates.

Another possible source of error is that we did not experimentally control for the age of the wood pieces, and thus an alternative explanation is that the observed correlations between fungal diversity and CO2 emissions rates (Fig. 1) might be caused by sampling along a successional gradient in which older wood pieces have less remaining wood to decompose (and thus lower emissions) and have also accumulated more fungal species. However, we found no relationship between decay class and emissions rates (Methods: Experimental setup and Statistical analyses) in our dataset, nor did we in the 320-piece superset from which our samples were drawn8, whereas this alternative explanation predicts that the least-decayed wood pieces should show the highest emissions. Also, we found mostly ‘pure-diversity’ effects of fungal communities on emissions (Table 2, S3, S4), whereas this alternative explanation invokes a successional sere and thus predicts compositional effects. We suggest a long-term experiment in which even-aged and sterilized wood pieces are allowed to be colonized and sampled for CO2 emissions rates and fungal diversity over many years in the field. We suggest that bacterial communities also be measured for correlations with CO2 emissions, although we caution that to estimate alpha diversity, sequencing effort will need to be much higher than for fungal communities. In addition to DNA-based metabarcoding, it might also be informative to sequence reverse-transcribed RNA from wood samples, in order to isolate the effect of living fungal species.

Conclusions

The slopes of our diversity-emissions relationships (Fig. 1) are steep enough to suggest that even modest declines in fungal diversity in dead wood could cause several-fold increases in CO2 emissions rates. For example, in June 2013, CO2 emissions varied by 5.6- and 14.4-fold among LC and LX wood pieces. These negative relationships between diversity and wood decomposition provide a strong justification to conduct large-scale surveys of the status of fungal biodiversity and its trajectories in the world’s forests. Global forest fragmentation, reduction of tree-species diversity by fires, logging and replanting, the removal of dead trees, and even increased rainfall could all reduce fungal biodiversity in forests644. All these changes could lead to faster wood decomposition. On the other hand, fragmented forests are drier, fungal distributions are being globalized45, and their local diversity is increased by rising CO246. The net effect of all these changes on carbon emissions from the world’s stock of decomposing wood is difficult to predict.

Methods

Site description

This study was conducted in the Ailao Mountains National Nature Reserve, Yunnan, China, which preserves the largest area of undisturbed, subtropical moist forest in China and has a substantial pool of woody debris (branches and logs, 74.9 × 103 kg ha−1, ref. 47). The study site was at an elevation of 2476 m, about 2 km north of the Ailao Field Station for Forest Ecosystem Studies (24.533 °N, 101.017 °E), and receives 1840 mm annual average precipitation. The climate is monsoonal with distinct cool/dry (November to April) and warm/wet (May to October) seasons48. Annual mean air temperature is 11.3 °C with monthly means ranging from 5 to 16 °C. Surface soils (0 to 10 cm) of the area are Alfisols with pH of 4.2 (in water). The surficial organic layer is 3 to 7 cm deep49. The study site is a broad-leaved evergreen subtropical forest, with the canopy dominated by Lithocarpus chintungensis, Rhododendron leptothrium, Vaccinium ducluoxii, Lithocarpus xylocarpus, Castanopsis wattii, Schima noronhae, Hartia sinensis, and Manglietia insignsis50.

Experimental setup

At our site, most woody debris comes from Lithocarpus chintungensis, (LC), Lithocarpus xylocarpus (LX), and Schima noronhae (SN), so we only examined those three species. In early 2010, branches from these three species, already decomposing on the forest floor, were identified to species by a botanist from the Xishuangbanna Tropical Botanical Garden, collected, and cut into a total of 320 wood pieces, sized to fit a field-respiration chamber (ca. 10 cm diameter and 20 to 30 cm length), tagged, weighed, and measured for size and decay class (further details in Liu et al.8). The three decay classes were DKC1 = a knife could not penetrate, DKC2 = a knife could slightly penetrate with appreciable resistance, DKC3 = a knife could deeply penetrate with little resistance51. We used similar-sized pieces to control for potential effects of wood size on fungal communities52. The pieces were placed on the forest floor within a 60 × 3 m belt transect following an elevation contour. We collected source wood from within 500 m of this transect, utilizing about 5% of downed woody debris from these species in these decay classes, potentially arising from about 6000 source trees (D.A. Schaefer, unpublished data).

Each wood piece was initially weighed with a GLL portable electronic balance (accuracy 0.5 g) and for moisture content with an Extech MO210 moisture meter (calibrated as in Liu et al.8). Their volumes were calculated as cylinders, based on length and the average of 5 circumferential measurements along their lengths. From those, initial weight, volume, and density were all calculated. Oven drying of these wood pieces was not done, because it would have altered microbial communities and wood chemistry.

CO2 emissions-rate measurements

Individual wood-piece CO2 release rates were measured in the field in a closed, ventilated chamber (10 L) connected to an infrared gas analyzer (Licor 820, Lincoln, NE, USA). After chamber closure and initial stabilization, linear CO2 concentration increase rates were logged for at least 5 min. Pieces remained in the field for CO2 measurements (within 5 m) and were handled carefully to limit fragmentation. Temperature and moisture were measured for each sample at each sampling time. The wood-piece CO2 release rate (RWD, μmol C g−1 h−1) was calculated as follows:

where ΔCO2 represents the measured CO2 concentration increase (ppm day−1), P is the internal pressure (kPa; measured by the Licor 820), V is the volume of the system (10.08 L, including the chamber volume and tubing volume and Licor optical path), Vs is the volume of the wood piece (L), R is the gas constant (8.314 L−1 kPa−1 K−1 mol−1), Ts is the wood temperature (°C), and WC is the carbon weight of each piece (g; 47% of its dry weight).

These measurements were made eight times from September 2010 to June 2013, approximately every four months (Fig. 1).

Sampling for genetic analysis

The larger ongoing study includes three species, each having three decay classes, thus producing nine strata. 65 wood pieces were selected from all nine strata. In addition, for each stratum, we calculated a mean CO2 emission rate during 2010–2012, and the pieces chosen for this experiment included ones that were consistently below, at, or above the mean rate, and also some that showed variation around the mean over time. Pieces for each of these subgroups were selected at random from the larger study. This stratified random sampling ensured that all nine groups and the full range of CO2 emission variation were represented in the metabarcoded samples. Those wood pieces were collected in September 2012, and 59 of those samples were recollected in June 2013 (6 samples were not relocated in June 2013, having lost their tags or been buried under new litter). We chose June and September because CO2 emissions from wood are typically higher in June and lower in September, despite similarly warm temperatures and high moisture availability8 (Fig. 1 top). Superficial litter and bark were removed from each wood piece with a flame-sterilized knife before drilling. An electric drill with a flame-sterilized, 11-mm drill bit was used to extract wood powder at three holes located near the ends and the middle of each wood piece. Each sample consisted of ~5 cm3 pooled material from those three holes, collected onto aluminum-foil sheets, and then immediately stored in 50 ml tubes and frozen at −20 °C for 2–3 days until transport on ice packs to our laboratory 10 hr away, where they were stored at −40 °C until extractions.

Comparing cumulative CO2 emissions and gravimetric weight losses

Sixty-four additional wood pieces (i.e. not used for metabarcoding) were retrieved from the field in April 2013 to test the extent to which CO2-emissions-estimated mass loss (averaged over the measurements taken from the year 2010 experiment start) accurately estimated directly measured mass loss (gravimetric weight loss). In the laboratory, these wood pieces were re-measured for volume (as above), and twenty-two wood pieces exhibiting >15% volumetric weight loss since the start of the experiment, indicating substantial fragmentation in the field, were excluded from the analysis. The remaining wood pieces were dried at 70 °C to constant weights and then individually weighed to the nearest 0.1 g on an electronic balance. CO2-estimated mass loss was positively correlated with gravimetric loss (linear regression, Gravimetric loss in grams over 3 years = 30.7 + 0.741 * CO2-based-decomposition Cd (R2 = 0.769, n = 17, p < 0.01); 28.7 + 0.56 * Cd (R2 = 0.762, n = 13, p < 0.01); and 40.3 + 0.468 * Cd (R2 = 0.592, n = 13, p < 0.01), for DKC1, 2, and 3, respectively. Inspection of scatterplots (Supporting Information S5) revealed that the main discrepancy was that CO2-emission mass-loss estimates slightly underestimated small mass losses.

DNA extraction, PCR, and 454 pyrosequencing of ITS2 amplicons

Total DNA was extracted from each sample of wood powder by adding 10 mL CTAB buffer (2% cetyl trimethyl ammonium bromide, 50 mM NaCl, 5 mM EDTA, 10 mM Tris, pH 8) and 20 μL β-mercaptoethanol per 5 cm3 of sample, homogenizing with a TH-02 homogenizer (Omni International, Kennesaw, GA USA) for 5 min at room temperature, incubating at 65 °C for 1 hour, centrifuging at 4,000 rpm for 1 min. After centrifugation, the supernatant was transferred to new Axygen® 2.0 mL microcentrifuge tubes and extracted using one volume of chloroform by vortexing for 20 min and centrifuged at 4 °C, 12,000 rpm for 10 min. The supernatant was then transferred to new microcentrifuge tubes and precipitated with 1.5 volumes of precooled isopropanol at −20 °C overnight. After centrifugation at 4 °C, 12,000 rpm for 20 min, the precipitate was washed with 70% ethanol and dissolved in 100 μL TE buffer (10 mM Tris, 1 mM EDTA, pH 8.0). DNA was purified by using QIAquick PCR Purification Kit. The quantity and quality of purified DNA was assessed with a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). Samples were PCR amplified using the forward primer ITS3_KYO2 5′-GATGAAGAACGYAGYRAA-3′ 18 and reverse primer ITS4 5′-TCCTCCGCTTATTGATATGC-3′ 53. The standard Roche A-adaptor and a unique 10 bp MID (Multiplex Identifier) tag for each sample were attached to the forward primer. PCRs were performed using approximately 10 ng DNA in a 20 μL reaction mixture containing 2 μL of 10X buffer (Mg2+ Plus), 0.02 mM dNTPs, 20 μg Bovine Serum Albumin, 1 μL DMSO, 0.4 μM of each primer, and 0.5 U HotStart Taq DNA polymerase (TaKaRa Biotechnology Co., Dalian, China) under a temperature profile of 95 °C for 10 min, followed by 35 cycles of 94 °C for 20 sec, 47 °C for 30 sec, and 72 °C for 2 min, and final extension at 72 °C for 7 min. For pyrosequencing, PCR products were gel-purified using Qiagen QIAquick PCR purification kit, quantified using the Quant-iT PicoGreen dsDNA Assay kit (Invitrogen, Grand Island, NY, USA), pooled and A-amplicon-sequenced on a Roche GS FLX (Branford, Connecticut, USA) at the Kunming Institute of Zoology.

Bioinformatic analyses

The sequences obtained were run through a pipeline for quality control, denoising and chimera removal, OTU-picking and taxonomic assignment. Quality Control: Header sequences and low-quality reads were removed from the raw output in the QIIME 1.8.0 environment (split_libraries.py: -l 100 -L 500 -H 30)54. The 65 samples collected in 2012 had been sequenced on four 1/8 regions, producing 679,361 raw reads and 525,679 post-quality-control reads (mean read length 286 bp). The 59 samples collected in 2013 had been sequenced on three 1/8 regions, producing 327,226 raw reads and 256,996 post-quality-control reads (mean read length 308 bp). These amplicon lengths were consistent with the expected mean length for ITS2 (327.2 bp, SD = 40) reported by Toju et al.18. Denoising and chimera removal: Denoiser in QIIME55 was used to remove characteristic 454 sequencing errors. Next, ITSx 1.0.356 was used to extract the variable ITS2 region from the whole reads (i.e. conserved 5.8 S and LSU flanking sequences were stripped) and to remove non-fungal-ITS reads. The extracted sequences were clustered at 99% similarity with USEARCH v7.0.109057 to remove replicate sequences and chimeras. OTU-picking and taxonomic assignment: We used two methods to cluster the reads into OTUs. First, we used a reference-based method in QIIME (pick_open_reference_otus.py: max_accepts 20 max_rejects 500 stepwords 20 length 12 –suppress_align_and_tree) in which reads were first clustered by matching at 97% similarity to the UNITE 12_11 fungal database58, which itself had previously been clustered at 97% similarity for use within QIIME. Unassigned reads (the vast majority) were then clustered de novo using the uclust option at 97% similarity, producing 1,807 OTUs in total. For these latter OTUs, we attempted to assign taxonomies using QIIME’s assign_taxonomy.py against the UNITE database. Second, we performed de novo 97%-similarity clustering with CROP 1.3359, producing 1,565 OTUs. We assigned taxonomies using NNCAuto in Claident60 version 0.1.2014.02.06 with database version 20140205, which is derived originally from GENBANK. Sequence data are deposited at datadryad.org (doi: 10.5061/dryad.k6j0h) and in GENBANK’s Short Read Archive (Accession number: PRJNA252416). An example bioinformatic script is in Supplementary Information and also deposited at datadryad.org (doi: 10.5061/dryad.k6j0h).

Statistical analyses

Analyses were performed using vegan 2.0–1061 and mvabund 3.8.462 in R 3.1.063. The HTML outputs of the R scripts are in Supplementary Information and also deposited at datadryad.org (doi: 10.5061/dryad.k6j0h). From both the OTU tables generated (CROP and uclust, see Bioinformatic analysis), we deleted one SN wood piece that had only 35 reads and then split the tables by wood species (LC, LX, and SN) and sample time (September 2012 and June 2013). We then generated a second pair of OTU tables by rarefying to the lowest read number per wood piece in the dataset (rrarefy() in R). Thus, for each wood species and sampling time, we have four OTU tables: CROP/non-rarefied, CROP/rarefied, uclust/non-rarefied, and uclust/rarefied. Finally, for each table, we used vegan’s diversity() function to estimate Shannon and Simpson diversity in each wood piece.

For each wood species, we used lm() in R to linearly regress emissions against fungal diversity. Thus, the June 2012 and September 2012 CO2 estimates were tested against the September 2012 fungal diversity estimate, and the June 2013 CO2 estimate was tested against the June 2013 fungal diversity estimate. Residuals were all adjudged visually to be near normally distributed, but with small indications of nonlinearity. We ignored the nonlinearities because the residuals suggested accelerating CO2 emissions at the lowest fungal diversity, making our results conservative. In trial models, we also tested for significant effects of wood surface temperature and decay class (see Experimental setup), but they did not interact significantly with the fungal diversity term and mostly did not enter significantly as additive terms, and so have been omitted here for simplicity. Liu et al.8 also did not find a correlation between decay class and CO2 emissions rates for these wood pieces.

To test whether a ‘pure diversity’ effect is sufficient to explain the observed diversity-function relationships, we use a method by Sandau et al.24 where community similarities (1-Jaccard binary) between all wood pieces are used to create a variance-covariance matrix that is then included in the linear regressions, thus taking into account potential non-independence of wood pieces due to the fact that some communities are similar to each other (Supporting Information S3). In Supporting Information S4, we also use conventional community analyses to test for an effect of community composition on CO2 emissions. We limit our tests to June 2012 and 2013, as only these exhibited significant declines in emissions with fungal-species diversity (Fig. 1A,C,D,F).

Additional Information

How to cite this article: Yang, C. et al. Higher fungal diversity is correlated with lower CO2 emissions from dead wood in a natural forest. Sci. Rep. 6, 31066; doi: 10.1038/srep31066 (2016).