Microbial taxonomical composition in spruce phyllosphere, but not community functional structure, varies by geographical location

Sample sites and sampling

Leaves of spruce (Picea spp.) tree species were collected from three sites in Gansu Province, China (Fig. 1): in the Yuzhong campus of Lanzhou University in Lanzhou city (LDU, 104.16°E, 35.94°N), in Xiaolongshan forest farm in Tianshui City (XLS, 106.56°E, 34.12°N), and a site in the Xinglongshan forest in Majiasi (MJS, 104.01°E, 35.88°N). At each location, five spruce (Picea spp.) tree species of Picea abies, Picea crassifolia, Picea koraiensis, Picea likiangensis var. rubescens, and Picea wilsonii were selected. In total, 15 leaf samples were obtained at each location (three replicates per tree).

In the year 2006, spruce seedlings were bred in Xiaolongshan forest farm. The spruce species germinated from seeds, which were surface sterilized with a 0.50% potassium permanganate (KMnO₄) solution. Then the spruce seedlings were transplanted in the Yuzhong campus of Lanzhou University and the Xinglongshan forest in Majiasi in the year 2009 and 2014, respectively. All samples were collected in early October 2017 and all the trees were 11 years old. Weather conditions were warm and moist during the sampling.

The Yuzhong campus of Lanzhou University is in the northwest of the Loess Plateau, the middle part of Gansu Province, which is 46 km from Lanzhou city. The climate of this location is semiarid, and the mean annual precipitation is low (381.80 mm), nearly 78% of rainfall events are below 10 mm and 36% of precipitations are below five mm, and the mean annual temperature is 6.57 °C and droughts commonly occur. The soil in the region, which is two m deep, is the clay loam of Loess origin with a bulk density ranging from 1.38 to 1.45 g cm⁻³, and the water holding capacity of the soil is 21.18%. The results are supported by Wang et al. (2008).

The Xinglongshan forest in Majiasi is a unique island-like area where primitive forest is well preserved at the western edge of the Chinese Loess Plateau. The mean annual temperature is 4.10 °C, with a mean annual precipitation of 625 mm (Lei et al., 2010). Leached brown soil covers the area and was developed under semiarid climate under forest cover (Liu et al., 1990).

Xiaolongshan forest farm is located in the southeast of Gansu Province. The climate of this location is warm with humid and semi-humid continental monsoon. It is characterized by a long-term average annual temperature of 7–12 °C and an average annual precipitation of 600–900 mm. The soil in the area is dominated by mountain cinnamon soils and brown forest soils listed in the Chinese Soil Taxonomy. The information was obtained from the Xiaolongshan Forestry Bureau of Gansu Province.

Only “healthy” leaves that with no signs of disease or decay (such as browning or spotting) were selected assessed with the naked eye. Ten grams of fresh needles and shoots were collected from multiple canopy positions and mixed into a 100 ml plastic centrifuge tube containing 50 ml sterile phosphate buffer. The tube was agitated on vortex for 10 s and repeated for six times, and then the buffer was filtered through an EMD Millipore Sterivex-GV Polyvinylidene Fluoride 0.22 µm filter (Millipore, Billerica, MA, USA). Then, the filter was placed in a sterile screw-cap tube and frozen at −80 °C before processing.

DNA extraction, PCR amplification, and sequencing processing

