Seasonal dynamics of microbial diversity in the rhizosphere of Ulmus pumila L. var. sabulosa in a steppe desert area of Northern China

Sampling sites

The spatial heterogeneity of soil physical and chemical properties and soil enzyme activities are low in OSL (Table S1; Fig. S1). The study was conducted at Sanggen Dalai (latitude, 42°40′N; longitude, 115°57′E; 1,300 m AMSL), which is located in the center of OSL. The open elm woodlands of this area are the most well-preserved communities within the OSL (Wang et al., 2015). The area is characterized as a temperate semi-arid continental climate zone, with an annual average temperature of 1–4 °C, an average frost-free period of 105 days, annual sunshine time of 3,200 h, an annual average wind speed of 4.2 m/s, and annual precipitation of 320 mm. In this region, 90% of the dunes are fixed and semi-fixed dunes, while shifting dunes only account for a small proportion. However, increasing numbers of semi-fixed dunes have transitioned to shifting dunes due to drought and overgrazing. The zonal vegetation is that of a steppe desert. The vegetation on fixed and semi-fixed dunes is more productive, and they have been thus used as sandy pastures. Elms in the area are distributed on the windward slopes, leeward slopes, and lowlands of the sand dunes. The present study selected four lowlands with elms as fixed sampling sites (Fig. 1). Sampling on the private pasture was approved by the owner of the pasture, Ying Tao.

Soil sample collection

Soil sampling was conducted in May (spring), August (summer), and October (autumn) of 2017. Low temperatures and thick snow cover in winter precluded soil sampling in the winter. Three individual elms were randomly selected for sampling from each of the four sampling sites (Fig. 1). At each site, the three sampled trees were 10 to 20 m apart from one another. Rhizosphere soils were collected at depths of 5–30 cm extending east, south, west, and north, 1m from the trunk of each tree. Root systems were carefully excavated using a spade, and loosely adhering soils on the roots were shaken off and discarded. The root-adherent soil particles were then collected. During each season, the soil samples from each sampling site (3 trees × 4 direction = 12 samples per site) were pooled in equal proportions to obtain a composite soil sample for each site. At each location, the same three trees were sampled at the three sampling periods. Thus, a total of 12 soil composite soil samples (i.e., 4 sites ×3 seasons) were collected. Each soil sample was subdivided into three components, with one stored at −80 °C for use in determining soil microbial community diversity, while the other was stored at 4 °C for use in soil enzyme activity determination, and the third was dried for use in determining soil physicochemical properties.

Determination of soil physicochemical properties

Soil physiochemical properties were measured by referencing protocols described in Soil Agricultural Chemistry analysis (Bao, 2000). Total N was determined with the Kjeldah method, while available N was determined using alkaline hydrolysis diffusion. Total P was determined using NaOH-Mo-Sb colorimetry and available P was determined via NaHCO₃ extraction-Mo-Sb colorimetry. Available K was determined by flame photometry. Organic matter content was determined using a K₂Cr₂O₇ oxidation-external heating titration. pH was determined by potentiometry while moisture content was determined by measurement before and after oven drying. Lastly, soil temperature was determined using a miniature electronic temperature recorder (DS1921G; WDS). During the sampling month, continuous observation was made for one month.

DNA extraction

