Hydrological and soil physiochemical variables determine the rhizospheric microbiota in subtropical lakeshore areas

Sample collection

A detailed overview of the studied lakes and water level features of the sampling area was given in our previous report (Zhang et al., 2018). The RMs were collected from lakeshore areas of three lakes: Lake Chaohu (LC), Lake Wuchang (LW), and Lake Dahuchi (LD) in summer (August 28–31, 2016; SU), autumn (November 17–20, 2016; AU), winter (February 18–20, 2017; WI), and spring (May 13–15, 2017; SP). One transect free of any artificial disturbance was established at each lake (Fig. 1). Along each transect, five sites numbered I, II, III, IV, and V were sequentially set perpendicular to the lakeshore from the mean annual lowest water level to its highest water level (Fig. 1). The elevation differences were equal between any two adjacent sites in each transect. Three rhizosphere soil samples were randomly collected from each site using a portable root-soil core sampler (4.0 cm diameter with a 17 cm depth), and these samples mixed to form a composite sample. The latter was named according to the season, lake, and location of sampling: for instance, “SPLWIII” refers to a sample obtained from site III of Lake Wuchang in spring. A total of 58 samples were collected from the three lakes spanning four seasons; the samples SULDV and WILCV were not collected due to logistical limitations in the field.

Measurement of physiochemical and hydrological variables