The phylloshperic microbial DNA was extracted using the E.Z.N.A.^® Soil/Stool DNA Kit (Omega Bio-tek, Norcross, GA, USA) according to manufacturer’s protocols. Quantity and quality of the genomic DNA were checked on 0.80% (w/v) agarose gels. The V4 hypervariable region of the bacterial 16S rRNA gene was amplified with TransGen AP221-02 (TransStart FastPfu DNA Polymerase) using primers 515F-806R (Caporaso et al., 2011); PCRs contained 4 µl 5 × FastPfu buffer, two µl 2.50 mM deoxyribonucleoside triphosphates mix, 0.8 µl forward primer (five μM), 0.8 µl reverse primer (five μM), 0.40 µl FastPfu Polymerase, 0.20 µl 20 mg/mL bovine serum albumin solution and 10 ng Template DNA in a volume of 20 µl. The ITS1-ITS2 (White et al., 1990) was used to target partial fragments of the fungal ITS gene with TAKARA rTaq DNA polymerase (Takara Shuzo, Kyoto, Japan); PCRs contained 2 µl 10 × PCR buffer, 2 µl 2.50 mM deoxyribonucleoside triphosphates mix, 0.8 µl forward primer (five μM), 0.80 µl reverse primer (five μM), 0.2 µl rTaq Polymerase, 0.20 µl 20 mg/ml bovine serum albumin solution and 10 ng Template DNA in a volume of 20 µl. All samples were heated at 94 °C for 3 min in the first round of denaturation and were then subjected to 32 cycles of PCR consisting of 30 s at 95 °C, 30 s at 62 °C, and 45 s at 72 °C. Cycling was performed in an automated DNA thermal cycler (ABI GeneAmp^® PCR System 9700). After the last cycle, the samples were incubated for an additional 10 min at 72 °C. Amplicons were extracted from 2% (w/v) agarose gels and purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) according to the manufacturer’s instructions and quantified using QuantiFluor™ -ST (Promega, Madison, WI, USA). Sample libraries were pooled in equimolar and paired-end sequenced (2 × 300 bp) on an Illumina MiSeq (Illumina, San Diego, CA, USA) platform. Raw FASTQ files were quality-filtered by Trimmomatic and merged by FLASH with the following criteria: (i) The reads were truncated at any site receiving an average quality score <20 over a 50 bp sliding window. (ii) Sequences whose overlap being longer than 10 bp were merged according to their overlap with mismatch no more than two bp. (iii) Sequences of each sample were separated according to barcodes (exactly matching) and Primers (allowing for two nucleotides mismatching) and reads containing ambiguous bases were removed. Then the dataset was analyzed using QIIME (version 1.9.1, http://www.qiime.org). Operational taxonomic units (OTUs) with 97% similarity cutoff were clustered using UPARSE (version 7.1, http://drive5.com/uparse/) (Edgar, 2013). Chimeric sequences were identified and removed using UCHIME (Edgar et al., 2011). The taxonomy of each bacterial 16S rRNA gene sequence and fungal ITS sequence was analyzed by the RDP Classifier algorithm (http://rdp.cme.msu.edu/) against the Silva (Release119, http://www.arb-silva.de) and UNITE (Release 7.0, http://unite.ut.ee/index.php) databases at a confidence threshold of 70%, respectively. Sequences shorter or longer than the expected amplicons size and chimeras were removed. Cyanobacteria, mitochondria, archaeal and non-fungal reads were also filtered. After quality control, 1,921,840 and 2,915,075 valid sequences were clustered into 3,588 bacterial and 3,694 fungal OTUs at 97% sequence identity level, respectively (Table S1). There are 40,654 and 120 bacterial and fungal sequences could not be classified to any known bacterial and fungal phyla, respectively. Here, we defined OTUs with abundance higher than 5% in at least one of the replicates as abundant OTUs. The rarefaction analysis based on Mothur v.1.21.1 software was conducted to reveal the diversity indices, including the Chao1, Shannon, Simpson, and coverage indices. All sequences can be downloaded from the NCBI Sequence Read Archive database under the accession BioProject number PRJNA506625.

Tax4Fun and FUNGuild analyses

Tax4Fun is an inexpensive method to estimate the functional capacity of a microbial community (Aßhauer et al., 2015). This tool uses 16S rRNA gene profiles and then indirectly infers the abundance of functional genes. The results are KEGG orthologs (KOs) which wholly depend on the Silva ontology. Formatted sequences were clustered at 0.97 sequence similarity described above. Then, taxonomic information was assigned using Silva 123 downloaded from the Tax4Fun website (http://tax4fun.gobics.de/). An OTU table was created and fed to Tax4Fun R package. The Tax4Fun function was run with all default settings as described on the Tax4Fun website.

The fungal functional guild of the OTUs was inferred using FUNGuild v1.0 (http://funguild.org). This tool will infer organism trophic mode (Nguyen et al., 2016).

Data analysis

Alpha-diversity parameters were estimated (for a complete list of the calculated indexes see Table S2 in Supporting Information), and distance matrices (Bray–Curtis) of the samples were computed. A two-way ANOVA on the calculated alpha-diversity indices were performed in order to test differences in microbiome diversity between the different samples. Variations between different groups in Principal Coordinate Analysis (PCoA) were calculated by permutational multivariate analysis of variance (PERMANOVA) in R package “vegan,” using the function “adonis,” permutations = 999 (R Core Team, 2008). Bray–Curtis distance matrices used in the PCoA analysis were calculated by “dist.shared” command to test differences in the microbiome composition between different samples. A permutational multivariate analysis of dispersion (PERMDISP) for vegan’s betadisper was then run on this Bray–Curtis dissimilarity matrix to test the homogeneity of variance in community data (Anderson & Walsh, 2013). The statistical significance of the relationship between species dissimilarities and geographic locations was tested by ANCOVA analysis using SPSS. Analyses for the Venn diagram generation in microbial communities were performed using the Mothur v.1.21.1 suite of programs. The species typical of the different locations were determined using indicator species analysis (ISA). ISA was performed using the multipatt function implemented in the “indicspecies” package in R with 99,999 permutations and allowing combinations between habitats (De Cáceres, Legendre & Moretti, 2010). In this study, ISA were conducted with taxa grouped at the class and order level of taxonomic resolution, respectively. The p.adjust command was used to define a corrected critical p-value corresponding to a “false discovery rate” (FDR) of 0.05 (Benjamini & Hochberg, 1995). An FDR of 5% means that among all features called significant, 5% of these are truly null on average. Map created in ArcMap 9.3.

Phyllosphere microbial communities associated with spruce trees

To identify the phyllosphere microbiota associated with spruce trees, 16S rRNA and ITS sequences of 45 replicates across five species and three locations were obtained.

A heatmap representation of the abundant OTUs abundances is presented in Fig. S1. A total of 19 abundant OTUs were identified in 16S datasets. A total of 18 identified abundant OTUs were ascribed to four bacterial phyla (Proteobacteria, Actinobacteria, Bacteroidetes, and Firmicutes), while one unclassified at phylum level. A total of 24 abundant OTUs were identified in ITS datasets. At a finer taxonomic level, 20 identified abundant OTUs were ascribed to five fungal orders (Pleosporales, Chaetothyriales, Taphrinales, Capnodiales, and Dothideales), while four were unclassified at the order level. Some OTUs were shared between different species from geographically distant samples, suggesting the associated microbial taxa may be widespread.

In the bacterial communities, the most abundant bacterial phylum was the Proteobacteria, with an average abundance of 54.15% of the reads, followed by Bacteroidetes (20.32%), and Actinobacteria (10.92%). The remaining 14.62% of reads belonged to other bacterial phyla or were unclassified (unable to be taxonomically assigned, given the training set). In the fungal communities, the most abundant fungal order was the Pleosporales, with an average abundance of 33.2% of the reads, followed by Capnodiales (17.1%), Chaetothyriales (8.5%), and Taphrinales (7.8%). The four orders belonged to Ascomycota, which accounted for 66.6% of the reads.

Community dissimilarity estimation by diversity indices

The two-way ANOVA on the alpha-diversity measures (Chao 1 index (C), Shannon diversity index (H′), and Simpson diversity index (λ); Fig. 2; Table S3) of both bacterial and fungal communities evaluating all the factors considered in the experimental design and their interactions. The results revealed that the geographical location, the host plant species and the interaction of these two factors had the most influence on microbial alpha diversities. Comparing different samples (LDU, MJS, and XLS) in the bacterial communities, only the richness was highly significant with the location effect (F = 283.34, p < 0.001; Table S3), the host tree effect (F = 10.79, p < 0.001; Table S3), and the interaction effect of the location and host tree (F = 3.97, p = 0.003; Table S3), respectively, whereas for the fungal communities, the richness, the evenness, and the diversity were all significant with these factors at different degrees (Table S3). Host tree species had significant, albeit small, influences on the fungal community diversity (F = 2.89, p < 0.05; Table S3) and evenness (F = 3.50, p < 0.05; Table S3). The interaction between the sample location and the host tree species showed a highly significant effect of the fungal community for all the analyzed indexes (C, F = 7.33, p < 0.001; H′, F = 5.57, p < 0.001; λ, F = 5.04, p = 0.001; Table S3), while only the OTUs index was significant in the bacterial community (D, p < 0.001; H′, p = 0.78; λ, p = 0.15; Table S3).

Similarity in community composition depends upon taxonomical OTUs and functional groups

Differences in the microbial OTUs datasets both for the bacterial (Fig. 3A) and fungal (Fig. 3B) communities between all samples were visualized in a PCoA plot, where a clear separation of samples attributed to their geographic locations was evidenced. The significance of these differences was proven by PERMANOVA (Table S4). PERMDISP test demonstrated that variance is homogeneous in different sample groups for bacterial community, bacterial KO, fungal community, fungal guilds (p > 0.05) similar in all groups (Table S5). Leaves sampled as little as 14 km apart were found to harbor significantly different taxonomical communities (between LDU and MJS), and the geographic variability on same tree species exceeded the variability in community composition between different trees sampled at one site. Regarding the functional datasets of 16S (Fig. 3C), there is no clear cluster by the geographic location than there was by the taxonomic composition. But the functional composition of the fungal assemblage via FUNGuild designations appears to be more closely tied to the geographical location (Fig. 3D). It could be that this clustering phenomenon is caused by the shortage of the software itself. Because the FUNGuild is a tool of functional prediction which is a simple way to sort large sequence pools into ecologically meaningful categories, in an organism trophic mode not functional gene. Unlike Tax4Fun, the FUNGuild-based methods are inefficient to predicate various functional genes of microorganisms due to the technical shortage. The shortage of FUNGuild can be improved by metagenomic analysis, which is the key to know more about the energy and material metabolism of microbes and communities. But it is less cost-ineffective than FUNGuild, a tool that can be used to taxonomically parse fungal OTUs by the ecological guild independent of the sequencing platform or the analysis pipeline.

Venn diagrams of taxonomical OTUs and functional groups

To evaluate the distribution of OTUs and functional groups among the different sampling locations, Venn diagrams were constructed (Figs. 4A, 4B, 4E and 4F). This showed that 40.74%, 67.65%, and 54.29% of the bacterial OTUs (Fig. 4A) and 37.3%, 24.6%, and 30.2% of the fungal OTUs (Fig. 4E) were common (area in gray) to the three locations (LDU, XLS, and MJS, respectively). However, the proportion of the overlapped functional groups was much higher than those shared by OTUs. In 16S datasets, the overlapped KOs were 99.44% (LDU), 99.78% (XLS), and 99.84% (MJS), respectively (Fig. 4B). Similarly, the overlapped functional guilds were 81.8% (LDU), 73.8% (XLS), and 78.9% (MJS), respectively (Fig. 4F). Among the five-tree species, the Venn Diagrams showed that 67.11% ± 4.62% of the bacterial OTUs (Fig. 4C) and 56.73% ± 2.96% fungal OTUs (Fig. 4G) were shared, respectively. While the overlapped bacterial KOs was up to 99.77% ± 0.06% (Fig. 4D), and the overlapped functional guilds was up to 86.56% ± 1.57% (Fig. 4H).

Indicator species analysis

Indicator species analysis identifies the bacterial and fungal taxa that are significantly unique to each location. A total of 54 bacterial classes and 58 fungal orders were found with a significant preference (q < 0.05) for the three sampling sites, respectively (Tables S6 and S7). In 16S datasets, the 54 classes represented 2.36% of the total bacterial sequences. In ITS datasets, the 58 orders represented 4.09% of the total fungal sequences. There were 20, 12, and 40 indicator fungal taxa in LDU, MJS, and XLS, respectively (Fig. S2B). And there were 53, eight, and five indicator bacterial taxa in LDU, MJS, and XLS, respectively (Fig. S2A). Nearly all the relative abundance of the indicator taxa was low; except two taxa (the classes of Acidobacteriia and Gemmatimonadetes) were more than 1% in bacterial community (Fig. S2A) and two taxa (the orders of Sebacinales and Pezizales) were more than 2% in fungal community (Fig. S2B). Taxonomic names such as “other” and “unidentified” (has no assignment of affiliation at the relevant taxonomic level) were not shown in the bubble plot (Fig. S2).