Published Date
Diem Nguyen ,
Johanna Boberg
Katarina Ihrmark
Elna Stenström
Jan Stenlid
Diversity gradient
FunDivEUROPE
Forests
Fungal community
Norway spruce needles
Tree species richness
For further details log on website :
http://www.sciencedirect.com/science/article/pii/S1226086X16302301
doi:10.1016/j.funeco.2016.07.003
Open Access, Creative Commons license
Author
Received 18 November 2015. Revised 1 July 2016. Accepted 6 July 2016. Available online 5 August 2016. Corresponding Editor: Marie Louise Davey
Abstract
Foliar fungal species are diverse and colonize all plants, though whether forest tree species composition influences the distribution of these fungal communities remains unclear. Fungal communities include quiescent taxa and the functionally important and metabolically active taxa that respond to changes in the environment. To determine fungal community shifts along a tree species diversity gradient, needles of Norway spruce were sampled from trees from four mature European forests. We hypothesized that the fungal communities and specific fungal taxa would correlate with tree species diversity. Furthermore, the active fungal community, and not the total community, would shift along the tree diversity gradient. High-throughput sequencing showed significant differences in the fungal communities in the different forests, and in one forest, tree diversity effects were observed, though this was not a general phenomenon. Our study also suggests that studying the metabolically active community may not provide additional information about community composition or diversity.
Keywords
1 Introduction
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, 1998, Hantsch 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, 1978, Petrini, 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, 2000, Arnold 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, 1995, Vacher 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., 2005, Jumpponen and Jones, 2009, Jumpponen et al., 2010, Cordier et al., 2012b, Menkis 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., 2004, Baldrian 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. 1, Table 1). Sampling was conducted over approximately 2 weeks in 2012 or 2013 for each of the four forests during the vegetation period.
Table 1. Sampling site locations, forest type and composition of forest mixture with respect to Norway spruce.
Forest site | Râşca | Hainich | Białowieża | North Karelia |
---|---|---|---|---|
Romania | Germany | Poland | Finland | |
Forest type | Mountainous beech | Temperate Beech | Hemi-boreal | Boreal |
Species pool | Picea 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 |
Location | 47.3° N, 26.0° E | 51.1° N, 10.5° E | 52.7° N, 23.9° E | 62.6° N, 29.9° E |
Study area size (km × km) | 5 × 5 | 15 × 10 | 30 × 40 | 150 × 150 |
Mean annual temperature, °C | 6.8 | 6.8 | 6.9 | 2.1 |
Mean annual precipitation, mm | 800 | 775 | 627 | 700 |
Sampling period | July 2013 | July 2012 | July–August 2013 | August 2012 |
Number of plots (trees) with Picea abies | ||||
Monoculture | 2 (12) | 2 (12) | 2 (12) | 4 (24) |
Two species | 5 (15) | 3 (8) | 5 (15) | 8 (24) |
Three species | 5 (15) | 3 (8) | 6 (18) | 4 (12) |
Four species | 3 (9) | 2 (6) | 8 (24) | |
Five species | 2 (6) | |||
Total | 15 (51) | 10 (34) | 23 (75) | 16 (60) |
Number of trees sampled and sequenced | ||||
51 | 34 | 75 | 60/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 (http://scata.mykopat.slu.se). 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, http://blast.ncbi.nlm.nih.gov). 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., 2006, Hibbett et al., 2007, Boehm et al., 2009, Hodkinson and Lendemer, 2011, Rosling et al., 2011, Zhang et al., 2011, Boonmee et al., 2014, Ertz 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, 1987; Hosmer 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 skagwayensis, Aureobasidium pullulans(ubiquitous yeast-like fungus), Celosporium larixicola, Ceramothyrium carniolicum, Chrysomyxa ledi (rust pathogen of Norway spruce), Epicoccum nigrum, Exobasidium arescens, Heterobasidion parviporum (root and butt rot pathogen of Norway spruce), Phaeocryptopus gaeumannii (pathogen of Douglas-fir (Pseudotsuga menziesii)), Ramularia vizellae (plant pathogen), Sydowia polyspora, Taphrina 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 | ||||
---|---|---|---|---|---|
RNA | DNA | ||||
Ascomycota | 215 | 11991 | 34967 | ||
Basidiomycota | 68 | 5077 | 14092 | ||
Fungus | 2 | 141 | 258 | ||
Glomeromycota | 1 | 7 | 8 | ||
Total | 286 | 17216 | 49325 | ||
Number of OTUs | Number of reads | ||||
Only in RNA | Only in DNA | Only in RNA | Only in DNA | ||
Ascomycota | 4 | 13 | 71 | 297 | |
Basidiomycota | 2 | 4 | 26 | 49 | |
Fungus | 0 | 0 | |||
Glomeromycota | 1 | 0 | 13 | ||
Total | 7 | 17 | 110 | 346 | |
Number of OTUs | Number of reads | ||||
RNA total | DNA total | RNA total | DNA total | ||
Ascomycota | 219 | 228 | 12062 | 35264 | |
Basidiomycota | 70 | 72 | 5103 | 14141 | |
Fungus | 2 | 2 | 141 | 258 | |
Glomeromycota | 2 | 1 | 20 | 8 | |
Total | 293 | 303 | 17326 | 49671 |
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. 2, Supplementary 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).
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).
Groupa | Putative taxonomic assignmentb | Parameter estimate for tree species richness | ||||
---|---|---|---|---|---|---|
Estimate | StdError | Z-value | P | |||
Râşca | ||||||
OTU_4 | A | Sydowia polyspora | 0.02 | 0.20 | 0.09 | 0.927 |
OTU_5 | A | Arthoniomycetes sp. | 0.08 | 0.14 | 0.60 | 0.547 |
OTU_12 | A | Tumularia sp. | −0.10 | 0.12 | −0.83 | 0.405 |
OTU_8 | A | Celosporium larixicola | 0.00 | 0.23 | −0.02 | 0.988 |
OTU_30 | A | Arthoniomycetes sp. | 0.53 | 0.35 | 1.53 | 0.125 |
OTU_23 | A | Leotiomycetes sp. | 0.36 | 0.14 | 2.52 | 0.01 |
OTU_3 | A | Cladosporium sp. | 0.03 | 0.14 | 0.21 | 0.84 |
OTU_24 | A | Dothideomycetes sp. | −0.72 | 0.16 | −4.40 | 1.07E-05 |
OTU_26 | A | Cyphellophora sp. | 0.08 | 0.17 | 0.47 | 0.64 |
OTU_2 | A | Helotiales sp. | 0.15 | 0.21 | 0.71 | 0.48 |
OTU_7 | A | Aureobasidium pullulans | 0.45 | 0.19 | 2.44 | 0.01 |
Hainich | ||||||
OTU_2 | A | Helotiales sp. | −0.34 | 0.08 | −4.38 | 1.19E-05 |
OTU_3 | A | Cladosporium sp. | −0.09 | 0.07 | −1.25 | 0.213 |
OTU_7 | A | Aureobasidium pullulans | 0.95 | 0.20 | 4.83 | 1.35E-06 |
OTU_4 | A | Sydowia polyspora | −0.31 | 0.20 | −1.59 | 0.111 |
OTU_6 | A | Ceramothyrium sp. | 0.94 | 0.00 | 413.70 | <2e-16 |
OTU_8 | A | Celosporium larixicola | −0.16 | 0.21 | −0.75 | 0.451 |
OTU_11 | A | Taphrina carpini | 0.07 | 0.17 | 0.41 | 0.680 |
OTU_18 | A | Epicoccum nigrum | 0.20 | 0.11 | 1.92 | 0.055 |
OTU_12 | A | Tumularia sp. | 0.06 | 0.21 | 0.29 | 0.770 |
OTU_49 | A | Ascochyta skagwayensis | −0.07 | 0.35 | −0.19 | 0.852 |
OTU_33 | A | Ramularia vizellae | −0.25 | 0.21 | −1.17 | 0.241 |
OTU_35 | A | Trichomerium sp. | 0.77 | 0.23 | 3.44 | 0.001 |
OTU_38 | A | Taphrina vestergrenii | −0.03 | 0.19 | −0.18 | 0.858 |
OTU_17 | A | Chaetothyriales sp. | 0.02 | 0.17 | 0.11 | 0.913 |
OTU_39 | B | Sporobolomyces sp. | −0.43 | 0.00 | −136.60 | <2e-16 |
OTU_13 | B | Exobasidium sp. | −0.09 | 0.30 | −0.30 | 0.767 |
OTU_60 | A | Hortaea sp. | −0.94 | 0.47 | −1.99 | 0.046 |
OTU_61 | A | Chaetothyriales sp. | 0.12 | 0.00 | 49.70 | <2e-16 |
Białowieża | ||||||
OTU_3 | A | Cladosporium sp. | 0.07 | 0.09 | 0.80 | 0.422 |
OTU_7 | A | Aureobasidium pullulans | 0.45 | 0.29 | 1.56 | 0.119 |
OTU_5 | A | Arthoniomycetes sp. | 0.20 | 0.18 | 1.08 | 0.281 |
OTU_4 | A | Sydowia polyspora | 0.04 | 0.13 | 0.27 | 0.785 |
OTU_6 | A | Ceramothyrium sp. | 0.15 | 0.15 | 1.02 | 0.308 |
OTU_40 | A | Metschnikowia sp. | −0.49 | 0.37 | −1.33 | 0.183 |
OTU_12 | A | Tumularia sp. | −0.06 | 0.13 | −0.47 | 0.642 |
OTU_15 | B | Exobasidium arescens | −0.33 | 0.14 | −2.34 | 0.019 |
North Karelia | ||||||
OTU_1 | B | Chrysomyxa ledi | 0.34 | 1.07 | 0.32 | 0.751 |
OTU_4 | A | Sydowia polyspora | −0.78 | 0.44 | −1.78 | 0.075 |
OTU_3 | A | Cladosporium sp. | 0.06 | 0.38 | 0.15 | 0.882 |
OTU_8 | A | Celosporium larixicola | −1.04 | 0.50 | −2.08 | 0.038 |
OTU_7 | A | Aureobasidium pullulans | −0.80 | 0.49 | −1.63 | 0.103 |
OTU_22 | A | Phaeocryptopus gaeumannii | −0.34 | 0.87 | −0.39 | 0.698 |
OTU_9 | A | Dothideomycetes sp. | −0.60 | 0.56 | −1.06 | 0.287 |
OTU_5 | A | Arthoniomycetes sp. | −0.03 | 0.81 | −0.04 | 0.966 |
- aGroup refers to taxonomic level of phylum, where A is Ascomycota and B is Basidiomycota.
- bPutative 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. 3, Supplementary 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).
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).
For further details log on website :
http://www.sciencedirect.com/science/article/pii/S1226086X16302301
No comments:
Post a Comment