Rhizosphere soil pH and moisture were measured using the potentiometric method and the oven-drying method, respectively (Lu, 2000). Organic content (OC) was measured using the K₂Cr₂O₇ titration method, while the total nitrogen (TN) and total phosphorus (TP) contents were measured respectively with the Kjeldahl method and molybdenum blue colorimetry (Lu, 2000). Daily average water level data from 2012–2016 was used to calculate the fluctuation range of the water level (FRWL), submerged duration (SD), and average submerged depth (ASD); the data were obtained from the Jiangxi Poyang Lake National Nature Reserve and relevant hydrological website (http://yc.wswj.net/ahsxx/LOL/public/public.html). The FRWL was calculated as the difference between the highest and lowest water level for a given lake within a calendar year. The SD was calculated as the sum of days the site was under water, and the ASD was calculated according to the relative elevation of each sampling site and daily hydrological data.

DNA extraction and sequencing

Microbial DNA was extracted from samples using the SDS-based DNA extraction method, as done in earlier studies (Ni et al., 2010; Ni et al., 2012; Natarajan et al., 2016). The DNA integrity was checked by 1.0% agarose gel electrophoresis at 120 V for 30 min. The DNA concentration and purity were evaluated using a Nanodrop 2000 spectrophotometer (Thermo Scientific, USA) and then diluted to 1 ng/µl using sterile water. The V4 hypervariable region of prokaryotic 16S rRNA gene was amplified by using the prokaryotic specify primer set 515F and 806R with sample-specific barcode sequences, as done in earlier studies (Xiang et al., 2018). In brief, each 25-µl reaction mix contained 1 × PCR buffer, 0.25 U of Taq polymerase (Transgen Biotech, China), 0.2 mM of each dNTP, 1.0 µM of each primer, and 10 ng of microbial genomic DNA. The thermal cycling procedure was predenaturation at 94 °C for 10 min, followed by 30 cycles of 94 °C for 30 s, 56 °C for 30 s, and 72 °C for 30 s, with a final extension at 72 °C for 10 min. After this amplification, the PCR products were subjected to electrophoresis using a 2% agarose gel and quantified using a Nanodrop 2000 spectrophotometer (Thermo Scientific, USA). All amplicons were then pooled together with an equal molar amount from each sample (Huang et al., 2018) and purified using a DNA gel extraction kit (QIAGEN, Germany). Next, the pooled amplicons were sequenced using an Illumina HiSeq 2500 system at Beijing Novogene Technology Co., Ltd.

Data processing and analysis

Raw data were merged with tags by using FLASH v1.2.7 (Magoč & Steven, 2011), and divided among the samples according to the barcode sequences using QIIME v1.9.0 (Caporaso et al., 2010). After removing the barcode and primer sequences, the reads of low-quality sequences were detected and removed in QIIME v1.9.0, as recently described by Huang et al. (2018). Next, any chimeras present were sought and filtered out by UCHIME software with the “Gold” database (http://drive5.com/uchime/uchime_download.html). The sequences were classified into operational taxonomic units (OTUs) by setting a threshold of 97% identified sequence by using UPARSE v7.0.1001 software (Edgar, 2013); the highest frequency sequence in each OTU was then selected as the representative sequence of the OTU.

Phylogenetic information for each OTU was annotated, according to the representative sequence, by using the Mothur pipeline-referenced SILVA SSUrRNA database (Quast et al., 2013). A phylogenetic tree of the OTUs was constructed in MUSCLE v3.8.31 (Edgar, 2004). Four alpha-diversity indexes—observed OTU counts, Chao1 index, Shannon index, and Good’s coverage—were calculated using QIIME v1.9.0. Rarefaction curve, and rank-abundance curve were drawn using R v2.15.3. The weighted UniFrac distances were calculated using QIIME v1.9.0.

All the bacterial sequences have been deposited in the NCBI SRA database under accession number SRP161734.

Statistical analysis

Results for each variable are presented as the mean ± standard error for each group (Huang et al., 2018). The non-parametric Kruskal-Wallis rank sum test and post-hoc tests were used to identify significantly different taxa among different groups with STAMP software (Parks et al., 2014). Correspondence analysis (CA), canonical correspondence analysis (CCA), Non-metric multidimensional scaling (NMDS), and non-parametric multivariate analysis of variance (MANOVA) were conducted using the “vegan” package in the R platform. A heatmap profile was drawn using HemI software. Separate Wilcoxon tests were used to compare the α-diversity indexes between different groups, by using the agricolae package in R.

Environmental variables and taxonomic composition of the RMs

Both OC and TN differed significantly between spring and winter (Figs. 2A and 2C, and Table S1), as did soil pH, but it also was significantly different between summer and autumn, and likewise between summer and winter (Fig. 2B, and Table S1). Soil moisture, pH, elevation, and FRWL were significantly different among the three lakes. OC and TN were significantly different between LC and LD, and between LD and LW, while TP was significantly different between LC and LD, and between LC and LW (Figs. 2D–2J, and Table S1). The SD and ASD were significantly different among sampling sites; specifically, the OC of site I differed from all other sites; TN differed between sites I and V; moisture of site V differed from all other sites except I (Figs. 2K–2O, and Table S1).

After removing low quality sequences, a total of 2 488 433 (42 904.02 ± 1201.19 per sample) high-quality sequences were obtained for the 58 samples. To eliminate the influence of sequencing depth upon our results, all samples were randomly resampled to 24 038 sequences for further analysis, which was the lowest number of sequences per sample. Besides a few unclassified sequences (1.11 ± 0.14%), other sequences could be classified into 14 archaeal and 58 bacterial phyla. However, only four archaeal—i.e., Euryarchaeota (2.41 ± 0.47%), and Thaumarchaeota (2.31 ± 0.55%), DHVEG-6 (1.04 ± 0.19%), Bathyarchaeota (0.59 ± 0.15%)—and 21 bacterial phyla—Proteobacteria (32.64 ± 1.12%), Acidobacteria (23.79 ± 1.76%), Nitrospirae (8.16 ± 0.61%), Firmicutes (4.76 ± 0.54%), Chloroflexi (4.51 ± 0.33%), Bacteroidetes (2.98 ± 0.34%), Gemmatimonadetes (2.60 ± 0.22%), Verrucomicrobia (2.39 ± 0.27%), Actinobacteria (1.84 ± 0.19%), Ignavibacteriae (1.52 ± 0.11%), Latescibacteria (1.37 ± 0.15%), Spirochaetes (0.83 ± 0.09%), Aminicenantes (0.75 ± 0.20%), Planctomycetes (0.66 ± 0.15%), Zixibacteria (0.48 ± 0.05%), Parcubacteria (0.42 ± 0.06%), Cyanobacteria (0.41 ± 0.07%), Omnitrophica (0.20 ± 0.07%), GAL15 (0.17 ± 0.06%), AC1 (0.12 ± 0.03%), and Chlamydiae (0.08 ± 0.04%) —were found to be dominant, in that their relative abundances exceeded 1% in at least one of our samples (Fig. 3A). In Proteobacteria, Deltaproteobacteria was the most abundant class, followed by Betaproteobacteria, and Alphaproteobacteria (Table S2).

Comparing the RMs among different lakes, seasons, and sites

Within these four archaeal and 21 bacterial phyla, the relative abundances of Acidobacteria, Actinobacteria, Bacteroidetes, Bathyarchaeota, Gemmatimonadetes, and Proteobacteria were significantly different among the three lakes. More than half of these dominant phyla also differed significantly among seasons in their relative abundance, while the relative abundance of DHVEG-6, Ignavibacteriae, Nitrospirae, Spirochaetes, and Zixibacteria were significant differences among the sampling sites (Fig. 3B).

Almost all the sequences (98.89 ± 0.14%) were classified into phyla, for a total of 21 361 OTUs detected. The rarefaction curve showed that most samples reached the platform, which implied that the sequences could represent the microbiota structures of the RMs (Fig. S1A). The rank-abundance curve showed that the OTUs in the microbiota were very uneven (Fig. S1B). Of the 21 361 OTUs, only 182 of them dominated the RMs (i.e., relative abundance >1% in at least one sample; Fig. S2), harboring 30.87 ± 1.53% of all the high-quality sequences, which were consistent with the result of the rank abundance curve (Fig. S1B). The microbiota first tended to cluster according to lakes, and cluster according to seasons within the same lake. However, no evidence was found to indicate the microbiota clustered according to sampling sites of the lakeshore areas (Fig. S2).

Spatially, there were 118 dominant OTUs for which we detected a significant difference in their relative abundance among the three lakes. The relative abundances of the dominant OTUs classified into Nitrospira cf. moscoviensis SBR1015, Gaiella sp., Tenderia sp., Steroidobacter sp., Steroidobacter sp. WWH78, Methylomirabilis sp., Sphingomonas sp., Methylotenera sp., Sulfuricurvum sp., Methanosaeta sp., Nitrosoarchaeum sp., Clostridium beijerinckii, Thiobacillus sp., Eubacterium sp., Sulfurifustis sp., family Gemmatimonadaceae, Gallionellaceae, Rhodospirillaceae, Nitrospiraceae, Desulfurellaceae, MIZ17, Sh765B-TzT-35, 0319-6A21, order Xanthomonadales, Holophagae, NB1-j, class Acidobacteria, and SAGMCG-1 were all significantly higher in LC than in the other two lakes. The relative abundances of the dominant OTUs classified into Koribacter sp., Nitrotoga sp., Acidibacter sp., Solibacter sp., Terracidiphilus sp., family Acidobacteriaceae, DA111, FW13, ASC21, Gallionellaceae, order Sva0485, Holophagae, class Acidobacteria, JG37-AG-4, and SAGMCG-1 were all significantly higher in LD than in the other two lakes. Finally, the relative abundances of the dominant OTUs classified into Koribacter sp., Rhodanobacter sp., Geobacter sp., Methanoperedens sp., Nitrosotalea sp., family Nitrospiraceae, Nitrosomonadaceae, Acetobacteraceae, Gemmatimonadaceae, order Sva0485, and MSBL5 were all significantly higher in LW than in the other two lakes.

Temporally, 112 of the 182 dominant OTUs detected showed significant differences in their relative abundance among the four seasons. In the summer, the dominant OTUs classified into Clostridium beijerinckii, Thiobacillus sp., family 0319-6A21, FW13, Acidobacteriaceae, Sh765B-TzT-35, order Holophagae, Sva0485, and Xanthomonadales were significantly increased, while those classified into Nitrosoarchaeum sp., and family Acidobacteriaceae significantly decreased. In autumn, the dominant OTUs classified into Anaerostipes hadrus, Acetobacter pasteurianus, Acidibacter sp., Eubacterium sp., Faecalibacterium sp., Lactobacillus vini, Lactobacillus sp., Roseburia inulinivorans, Ruminococcus bicirculans, Ruminococcus sp., Subdoligranulum sp., Serratia marcescens, Methylomirabilis sp., Koribacter sp., Nitrosotalea sp., Methanoperedens sp., family Nistrospiraceae, Gallionellaceae, and FW13 were significantly increased, while those classified into Sideroxydans sp., Haliangium sp., Sulfuricurvum sp., and family Nitrosomonadaceae, order NBi-j, class Acidobacteria, and ML635J-21 significantly decreased. In winter, the dominant OTUs classified into Carnobacterium maltaromaticum, Telluria mixta, Enterobacter sp., Methylotenera sp., Solibacter sp., family Methylococcaceae, Gemmatimonadaceae, order MSBL5, Holophagae, NB1-j, and Holophagae were significantly increased, while those classified into family Nitrospiraceae, MIZ17, 0319-6A21, order TRA3-20, and class Acidobacteria significantly decreased. Finally, in spring, the dominant OTUs classified into Clostridium sp. ND2, Sideroxydans sp., Terracidiphilus sp., Geobacter sp., Arthromitus sp., Gallionella sp., Sulfuricurvum sp., Sphingomonas sp., Nitrotoga sp., family Nitrospiraceae, NIZ17, Acidobacteriaceae, class Acidobacteria, ML635J-21, and SAGMCG-1 were significantly increased, while those classified into Koribacter sp., and class JG37-AG-4 significantly decreased.

Among sampling sites, we detected 43 dominant OTUs with significant differences across this elevation gradient. The relative abundances of the most significantly different dominant OTUs gradually changed with elevation: those classified into Methanoperedens sp., family ASC21, Gallionellaceae, Nitrospiraceae, order 43F-1404R, Sva0485, and MSBL5 were diminished from sampling site I through V, whereas those classified into Solibacter sp., Sphingomonas sp., Terracidiphilus sp., and family Acidobacteriaceae gradually increased.

More OTUs were detected from the lakeshore RMs in spring than in the other seasons (Kruskal-Wallis test, χ² = 11.935, p = 0.008; Fig. S3A). The OTU counts at LC were more diverse than those occurring in the other two lakes (Kruskal-Wallis test, χ² = 13.611, p = 0.001; Fig. S3B). Although there was a trend for the OTU counts to decline from sampling site I to V, no significant difference was detected among the sampling sites (Kruskal-Wallis test, χ² = 4.878, p = 0.30; Fig. S3C). The Shannon indexes were not significantly different among the four seasons (Kruskal-Wallis test, χ² = 7.396, p = 0.060; Fig. S3D). However, the Shannon indexes of the RMs at LC exceeded those of the other lakes (Kruskal-Wallis test, χ² = 21.112, p < 0.001; Fig. S3E). No significant difference was detected among the five sampling sites in the Shannon index of the lakeshore RMs (Kruskal-Wallis test, χ² = 2.619, p = 0.624; Fig. S3F). Chao 1 index of the lakeshore RMs in SP was significantly higher than other seasons (Fig. S3G), while the Goods’ coverage of the lakeshore RMs in SP was significantly lower than other seasons (Fig. S3J). No significant difference was detected among the lakes and sampling sites in the Chao1 index and Goods’ coverage (Figs. S3H, S3I, S3K, and S3L). Because the sampling sites were closely related to ASD, we also analyzed the correlation between the α diversity indexes and the ASD. Our results showed that except Shannon index significantly positively correlated with ASD (F-test, F = 4.35, p = 0.04; Fig. S4), OTU count, Chao1 index, and Goods’ coverage did not significantly correlate with ASD (F-test, p > 0.05; Fig. S4).

Although the CA based on all OTUs did not absolutely distinguish the RMs from different lakes, seasons, or sampling sites (Figs. 4A–4C), the CCA with Monte Carlo testing and MANOVA revealed significant differences in the RMs among lakes (MANOVA, F = 5.99, p < 0.01), seasons (MANOVA, F = 3.51, p < 0.01), and sampling sites (MANOVA, F = 4.03, p < 0.01) (Fig. 4D). NMDS result also showed that the RMs had trends to distinguish according to lakes and sampling sites (Fig. S5).

As dispersal limitation caused by geographic distance is one of the major mechanisms that maintain β-diversity of microbial communities (Eisenlord, Zak & Upchurch, 2012; Ni et al., 2014; Cao et al., 2016; Shirani & Hellweger, 2017), we analyzed the variation of the weighted UniFrac distances between the sampling sites with the geographic distances. Our results showed that the weighted UniFrac distances significantly increased with the increases of the geographic distances (Fig. 5). However, the R² of the regression equations were very low (<0.20; Fig. 5).

Influence of environmental factors on the taxonomic composition of RMs

CCA was used to analyze which environmental factors best explained the microbiota structure in lakeshore rhizospheres. Since SD and ASD, pH and elevation, and OC and TN were significantly correlated with each other (the absolute value of Pearson correlation R-values >0.9 and p-values < 0.001; Fig. S6), the elevation, OC, and ASD variables were deleted before conducting the CCA. When performed with the Monte Carlo test, it showed that FRWL and pH (or elevation) were the most important factors, followed by TN (or OC), moisture and TP, for significantly influencing the lakeshore RMs (Fig. 6A) and their dominant microbiota (Fig. 6B).

To analyze which dominant OTUs were influenced by hydrological and soil physiochemical variables, linear regressions were used. The FRWL was a positive predictor of the dominant OTUs in the genus Candidatus Koribacter, genus Candidatus Solibacter, family Acidobacteriaceae, and class Acidobacteria, but it was negatively related to the dominant OTUs in the genus Candidatus Methanoperedens and the phylum Bathyarchaeota (Table S3). A greater soil pH negatively affected the dominant OTUs in the genus Candidatus Koribacter, genus Candidatus Solibacter, genus Sphingomonas, genus RB41, family Rhodospirillaceae, family Acidobacteriaceae, family FW13, order Sva0485, and class Acidobacteria, but positively influenced those OTUs in the genus Geobacter, family Nitrospiraceae, order 43F-1404R, order Sva0485, and phylum Bathyarchaeota. Both soil TN and TP largely had a positive influence on the dominant OTUs in the family 0319-6A21 and family Rhodospirillaceae, respectively (Table S3).