An SDS (sodium dodecyl sulfate)-based method was used to extract total DNA from the rhizosphere soil samples (Zhou, 1996). Soil samples (5 g) were weighed and ground to a powder in a ceramic mortar using liquid nitrogen. Samples were then transferred to 50 mL centrifuge tubes. DNA extraction buffer was added (13.5 mL; 100 mM Tris–HCl [pH 8.0], 100 mM sodium EDTA [pH 8.0], 10 mM sodium phosphate [pH 8.0], 1.5 M NaCl, and 1% CTAB) in addition to 100 µL of proteinase K (10 mg/mL). The tubes were then placed on a shaker with rotation at 225 r/min for 30 min at 37 °C. Following incubation, 1.5 mL of 20% SDS was added, and the samples were incubated in a 65 °C water bath for 2 h with shaking of tubes every 15 min to ensure adequate mixing. After incubation, the mixtures were centrifuged at room temperature at 6,000 r/min for 10 min. The supernatants were then transferred to new 50 mL centrifuge tubes, while the remaining soil pellets were subjected to two additional DNA extraction rounds using 4.5 mL of DNA extraction buffer and 0.5 mL of 20% SDS, with mixing for 10 s, incubation at 65 °C for 10 min, and centrifugation as described above. Supernatants from triplicate extractions were pooled together and mixed with an equal volume of chloroform-isoamyl alcohol (24:1). The samples were mixed by inversion and then centrifuged at 6000 r/min for 15 min. The aqueous phase was then transferred to a new 50 mL centrifuge tube. Isopropyl alcohol (0.6 × volume) was added to the samples, mixed by inversion, and allowed to precipitate at −20 °C overnight. The mixture was then centrifuged at 11,000 r/min for 20 min at 4 °C. Visible black or brown crude DNA extracts precipitated during this process. The precipitate was then transferred to a sterile 2.5 mL centrifuge tube. The crude DNA extract was washed with 70% cold ethanol, placed in a fume hood to allow ethanol volatilization, and 100 µL of sterile deionized water were finally added to dissolve DNA in water. DNA concentrations and purities were evaluated with 1% agarose gel electrophoresis.

PCR amplification and purification

The V4 hypervariable regions of bacterial 16S rRNA genes were amplified using the forward primer 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and reverse primer 806R (5′-GTGCCAGCMGCCGCGGTAA-3′) to construct bacterial community libraries for HiSeq Illumina sequencing (Evans et al., 2014). Taxonomic coverage of 16S rRNA genes was evaluated in https://www.arb-silva.de/search/testprime/ (Klindworth et al., 2013), which yielded an overall taxonomic coverage of 92.4%. Likewise, the ITS1 region of the fungal ITS was amplified using the forward primer ITS5-1737F (5′-GGAAGTAAAAGTCGTAACAAGG-3′) and reverse primer ITS2-2043R (5′-GCTGCGTTCTTCATCGATGC -3′) (Lu et al., 2013). ITS1 is the common prime target for the evaluation of fungal diversity through deep sequencing, but might over estimate the fungal diversity because of its variable length (Yang et al., 2018). All PCR reactions were carried out in 30 µL reaction volumes with 15 µL of Phusion^® High-Fidelity PCR Master Mix (New England Biolabs, Ipswich, MA, USA), 0.2 µM each of forward and reverse primers, and about 10 ng of template DNA. The following PCR procedure was used for bacterial 16S rRNA gene amplification: an initial denaturation at 98 °C for 1 min, 30 cycles of (98 °C for 10 s, 50  °C for 30 s, 72 °C for 30 s), and a final extension at 72 °C for 5 min. The PCR protocol for fungal ITS comprised an initial denaturation at 95  °C for 2 min, and then 30 cycles of (95  °C for 30 s, 55 °C for 30 s, and 72  °C for 30 s), with a final extension at 72  °C for 5 min. An equal volume of 1 × loading buffer (containing SYBR green) was mixed with the PCR products and the mixtures were subject to electrophoresis on a 2% agarose gel for detection. PCR products were then mixed in equimolar ratios, and the resultant pooled PCR products were purified using a GeneJETTM Gel Extraction Kit (Thermo Scientific, Waltham, MA, USA). After purification, the PCR product pools were sequenced on the Illumina Hiseq PE250 platform.

Bioinformatics and statistical analysis of sequence data

