Monday, 17 October 2016

Do foliar fungal communities of Norway spruce shift along a tree species diversity gradient in mature European forests?

Published Date
Open Access, Creative Commons license

  • Diem Nguyen ,
  • Johanna Boberg
  • Katarina Ihrmark
  • Elna Stenström
  • Jan Stenlid
    • Department of Forest Mycology and Plant Pathology, Swedish University of Agricultural Sciences, Box 7026, 75007 Uppsala, Sweden
    Diversity of species is considered beneficial for most ecosystems (Cardinale et al., 2012 and Gamfeldt et al., 2015). Mixed forests have been associated with higher levels of ecosystem services (Gamfeldt et al., 2013 and Carnol et al., 2014) and are thought to reduce the risk of fungal pathogen disease susceptibility as compared to monoculture stands (Pautasso et al., 2005 and Felton et al., 2016). High tree species diversity, as proposed by the insurance hypothesis, may maintain the overall integrity of a forest ecosystem by reducing this risk (Yachi and Loreau, 1999). Both foliar pathogens and endophytes may be affected by the tree diversity in the stand (Müller and Hallaksela, 1998Hantsch et al., 2013 and Nguyen et al., 2016), although the mechanism for such an effect is not yet clear.
    Tree leaves and needles host a range of organisms, including fungi. These foliar fungal species may have beneficial, antagonistic or no apparent impact on the tree (Carroll, 1988 and Rodriguez et al., 2009). Fungi can be found on the leaf surface as epiphytes (Legault et al., 1989 and Santamaría and Bayman, 2005) or inside leaves as endophytes, not causing any obvious symptoms (Carroll and Carroll, 1978Petrini, 1992 and Müller and Hallaksela, 2000). The distribution and abundance of foliar fungi vary not only among host species (Deckert and Peterson, 2000) but also among genotypes of one species (Bálint et al., 2013). At the individual plant level, the fungal communities are diverse and can differ within a plant (Müller and Hallaksela, 2000Arnold et al., 2003 and Cordier et al., 2012a) and also within a leaf (Lodge et al., 1996and Arnold et al., 2000). Across large geographical scales, variation in the composition of foliar fungi has been observed in conifers (Terhonen et al., 2011 and Millberg et al., 2015). Such differences have often been attributed to variations in temperature, precipitation patterns, vegetation zones, geographic distance and other environmental factors (Helander, 1995Vacher et al., 2008 and Zimmerman and Vitousek, 2012). Fungal communities of coniferous hosts have higher diversity at higher latitudes (Arnold and Lutzoni, 2007 and Millberg et al., 2015). At smaller spatial scales, differences in the foliar fungal community may be caused by forest structure, host density, microclimate and the surrounding environment (Saikkonen, 2007). Given that foliar fungal communities respond to all of these factors, they might be expected to vary with the surrounding vegetation, although this has not been studied in detail.
    The fungal community can be studied in a number of ways. Molecular-based approaches, in contrast to culture-based methods, allow the detection of many more species (Amann et al., 1995), including species that cannot be obtained in culture. The recent advances in high-throughput sequencing technology have increased both the resolution and scope of fungal community analyses and have revealed a highly diverse and complex mycobiota of plant foliage (O'Brien et al., 2005Jumpponen and Jones, 2009Jumpponen et al., 2010Cordier et al., 2012bMenkis et al., 2015 and Millberg et al., 2015). Most studies have so far described the fungal community by sequencing the ribosomal RNA genes (rDNA), which provide a description of all members of the community, regardless of activity level. For example, dead organisms with intact genetic material, resting spores and vegetative mycelia may be detected by sequencing methods (England et al., 1997 and Demanèche et al., 2001). However, sequencing the ribosomal RNA transcripts (rRNA) instead will reveal the metabolically active and functionally important taxa of the community and provide insights into the activity of these fungi in environmental samples (Pennanen et al., 2004Baldrian et al., 2012 and Delhomme et al., 2015). The metabolically active fungi would presumably reflect the portion of the fungal community that is responding to variation in the environment, such as vegetation gradients. By focusing on these fungi, the influence of external factors shaping foliar fungal communities can be elucidated.
    The overall aim of the study was to determine whether fungal communities associated with Norway spruce (Picea abies) needles correlate with variation in tree species richness along a tree species diversity gradient. Furthermore, we investigated whether the metabolically active community responds differently than the total community to this tree diversity gradient. To that end, current-year needles were collected from Norway spruce trees across a tree species diversity gradient in four mature European forests that represent different mature forest types (Baeten et al., 2013), and the fungal communities were analyzed using 454 pyrosequencing. Additionally, at one forest site, both RNA and DNA were sequenced from the same foliar samples to compare the metabolically active fungal community (RNA) relative to the total fungal community (DNA), respectively. We hypothesized that the fungal communities and specific fungal taxa would correlate with varying mixtures of tree species along a tree species diversity gradient. Further, we expected that the metabolically active fungal community (RNA), and not the total fungal community (DNA), would shift along the tree species diversity gradient, as they are likely to actively respond to changes in the microenvironment created by changes in tree species composition in the mixtures.

    2 Material and methods

    2.1 Sampling sites and collection

    The study was conducted in mature forests in four countries (i.e. Finland, Romania, Germany and Poland), spanning four major European forest types, established within the FunDivEUROPE Exploratory Platform (Baeten et al., 2013) (Fig. 1Table 1). Sampling was conducted over approximately 2 weeks in 2012 or 2013 for each of the four forests during the vegetation period.
    Fig. 1. Map of European sampling sites in different forests utilized to study the foliar fungal community of current-year Norway spruce needles.
    Table 1. Sampling site locations, forest type and composition of forest mixture with respect to Norway spruce.
    Forest siteRâşcaHainichBiałowieżaNorth Karelia
    Forest typeMountainous beechTemperate BeechHemi-borealBoreal
    Species poolPicea abies
    Abies alba
    Fagus sylvatica
    Acer pseudoplatanus
    Picea abies
    Acer pseudoplatanus
    Fagus sylvatica
    Fraxinus excelsior
    Quercus robur
    Picea abies
    Pinus sylvestris
    Betula pendula
    Carpinus betulus
    Quercus robur
    Picea abies
    Pinus sylvestris
    Betula pendula
    Location47.3° N, 26.0° E51.1° N, 10.5° E52.7° N, 23.9° E62.6° N, 29.9° E
    Study area size (km × km)5 × 515 × 1030 × 40150 × 150
    Mean annual temperature, °C6.
    Mean annual precipitation, mm800775627700
    Sampling periodJuly 2013July 2012July–August 2013August 2012
    Number of plots (trees) with Picea abies
    Monoculture2 (12)2 (12)2 (12)4 (24)
    Two species5 (15)3 (8)5 (15)8 (24)
    Three species5 (15)3 (8)6 (18)4 (12)
    Four species3 (9)2 (6)8 (24)
    Five species2 (6)
    Total15 (51)10 (34)23 (75)16 (60)
    Number of trees sampled and sequenced
    51347560/55 (all sites/Finland RNA)
    Standardized plots of 30 × 30 m were delimited within each forest, where different compositions of tree species were targeted to create a tree species diversity gradient with richness levels ranging from monoculture to three- (North Karelia, Finland), four- (Râşca, Romania and Hainich, Germany) or five-species (Białowieża, Poland) mixtures and different tree species assemblages at each level of species richness (Table 1). Focal Norway spruce trees were randomly selected from a pool of those trees with the largest diameter at breast height (19–79 cm). Within each plot, six trees in monoculture plots and three trees in mixtures were sampled. In total, 220 trees in 64 plots were sampled from the four forests (Table 1).
    From each tree, two branches were cut from the southern exposure: one from the sun-exposed upper part of the canopy, and one in the lower third of the canopy. Shoots were collected from these branches from each tree. Per branch, five current-year shoots were sampled resulting in a total of 10 shoots per tree that were placed into one paper bag per tree. Each sample represented a tree. To prevent changes in the fungal community that could result from growth of opportunist organisms, shoots with needles still attached were immediately dried at 60 °C for 3 d (samples from forests in North Karelia and Hainich). Dried needles, detached from their shoots, were mixed in their respective paper bags. When it was not possible to dry the shoots within 24 h, which was the case for samples from forests in Râşca and Białowieża, samples were stored at 4 °C for a maximum of 2 weeks, and at −20 °C until further processing. These samples were then freeze dried for 3 d. Needles were removed from the shoots and mixed in their respective papers bags. For all samples, a subsample of 20 needles was removed randomly from each bag for further analysis of the total fungal community.
    In addition to the total fungal community from these four forests, an assessment of the metabolically active and the total fungal community was studied in current-year needles from the forest in North Karelia, Finland. The same branches and shoots collected as described earlier were used. Needles were collected immediately after cutting down branches from each of the 60 trees. Two needles from each of the five shoots from the top branch and two needles from each of the five shoots from the bottom branch, 20 needles in total, were collected directly into 2 mL screw-cap centrifuge tubes with 1 mL of RNAlater (Thermo Fisher Scientific, Waltham, USA), without any sterilization procedure, and stored at 4 °C for a maximum of 2 weeks, until longer term storage at −20 °C was possible.

    2.2 Molecular detection of the needle-associated fungal community

    Samples originating from Râşca, Hainich, Białowieża and North Karelia that were used to study the total fungal community (all sites study) were prepared separately from those samples from North Karelia that were used for the study of the metabolically active and total fungal community (Finland RNA study). The fungal communities from both studies were determined with high-throughput sequencing of the internal transcribed spacer (ITS) region of the ribosomal RNA genes.
    Prior to DNA extraction for the all sites study, desiccated needles were washed with 0.1% Tween-20 solution to remove ephemerally attached organisms and particles. Subsequently, the needles were dried on clean filter paper and transferred to a 2 mL screw-cap centrifuge tube together with two metal nuts that fit into the tube, and homogenized using a bead beater (Precellys 24, Bertin Technologies, Rockville, USA) at 5500 RPM for 20 s twice, with a 10 s pause in between, until a powder was produced. Genomic DNA was extracted with CTAB buffer (3% cetyltrimethylammonium bromide (CTAB), 2 mM EDTA, 150 mM Tris–HCl, 2.6 M NaCl, pH 8). Chloroform was used to remove protein contaminants. DNA was precipitated with 2-propanol, washed with 70% ethanol, and resuspended in water. Extraction negative controls (i.e. centrifuge tubes that did not include samples) were also included.
    The fungal ITS2 region was amplified with primers gITS7 (Ihrmark et al., 2012) and ITS4 (White et al., 1990), with ITS4 extended with a unique 8-base pair (bp) sample identification barcode for each sample. The resulting amplicons were 250–400 bp in length. Amplification of each sample occurred in 50 μL reactions [0.025 U μL−1DreamTaq DNA Polymerase and buffer (Thermo Fisher Scientific, Waltham, USA) 200 μM dNTP, 500 nM gITS7, 300 nM ITS4, 2.75 mM MgCl2, and 0.125 ng μL−1genomic DNA template (or 1:10 dilution of extraction controls) was performed using an Applied Biosystems 2720 Thermal Cycler (Applied Biosystems, Carlsbad, USA). PCR negative controls (no samples added) were used to evaluate that there was no contamination during the preparation for PCR. The PCR cycle parameters consisted of an initial denaturation at 95 °C for 5 min, 26–35 cycles of denaturation at 95 °C for 30 s, annealing at 56 °C for 30 s and extension at 72 °C for 30 s, followed by a final elongation step at 72 °C for 7 min. The number of amplification cycles was determined individually for each sample to preserve the fungal genotype composition (Lindahl et al., 2013). Thus, PCR was interrupted while in the exponential phase to yield weak to medium-strong amplicons as visualized by gel electrophoresis on 1% agarose gels. The amplicons were purified with the Agencourt AMPure XP kit (Beckman Coulter, Brea, USA) and the concentration was determined using the Qubit Fluorometer and dsDNA High Sensitivity Assay Kit (Thermo Fisher Scientific, Waltham, USA). Each sample was amplified in triplicate and each replicate was individually purified and concentration determined. PCR amplicons were mixed in equal mass proportion into a general sample that was further purified using E.Z.N.A. Cycle-Pure Kit (Omega Bio-tek, Norcross, GA, USA) for samples from Hainich and North Karelia, or JETQUICK DNA Clean-Up Spin Kit (Genomed, Löhne, Germany) for samples from Râşca and Białowieża. There were two libraries, namely one consisting of samples from Hainich and North Karelia, and another with samples from Râşca and Białowieża. Each library was subjected to ligation of sequencing adaptors and separately sequenced with the 454 sequencing technology using GS FLX Junior (Roche, Switzerland) that is equivalent to an eighth of a plate, with Titanium series chemistry, and carried out by Eurofins Genomics (Ebersberg, Germany). DNA extraction controls were sequenced, while PCR negative controls were not sequenced, as nothing was amplified from them.
    Samples for the Finland RNA study, i.e. those from North Karelia, Finland that were used for the analysis of the metabolically active and total fungal community, were subjected to RNA and DNA extraction as follows. Needle samples were removed from RNAlater and ground with nucleic acid-free mortar and pestle, both of which were pre-treated to above 450 °C, with sand and liquid nitrogen to maintain below freezing conditions that ensured RNA stability. Needles were pulverized and transferred to 15 mL plastic tubes and stored at −70 °C until further processing. RNA extraction was performed according to Chang et al. (1993). To each tube, 4 mL of extraction buffer (CTAB extraction buffer supplemented with PVPP [2% CTAB, 25 mM EDTA pH 8.0, 0.1 M Tris-HCl pH 8.0, 2 M NaCl, 2% mercaptoethanol and 2% PVPP]) was added and incubated at 65 °C for 10 min. Two chloroform purification steps to remove proteins were performed with centrifugation at 6800g for 5 min. Overnight precipitation at 4 °C with 2 M lithium chloride and centrifugation at 6800g for 40 min to pellet the RNA were carried out. The RNA pellet was resuspended in 100 μL RNase-free water, transferred to 1.5 mL tubes and precipitated at −20 °C with 2.5 vol ethanol and 1/10 volume of 3 M NaOAc, pH 5.2. The precipitated RNA was centrifuged for 15 min at 16200g, followed by pellet resuspension with 16 μL RNase-free water. DNase I (Sigma, St. Louis, USA) digestion of the RNA removed any remaining DNA contamination. To confirm that there was no DNA remaining in the extracted RNA, amplification with primers gITS7 and ITS4 of the treated material was conducted as above. RNA concentrations were determined fluorometrically with Qubit RNA High Sensitivity Assay Kit (Thermo Fisher Scientific, Waltham, USA). The cDNA synthesis was carried out by using the iScript Select cDNA synthesis kit (Bio-Rad, Hercules, USA) with 5 μL of DNase treated RNA (approximately 1–500 ng of RNA) and 1 μL ITS4 primer, followed by 90 min incubation at 42 °C. PCR with iScript alone resulted in no visible bands on 1% agarose gels.
    The lithium chloride supernatant was retained following the RNA pelleting step (where the DNA fraction was expected to be) and DNA was precipitated with one volume 2-propanol, followed by 70% ethanol wash and resuspended in 200 μL molecular-grade water. Half of the samples had residual lithium chloride, and these were further purified with JETQUICK DNA Clean-Up Spin Kit. Amplification of DNA and cDNA with gITS7 and barcoded ITS4 primers followed the above protocol. DNA concentration of 0.25 ng/μL and for the cDNA, 2.5 μL of the iScript reaction was used for the amplification of ITS2. Amplicon purification, concentration determination, pooling, purification with E.Z.N.A. Cycle-Pure Kit and 454 pyrosequencing occurred as above. There were no extraction controls for the RNA/DNA extraction, though all PCR controls were sequenced. Five samples from one of the three-species mixture plots could not be successfully amplified and thus were not sequenced. In total, there were 55 RNA samples that derived from RNA extraction and a matching 55 DNA samples that were from the same RNA samples.

    2.3 Parsing the sequence data

    The sequence reads (‘reads’) generated from three separate 454 pyrosequencing runs were subjected to quality control and single-linkage clustering in the SCATA bioinformatics pipeline ( All three datasets from each of the sequencing runs (two from the all sites study and one from the Finland RNAstudy) were combined in one SCATA analysis. Quality filtering of the sequences using the high-quality region (HQR) extraction option (HQR is the longest part of a read that fulfills all the quality thresholds) included: the removal of short sequences (<200 bp) and sequences with low mean read quality score (<20) or with a score below 10 at any position. Homopolymers were collapsed to 3 bp before clustering, and restored to their original length before final analyses and downstream sequence identification. Sequences that were missing primer sequences and/or identification barcode sequence were also excluded, allowing for two mismatches in each of the primer sequences. Primer sequences and barcodes were removed, but information on the read associated with the sample was retained as metadata. The sequences were then clustered into different operational taxonomic units (OTUs) using single-linkage clustering based on 1.5% dissimilarity (i.e. 98.5% similarity) using the USEARCH clustering engine. This clustering threshold approximately corresponds to species level, and was validated by including reference sequences. Clemmensen et al. (2015)further elaborates on the SCATA pipeline workflow and provides justification of our procedure here. The most abundant genotype for each cluster was used to represent each OTU. For OTUs containing two sequences, a consensus sequence was produced.
    The construction of a neighbor-joining tree (Lindahl et al., 2013), based on MUSCLE alignment algorithms of sample and reference sequences, aided in the initial removal of non-fungal reads (e.g. plants, algae, one bacterium and a retroelement). OTUs were taxonomically identified by BLAST searches (Altschul et al., 1997) against the GenBank nucleotide database, optimized for somewhat similar sequences (blastn, Uncultured fungal clones or environmental samples were not excluded from the initial BLAST searches. To describe the taxonomic affiliations of the OTUs, two procedures were then used. In the first method, for the most abundant OTUs (as defined for each forest, and with a minimum of 10 reads per sample), species-level identification was attempted. For each of these OTUs, up to 100 matches from different authors were considered, preferably from published studies by authorities in taxonomy. For conflicting matches, the lowest common rank level was used for taxonomic assignment (Peršoh et al., 2010). The sequence must align 100% over the entire length of the queried sequence. The ITS similarities for defining taxa in GenBank were delimited following Ottosson et al. (2015) as follows: species level – 98–100% similarity and e-values below e−100; genus level – 90–97% and e-values below e−90; and order level – 80–89% and e-values e−80. These OTUs were used for further analyses specified below.
    In the second method to describe the taxonomic affiliations of the OTUs, where the interest was to identify the taxonomic rank phylum, class and order for each OTU, several top GenBank matches providing taxonomic rank information were evaluated and phylum, class and order were recorded. Uncultured fungal clones or environmental samples do not provide taxonomic information and were thus excluded. If the top match was Phylum (or Class or Order) sp., for example Ascomycota sp. or Helotiales sp., the match was not selected for taxonomic rank identification. Rather the next best match to species, even if Genus sp. (e.g. Hymenoscyphus sp.) was included for taxonomic rank identification. Where there were multiple matches that had the same high bitscore and similarity, but provided conflicting taxonomic resolution, the lowest common rank level was used for taxonomic assignment. For conflicts at the phylum rank, i.e. one match was to a Basidiomycota species and another to an Ascomycota species, no taxonomic rank was assigned, and ‘Fungus’ was recorded. Phylogenetic classification was assigned after Schoch et al., 2006Hibbett et al., 2007Boehm et al., 2009Hodkinson and Lendemer, 2011Rosling et al., 2011Zhang et al., 2011Boonmee et al., 2014Ertz et al., 2014 and Chen et al., 2015.

    2.4 Statistical analyses of the fungal community diversity and composition

    The study had a hierarchical design. Each forest consisted of many plots that were associated with different tree species richness levels based on tree species composition, as noted earlier. There were multiple plots within a forest with the same tree species richness level (i.e. in North Karelia, there were four plots of Norway spruce in monoculture plots, and in Hainich, there were two plots containing Norway spruce in the three-species mixture plots), and multiple Norway spruce trees (either three or six sampled) within each plot. To avoid pseudoreplication errors and allow for downstream analyses, the reads for each tree sample were pooled at the plot level for each OTU that resulted in plot-level read values.
    Fungal diversity was analyzed using Hill numbers, which has been found to be a reliable metric with microbial sequence data from complex communities (Jost, 2006and Haegeman et al., 2013) and has previously been used for fungal community sequence data (Bálint et al., 2015). Hill numbers can be used to measure and incorporate richness (Hill 0), exponent of Shannon Index (Hill 1) and inverse Simpson (Hill 2) (Legendre and Legendre, 1998). A community is more diverse if all of its Renyi diversities, as represented by the Hill numbers, are higher than in another community. Fungal diversity was assessed at different levels: (i) between forests (that were sequenced and processed in the same way), (ii) across the tree species diversity gradient in each forest, with tree species richness specified as a discrete variable, and (iii) between community types, i.e. active (RNA) and total (DNA) fungal communities. As carried out by Bálint et al. (2015), linear models were used to explain Hill 0, Hill 1 and Hill 2 with the square root of total read numbers for a plot-level sample. The square root of the read numbers accounts for differential sequencing depth in different samples. The partial residuals of the Hill numbers, after accounting for differential sequencing depth, between forests, and among tree species richness levels within a forest or community type were compared with Tukey's HSD.
    Rarefaction curves were constructed separately for each forest, as well as for community type (i.e. active (RNA) and total (DNA) separately), using sample-based data in the software EstimateS version 9.1.0 (Colwell, 2013). In the construction of these curves, each tree sample was individually specified, rather than using the pooled plot-level read values.
    The relative abundance of reads was calculated for each OTU in each plot-level sample and was used for the following analyses. Non-metric multidimensional scaling (NMDS) ordination was used to visualize the multivariate community composition, based on the relative abundance of each OTU. The metaMDS function in the vegan package (Oksanen et al., 2013) was used to make NMDS plots. Bray-Curtis dissimilarity matrices, three dimensions and 100 random starts were specified. The first solution served as the starting point for a second NMDS analysis.
    Additionally, permutational multivariate analysis of variance (PERMANOVA, adonisfunction (Anderson, 2001)) allowed the partitioning of the variance contributed by the explanatory variables, and thus tested significance of the difference between forests, between community types in North Karelia, and among levels of tree species richness for each forest in the all sites study and for each community type in the Finland RNAstudy, respectively. For all PERMANOVA analyses, sample distances were calculated with Bray-Curtis dissimilarity matrices and statistical significance was estimated with 999 permutations.
    For the most abundant OTUs, as defined earlier, the relationship between their relative abundance as a function of tree species richness was investigated. Comparison of the relative read abundance for within-species comparisons is admissible, though between-species comparisons are problematic due to innate sequence structure (Amend et al., 2010). The most abundant OTUs were determined by setting a threshold of 10 reads per sample and calculating the minimum number of reads required. The threshold level for North Karelia was 580 reads, Białowieża 740, Hainich 340, and Râşca was 490.
    Generalized linear mixed-effect models (GLMMs) were used to determine significance of the distribution of individual OTUs that were most abundant along a tree species diversity gradient for the all sites and Finland RNA studies. Specifically, we tested which, if any, of the most abundant OTUs correlated with tree species richness. GLMMs take into account the hierarchical sampling design and account for non-independence among observations and allow for nested and crossed random-effect terms (Zuur et al., 2009 and Schielzeth and Nakagawa, 2013). Consequently, tree-level sequence reads were used for these analyses, and not plot-level reads. The response variable, which was the proportion of a specific OTU, had a binomial distribution. The explanatory variable tested was tree species richness (treated here as a continuous variable to aid in interpretation of the results). The random factors included plot, and to properly account for the variance (that does not vary freely), random residual was estimated by specifying the sampling unit, i.e. a unique identifier for each tree sample. All analyses, unless otherwise specified, were carried out in R version 3.1.3 (R Core Team, 2013). GLMMs were run with the glmer function, with the binomial distribution and logit link (logit = ln(response probability of OTU/probability of the non-response of OTU)) specified, in the lme4 package (Bates et al., 2015). Odds ratios for the effects of explanatory variable can be calculated from the estimated regression coefficients as exponent of the estimated coefficient. In logistic regression, the odds ratios are interpreted as a change in the probability of an OTU to be present when tree species richness increases by one unit (Breslow and Day, 1987Hosmer and Lemeshow, 1989).

    3 Results

    3.1 Fungal community diversity and composition

    Following 454 sequencing of spruce needle samples across all four forests (all sitesstudy) for the total fungal community, 209,267 sequence reads were generated. Of these, 142,578 (68%) passed quality filtering. In the part of the study aimed to compare the metabolically active fungal community to the total fungal community (i.e. community type) from trees in North Karelia, Finland (Finland RNA study), 454 sequencing of current year needles resulted in 108,144 reads of which 86,035 (80%) passed quality filtering. Across both studies, there were 1737 global OTUs (non-singletons) and 1127 global singletons. Non-fungal OTUs (30 OTUs, 29966 reads, 13% of all reads) were removed. The two datasets, namely the all sites study and the Finland RNA study, were analyzed separately. Following the removal of OTUs with fewer than 10 reads globally for each dataset, there were 513 fungal OTUs in the all sites study and 310 fungal OTUs in the Finland RNA study.
    Some OTUs could be putatively identified in both datasets to species level. The abundant species included: Ascochyta skagwayensisAureobasidium pullulans(ubiquitous yeast-like fungus), Celosporium larixicolaCeramothyrium carniolicumChrysomyxa ledi (rust pathogen of Norway spruce), Epicoccum nigrumExobasidium arescensHeterobasidion parviporum (root and butt rot pathogen of Norway spruce), Phaeocryptopus gaeumannii (pathogen of Douglas-fir (Pseudotsuga menziesii)), Ramularia vizellae (plant pathogen), Sydowia polysporaTaphrina carpini (pathogen of oak hornbeam (Carpinus betulus)), and Taphrina vestergrenii (pathogen of fern).
    For the all sites study, a majority of the 513 OTUs were putatively assigned to Ascomycota (394 OTUs, 105,665 sequence reads) while Basidiomycota accounted for 115 OTUs (18,083 reads) (Supplementary Figure 1A). Four OTUs were putatively assigned to Chytridiomycota, Glomeromycota and unidentified fungi, accounting for 298 reads. At the class level, most OTUs were putatively assigned to Dothideomycetes (169 OTUs, 49,736 reads), Eurotiomycetes (60 OTUs, 12,799 reads) and Leotiomycetes (56 OTUs, 21,027 reads). There were 18 OTUs (5011 reads) that could not be assigned to class level. Overall, there were 35 Ascomycota orders, 21 Basidiomycota orders, 9 unknown orders in both phyla, all spanning 18 classes (and one unassigned Ascomycota class).
    For the Finland RNA study, putatively assigned Ascomycota made up the majority of the 310 OTUs and reads (232 OTUs, 47,326 reads), followed by Basidiomycota (74 OTUs, 19,244 reads). In total, there were putatively 29 Ascomycota orders, 15 Basidiomycota orders, and 9 unknown orders in both phyla, in 21 classes (Supplementary Figure 1B). Two unidentified fungal OTUs (399 reads) and two putative Glomeromycota OTUs (28 reads) were also detected.
    Furthermore, in the Finland RNA study, the number of OTUs in the total fungal community (DNA) (303 OTUs, 49,671 reads) and the metabolically active fungal community (RNA) (293 OTUs, 17,326 reads) were similar, although the number of reads was three times higher in the DNA fraction (Table 2). A majority of the OTUs were likewise putatively assigned to Ascomycota (active community: 219 OTUs (12,062 reads) and total community: 228 OTUs (35,264 reads)). There were 70 putative Basidiomycota OTUs (5103 reads) for the active community and 72 OTUs (14,141 reads) for the total community (Table 2). There were 17 OTUs (346 reads) that were uniquely in the DNA fraction, and seven OTUs (110 reads) that were uniquely in the RNA fraction. In general, the unique OTUs were very uncommon and the most abundant of these unique OTUs accounted for less than 0.16% of the total reads. Those OTUs uniquely in the DNA were predominantly Ascomycota (13 OTUs, 297 reads) and the remaining four OTUs (49 reads) were Basidiomycota. Similarly, OTUs uniquely in the RNA were predominantly Ascomycota (four OTUs, 71 reads), two OTUs (26 reads) were Basidiomycota and one OTU (13 reads) putatively assigned to Glomeromycota.
    Table 2. Comparison between community type, i.e. the active (RNA) and total (DNA) fungal community, in the Finland RNA study. Number of OTUs and reads for each putatively assigned orders for the respective community type are presented.
    Number of OTUs in both RNA and DNA
    Number of reads in both RNA and DNA
    Number of OTUsNumber of reads
    Only in RNAOnly in DNAOnly in RNAOnly in DNA
    Number of OTUsNumber of reads
    RNA totalDNA totalRNA totalDNA total
    The distribution of reads for both the all sites and Finland RNA studies showed log-normal distribution, such that many OTUs had few sequence reads and a few OTUs had many reads, that is typical of this type of sequence data (data not shown). Rarefaction curves for the data separated by forest in the all sites study did not reach saturation, suggesting that deeper sequencing would uncover more taxa (Supplementary Figure 2A). On the other hand, rarefaction curves constructed for the data separated by community type in the Finland RNA study showed a nearly saturated curve for the total community (DNA), but not the active community (RNA), suggesting deeper sequencing of the total fungal community than the active fungal community (Supplementary Figure 2B).

    3.2 Fungal community from the all sites study (total community analysis)

    Fungal community composition differed strongly across the tree species richness levels in Hainich (R2 = 0.50, P = 0.023; Fig. 2Supplementary Table 1). However, fungal community diversity, as determined by Hill numbers, was not different along the tree species richness gradient in any of the four forests (data not shown).
    Fig. 2. NMDS results of the fungal community composition of plot-level OTUs for each plot (GERxx, where xx represents the plot number) in the Hainich forest in Germany. Axes are arbitrary, represent two NMDS dimensions and scaled in units of Bray-Curtis dissimilarity. A stable solution was reached (stress value = 0.076). The colored ellipses mark the 95% confidence interval of the group centroids for each tree species richness level. Lines link the samples pooled at the plot level to the group centroids. The blue triangles are plots from the monoculture, black circles enclosed within the orange ellipse represent the two-species mixture, yellow squares enclosed within the red ellipse represent the three-species mixture and the grey squares represent the four-species mixture.
    The most abundant OTUs from each of the four forests varied in their relative read abundance, which was between 1% and 28% at the extremes, and found in varying number of samples (Supplementary Table 2). Negative tree species richness effects were observed for OTU_24 (Dothideomycetes sp.) in Râşca (Table 3). Tree species richness effects were also observed in Hainich; OTU_2 (Helotiales sp.) and OTU_39 (Sporobolomyces sp.) exhibited negative relationships. However positive relationships were seen for OTU_7 (A. pullulans), OTU_6 (Ceramothyrium sp.), OTU_35 (Trichomerium sp.) and OTU_61 (Chaetothyriales sp.) (Table 3).
    Table 3. Generalized linear mixed-effects model results at the fungal taxa level for the most abundant OTUs. Effect of tree species richness on OTU abundance across all individual trees for each forest site is presented. Significant results (Bonferroni corrected a = 0.05/number of samples within each country) are indicated in bold print (Râşca α = 0.0045; Hainich α = 0.0056; Białowieża α = 0.00625; and North Karelia α = 0.008).
    GroupaPutative taxonomic assignmentbParameter estimate for tree species richness
    OTU_4ASydowia polyspora0.
    OTU_5AArthoniomycetes sp.
    OTU_12ATumularia sp.−0.100.12−0.830.405
    OTU_8ACelosporium larixicola0.000.23−0.020.988
    OTU_30AArthoniomycetes sp.0.530.351.530.125
    OTU_23ALeotiomycetes sp.0.360.142.520.01
    OTU_3ACladosporium sp.
    OTU_24ADothideomycetes sp.−0.720.16−4.401.07E-05
    OTU_26ACyphellophora sp.
    OTU_2AHelotiales sp.
    OTU_7AAureobasidium pullulans0.450.192.440.01
    OTU_2AHelotiales sp.−0.340.08−4.381.19E-05
    OTU_3ACladosporium sp.−0.090.07−1.250.213
    OTU_7AAureobasidium pullulans0.950.204.831.35E-06
    OTU_4ASydowia polyspora−0.310.20−1.590.111
    OTU_6ACeramothyrium sp.0.940.00413.70<2e-16
    OTU_8ACelosporium larixicola−0.160.21−0.750.451
    OTU_11ATaphrina carpini0.070.170.410.680
    OTU_18AEpicoccum nigrum0.200.111.920.055
    OTU_12ATumularia sp.
    OTU_49AAscochyta skagwayensis−0.070.35−0.190.852
    OTU_33ARamularia vizellae−0.250.21−1.170.241
    OTU_35ATrichomerium sp.0.770.233.440.001
    OTU_38ATaphrina vestergrenii−0.030.19−0.180.858
    OTU_17AChaetothyriales sp.
    OTU_39BSporobolomyces sp.−0.430.00−136.60<2e-16
    OTU_13BExobasidium sp.−0.090.30−0.300.767
    OTU_60AHortaea sp.−0.940.47−1.990.046
    OTU_61AChaetothyriales sp.0.120.0049.70<2e-16
    OTU_3ACladosporium sp.
    OTU_7AAureobasidium pullulans0.450.291.560.119
    OTU_5AArthoniomycetes sp.
    OTU_4ASydowia polyspora0.
    OTU_6ACeramothyrium sp.
    OTU_40AMetschnikowia sp.−0.490.37−1.330.183
    OTU_12ATumularia sp.−0.060.13−0.470.642
    OTU_15BExobasidium arescens−0.330.14−2.340.019
    North Karelia
    OTU_1BChrysomyxa ledi0.341.070.320.751
    OTU_4ASydowia polyspora−0.780.44−1.780.075
    OTU_3ACladosporium sp.0.060.380.150.882
    OTU_8ACelosporium larixicola−1.040.50−2.080.038
    OTU_7AAureobasidium pullulans−0.800.49−1.630.103
    OTU_22APhaeocryptopus gaeumannii−0.340.87−0.390.698
    OTU_9ADothideomycetes sp.−0.600.56−1.060.287
    OTU_5AArthoniomycetes sp.−0.030.81−0.040.966
    • a
      Group refers to taxonomic level of phylum, where A is Ascomycota and B is Basidiomycota.
    • b
      Putative taxonomic assignment of each OTU.
    Fungal community composition differed strongly between pairs of forests that were differentially processed and sequenced at different times, regardless of tree species diversity, namely between North Karelia and Hainich, and between Râşca and Białowieża (R2 = 0.28, P < 0.001; Supplementary Table 1). Tree species richness did not influence the fungal community composition across forests (Supplementary Table 1). Fungal community diversity, as determined by Hill numbers, was not significantly different between the forests either (data not shown).

    3.3 Fungal community from the Finland RNA study (metabolically active and total community analysis)

    Fungal community composition of each community type, i.e. active fungal community (RNA) and total fungal community (DNA), in North Karelia, Finland was distinct (R2 = 0.08, P = 0.009; Fig. 3Supplementary Table 1). However, the active (RNA) and total (DNA) fungal communities did not differ significantly among the tree species richness levels (P = 0.102 and P = 0.169, respectively, Supplementary Table 1) and this was further confirmed by the overlap in the community profile in the NMDS of the fungal communities from both the metabolically active and total community, respectively, along the tree richness gradient (data not shown). Fungal community diversity, as determined by Hill numbers, was tested between the community types, and separately for the metabolically active and total fungal communities. The active community and total community were equally diverse, and were not significantly different along the tree diversity gradient for either community (data not shown).
    Fig. 3. NMDS results of the fungal community composition of plot-level OTUs for each community type, i.e. active (RNA) and total (DNA) fungal community, in the North Karelia forest in Finland. Axes are arbitrary, represent three NMDS dimensions and scaled in units of Bray-Curtis dissimilarity. A stable solution was reached (stress value = 0.13). The colored ellipses mark the 95% confidence interval of the group centroids for either the active RNA community (orange) or total DNA community (blue). Lines link the samples pooled at the plot level to the group centroids.
    The most abundant OTUs from each of the community types varied in their relative read abundance, which was between 1% and 15% at the most extreme, and found in varying number of samples (Supplementary Table 3). Negative tree species richness effects were observed for OTU_16 (H. parviporum) in the active (RNA) community (Table 4). Tree species richness effects were also observed in the total (DNA) community. Here, OTU_11 (T. carpini) was negatively affected by the tree species richness, while positive effects were seen for OTU_20 (Eurotiomycetes sp.) and OTU_32 (Pleosporales sp.) (Table 4).
    Table 4. Generalized linear mixed-effects model results at the fungal species level for the dominant OTUs. Tree species richness effect on OTU abundance across all individual trees for the fungal community type, i.e. active (RNA) and total (DNA), in North Karelia, Finland. Significant results (Bonferroni corrected a = 0.05/number of samples in each community type) are indicated in bold print (RNA α = 0.0009; and DNA α = 0.0009).

    For further details log on website :

    No comments:

    Post a Comment