Samples were demultiplexed using unique barcode sequences for samples, followed by cleaving of barcode and PCR primer sequences. Paired end reads from PE250 two-terminal sequencing on the HiSeq platform were then joined using the FLASH software (V1.2.7, http://ccb.jhu.edu/software/FLASH/) (Magoč & Salzberg, 2011), and the joined sequences were considered as ‘Raw Tags’. The QIIME software suite (V1.9.1, http://qiime.org/scripts/split_libraries_fastq.html) (Caporaso et al., 2010) pipeline was used to quality control and filter raw tags to obtain ‘clean tags’ representing high-quality sequence reads. Clean tag sequences were compared against the Gold database (http://drive5.com/uchime/uchime_download.html) (Haas et al., 2011) using the UCHIME Algorithm (http://www.drive5.com/usearch/manual/uchime_algo.html) (Edgar et al., 2011) to remove chimeric sequences and finally obtain ‘effective tags’. The Uparse software suite (V7.0.1001, http://drive5.com/uparse/) (Edgar, 2013) was used to cluster effective tags into operational taxonomic units (OTUs) that were defined at 97% nucleotide identities. Taxonomic classification of representative bacterial OTUs was conducted using the SILVA SSUrRNA database (http://www.arb-silva.de/) based on the Mothur algorithm (Pruesse et al., 2007). Likewise, the BLAST method (http://qiime.org/scripts/assign_taxonomy.html) (Altschul, 1990) within QIIME (V1.9.1) was used to taxonomically annotate fungal representative OTU sequences by comparison against the UNITE database (https://unite.ut.ee/) (Johnson, Warburton & Mills, 2010). The MUSCLE (V3.8.31, http://www.drive5.com/muscle/) (Edgar, 2004) aligner was used to generate multi-sequence alignments for all representative OTU sequences. Finally, the number of sequences per sample were rarefied to equivalent levels (48,262 sequences for Bacteria and 48,027 for Fungi). Subsequent α-diversity and β-diversity analyses were conducted on the rarefied diversity tables. Alpha diversity indices were calculated using QIIME (V1.9.1) and rarefaction curves were drawn using the R software package (R Core Team, 2013).

One-way ANOVA (analysis of variance) and Bonferroni corrections for multiple comparisons were conducted in the SPSS 20.0 software suite. QIIME (V1.9.1) was also used to calculate Unifrac distances among samples to evaluate β-diversity. The Adonis function in the vegan package for R was used to conduct PERMANOVA analysis. This is a multivariate ANOVA method that seeks to evaluate the variation between samples, as indicated by distance matrices (Anderson, 2010). The vegan package for R was used to generate NMDS (non-metric multi-dimensional scaling) plots (Kruskal, 1964) and construct RDA (redundancy discriminant analysis) biplots and evaluate the sources underlying sample variation and their relative importance via exploration of ‘explanatory’ and ‘response’ variables (Oksanen et al., 2013; Kenkel et al., 2002). Lastly, the LEfSe software package was used to conduct LEfSe analysis, using an LDA score of 4 to filter results, for identifying taxonomic groups that are consistent with experimental treatments and estimating the associated effect size through class comparisons (Segata et al., 2011).

Functional analysis of bacterial and fungal communities

To predict variation in the inferred functions of bacterial communities within elm rhizospheres, PICRUSt analysis was employed. A closed-reference OTU table was first generated in QIIME and compared against the KEGG functional database to obtain predicted functions that were inferred for the bacterial communities. Individual analysis steps are outlined in the online analysis platform (http: //picrust.github.io/picrust/; Langille et al., 2013). FUNGuild (https://github.com/UMNFuN/FUNGuild; Nguyen et al., 2015) was used to identify the potential ecological roles of fungi within the communities. To evaluate the level of difference in inferred functionalities, one-way ANOVA (analysis of variance) analysis and Bonferroni corrections for multiple comparisons were used with the inferred function dataset using the SPSS 20.0 software package (SPSS, Chicago, IL, USA).

Seasonal dynamics of soil physicochemical properties

The soils were weakly alkaline, exhibiting a pH range of 7.38–7.43 without significant differences among seasons (Table 1). The moisture content (MC) of soils over all three seasons was lower than 10%. The maximum MC was observed in spring (5.51%) and the minimum was observed in autumn (1.27%), which was significantly lower (p < 0.017). Surface soil (5–30 cm) temperature dynamics were consistent with those of air temperatures, wherein temperatures in summer were greater than in spring and autumn. Total P did not significantly differ among seasons, while available phosphorus levels in summer were significantly higher than those of spring and autumn (p < 0.017). Available K content was highest in spring and significantly differed from K contents in the autumn (p < 0.017). Lastly, seasonal variation in total nitrogen, effective nitrogen, and organic matter was not apparent (Table 1).

Response of microbial community α-diversity in elm rhizosphere soils across seasons

A total of 909,174 and 945,531 16S rRNA gene and ITS sequences for bacterial and fungal communities were obtained, respectively (Table S2; Table S3). The distribution of OTUs that were defined at the 3% nucleotide dissimilarity level were then identified among samples. A total of 7,676 16S rRNA gene OTUs were observed among all samples from the three seasons. Among these, 4,933, 6,269, and 5,430 OTUs were observed in the spring, summer, and autumn samples, respectively. Of these, 3,627 OTUs were common to all samples, which accounted for 47.25% of the total OTUs (Fig. 2A). In addition, 3,582 ITS sequence OTUs were observed among all samples from the three seasons. Among these, 2,449, 2,489, and 1,618 OTUs were observed in the samples from spring, summer, and autumn, respectively. A total of 986 fungal ITS OTUs were common to all samples, which accounted for 27.53% of the total fungal OTUs (Fig. 2B).

After rarefying the number of sequences per sample to 48,262 and 48,027 for bacterial and fungal communities, respectively, an average of 3,579 and 1,043 OTUs were observed for the communities, respectively. The rarefaction curves and species accumulation boxplots for both bacterial and fungal communities generally reached asymptotic levels, indicating that the sample numbers and the sequencing effort applied here was adequate to capture most of the bacterial and fungal diversity in these communities (Fig. 3; Fig. S2).

Four indices of α-diversity were analyzed including the Shannon, Simpson, Chao1, and ACE. In addition, we also evaluated Good’s coverage index that provides an estimate of diversity captured by the sequencing effort, which yielded estimates of 98.1% and 99.5% for the bacterial and fungal communities, respectively. The four diversity indices followed similar trends wherein bacterial diversity was highest in spring rhizospheres, followed by summer and autumn (Table 2). One-way ANOVA (analysis of variance) of the Shannon index values indicated a significant difference in diversity among the three seasons. Subsequent multiple comparisons (with Bonferroni correction) indicated that autumn communities were significantly different from those in the spring (p < 0.017), while the diversity of communities in the summer were not significantly different from diversity in the spring and autumn communities. In addition, one-way ANOVA revealed significant differences for Simpson indices among soils from the three seasons. Subsequent multiple comparisons indicated that the spring community diversity was significantly different from that in the fall (p < 0.017), while summer community diversity did not significantly differ from community diversity in the spring and autumn. Lastly, the seasonal differences in the Chao1 and ACE diversity indices were similar, where communities in the spring exhibited significantly different diversities compared to those from the summer and autumn (p < 0.017), but significant differences were not observed between summer and autumn rhizosphere communities (Table 2). The α-diversity indices for fungal communities all followed trends wherein diversity was highest in summer, followed by diversity in the spring and then autumn (Table 2). One-way ANOVA of Shannon and Simpson indices indicated lack of significant differences among seasons. In contrast, one-way ANOVA indicated significant differences in the Chao1 and ACE indices among the communities from the three seasons. Subsequent multiple pairwise comparisons indicated that autumn community diversity was significantly different from the diversity of the spring (p < 0.017), while the community diversity from summer did not significantly differ from the diversity observed in the spring and autumn (Table 2).

Variation in microbial community taxonomic composition of elm rhizospheres among seasons

The rhizosphere-associated bacterial OTUs were associated with 40 total phyla, of which 12 were in >0.6% average relative abundance, and are shown in Fig. 4A. Most of the rhizospheric bacterial sequences (95.1%) belonged to the Proteobacteria, Actinobacteria, Acidobacteria, Firmicutes, Bacteroidetes, Verrucomicrobia, Gemmatimonadete, Chloroflexi, and Planctomycetes. Proteobacteria was the most abundant phylum in nearly all samples, with the exception of GSP4, GSU1, and GSU4. Acidobacteria was the most abundant phylum in GSP4, while Actinobacteria dominated the GSU1 and GSU4 communities. Of the above phyla, Proteobacteria, Actinobacteria, Acidobacteria, and Firmicutes represented the largest components of each rhizosphere soil community, together comprising at least 67% of the total bacterial communities within each sample (Fig. 4A). To assess taxonomic variation among soils at a finer taxonomic resolution, bacterial classes that were in >3% relative abundance of the rhizosphere communities were further analyzed (Fig. 4B). Four classes (Alphaproteobacteria, Gammaproteobacteria, Betaproteobacteria, and Deltaproteobacteria) were identified within the Proteobacteria, two (unidentified Actinobacteria, and Thermoleophilia) from the Actinobacteria, two (Subgroup_6, and Blastocatellia) from the Acidobacteria, two (Bacilli and Clostridia) from the Firmicutes, one (Spartobacteria) from the Verrucomicrobia, one (Spartobacteria) from the Bacterioidetes, and one (unidentified Gemmatimonadetes) from the Gemmatimonadetes. The most abundant classes among all rhizosphere communities were the Alphaproteobacteria, unidentified Actinobacteria and Bacilli, which together comprised at least 24.93% of the total bacterial communities in each sample. The next most abundant classes were the ‘Subgroup_6, of the Acidobacteria, Thermoleophilia, Gammaproteobacteria, and Betaproteobacteria, which comprised at least 16.97% of the total bacterial communities within each sample. Clostridia populations were noticeably more abundant in the spring samples, comprising at least 5.8% of the bacterial communities in these samples (Fig. 4B).

A total of 28 fungal classes were identified in the rhizosphere soils that belonged to six phyla including the Ascomycota, Basidiomycota, Zygomycota, Chytridiomycota, Glomeromycota, and Neocallimastigomycota. All of these phyla were present in every sample, with the exception of the Neocallimastigomycota that were only present in spring and summer samples. The Ascomycota overwhelmingly dominated the communities for all samples, and comprised at least 80.53% of the total fungal communities for every sample (Fig. 4C). To assess taxonomic variation among soils at a finer taxonomic resolution, fungal classes in >0.6% relative abundances in soils were analyzed (Fig. 4D). Six classes (Sordariomycetes, Dothideomycetes, Eurotiomycetes, Leotiomycetes, Incertae_sedis_Ascomycota, and Pezizomycetes) were identified from the Ascomycota, two (Agaricomycetes and Tremellomycetes) from the Basidiomycota, one (Incertae_sedis_Zygomycota) from the Zygomycota, one (Glomeromycetes) from the Glomeromycota, and one (Chytridiomycetes) from the Chytridiomycota. The most abundant classes among all rhizosphere samples were the Sordariomycetes, Dothideomycetes, and Eurotiomycetes, which together comprised at least 52.58% of the total fungal populations in each of the samples. In particular, the classes Leotiomycetes and Incertae_sedis_Zygomycota together comprised at least 6.85% of the total fungal populations in each of the samples (Fig. 4D).

β-diversity of elm rhizosphere microbial communities among seasons

The NMDS ordination of the elm rhizosphere bacterial communities exhibited a stress value of 0.033 (Fig. 5A), which indicated appropriate representation of community compositional dissimilarities (Kruskal, 1964). Bacterial communities segregated by season of origin, wherein NMDS axis 1 separated the spring communities from the fall and summer communities, and NMDS axis 2 separated the fall and summer communities. Overall the bacterial community structure in summer and autumn was relatively similar (Fig. 5A). The stress value of the fungal community NMDS ordination was 0.119, which similarly indicated that the plots were appropriate for representing community dissimilarities (Fig. 5B). As in the bacterial community analyses, the NMDS ordinations indicated a general segregation of fungal communities by season, albeit with less discrete clustering by season. NMDS axis 1 generally separated the spring rhizosphere communities from those in the summer and autumn, while axis 2 generally separated communities from the summer and autumn. Overall the fungal community structure in summer and autumn was similar (Fig. 5B).

PERMANOVA tests indicated significant differences among bacterial communities from the summer and autumn samples (F = 2.8234, p = 0.041), summer and spring samples (F = 6.2309, p = 0.001), and autumn and spring samples (F = 5.6672, p = 0.029). Likewise, significant differences in fungal communities were observed only from the autumn and spring samples (F = 1.6186, p = 0.024) while the summer fungal communities were not significantly different from the spring and autumn communities (p > 0.05) (Table 3).

No bacterial taxonomic groups were significantly associated with summer samples and the number of groups associated with spring samples were greater than for autumn samples (Fig. 6). In particular, Bacteroidetes, Firmicutes, and Verrucomicrobia were important taxa that distinguished spring bacterial communities. The Bacteroidia class of the Bacteroidetes group was responsible for their delineation of spring communities. Likewise, the Ruminococcaceae family of the Firmicutes phylum and Clostridiales order also discriminated spring communities, while the Clostridia, and Lactobacillaceae family of the Lactobacillales order also distinguished spring communities. In contrast, the Bacillaceae and Planococcaceae families of the Bacillales order differentiated autumn samples from others (Fig. 6).

Effects of soil physicochemical properties on elm rhizosphere microbiota

To evaluate the environmental factors influencing bacterial and fungal community structures of elm rhizospheres, redundancy discriminant analysis (RDA) was used to assess community-environment relationships. The RDA analyses indicated that different environmental factors exhibited different effects on bacterial and fungal communities, wherein the main influencing factors also varied among seasons (Fig. 7). The first axis of the bacterial biplot explained 30.4% of the variation in sample-environment relationships, while the first and second axes together explained 50.9%. MC was the most significantly correlated environmental variable to bacterial community composition, while AK was also significantly correlated. The environmental factors that primarily affected bacterial communities in the spring were MC and AK, while the factors that affected bacterial communities in the summer were primarily P, AP, and OM. Lastly, pH exhibited a minor correlation to bacterial community composition in autumn samples (Fig. 7A). For fungal community variation, the first axis explained 43.6% of the variation in sample-environment, relationships, while the first and second axes together explained 66.7% of the variation. AN was the environmental variable mostly explaining fungal community compositions, although P, AP, N, and ST were also influential. Lastly, pH was modestly influential towards community compositional variation (Fig. 7B).

Bacterial community functional prediction

A total of 41 level two KEGG Orthology (KO) groups were identified in the bacterial communities that were distributed across six metabolic pathways. Among these pathways, those involved in metabolism, genetic information processing, and environmental information processing were most prevalent, accounting for 52.15%, 15.7%, and 13.65% of the totals. At a finer level of resolution, the inferred relative abundances of gene families associated with membrane transport (11.41%), amino acid metabolism (11.14%), and carbohydrate metabolism (10.62%) were particularly high in the elm rhizosphere communities (Fig. 8). Nevertheless, multiple comparison analyses indicated that the predicted distribution of gene functions did not significantly differ among rhizosphere communities from different seasons (p > 0.05).

Fungal community functional analyses

Functional analysis of the fungal communities using the FUNGuild software program indicated that eight different ecological guilds were present among the fungal communities of the elm rhizospheres, in addition to unidentified guilds. Among these, saprotrophs and pathotrophs were most abundant, accounting for 32.1% and 11.1% of the total communities (Fig. 9). Multiple comparison analyses indicated that the distribution of fungal functional guilds did not significantly differ among rhizosphere communities from different seasons, as was observed for the bacterial communities (p > 0.05).