A comparative study of bacterial diversity based on effects of three different shade shed types in the rhizosphere of Panax quiquefolium L.

General situation of this experiment

The experiment was performed in the Weihai farmland ginseng experimental area of Shandong Province (122°27′E, 36°96′N). Weihai is among the leading American ginseng production areas of China. As a shade-loving plant, American ginseng must be cultivated under shading conditions. Therefore, farmers began to use shade shed to ensure the yield and quality of American ginseng. This study adopted traditional (CTP), high flag (PP), and high arch shade sheds (GP) (Fig. 1). The traditional shed is a flat, wooden structure of 1.6 m, height. The high flat shed is a flat alloy structure, 2.2 m in height, while the high arch shed is an alloy structure measuring 2.6 m in height. All three sheds have double-layer shading nets.

Basal fertilizer (∼3–4 t/mu) was applied once before the experimentation, conducted between 2016 and 2018. The shading shed was a completely random block with three replicates.

Soil collection

At the end of October 2018 when 3-year-old ginseng roots were manually harvested, 0–10 cm of surface soil was randomly sampled. Six soil samples from one type of shed were combined to form one sample. Each shed type had three composite samples, namely CTP1, CTP2, CTP3, PP1, PP2, PP3, GP1, GP2 and GP3. One part of each sample was stored at −80 °C for subsequent DNA extraction. The other portion was homogenized and sieved with a 2-mm mesh sieve to remove large particles, air-dried and used for chemical analysis. Soil pH was determined at 1:2 soil: water ratio using a combination glass electrode (FE20K, Mettler-Toledo, Swiss). Soil total organic matter (TOM) was assayed following Walkley-Black method (Walkley & Black, 1934). Hydrolyzable nitrogen (HN) and available phosphorus (AP) were determined by chemical analytical methods, as described by Liu et al. (2014) and Murphy & Riley (1962). The chemical characteristics of these soil samples are shown in Table S1. The yield of American ginseng, disease incidence, type and severity are shown in Table S2. The disease severity index was calculated using the method described by Chiang et al. (2017).

DNA extraction, PCR amplification and sequencing

Total DNA was extracted using the Z.N.ATM Mag-Bind Soil DNA Kit (Omega Bio-Tek, GA, USA). DNA was separated and visually tested for quality on a 0.8% (w/v) agarose gel electrophoresis. The DNA concentration was measured with a NanoDrop™ 2000 spectrophotometer (Thermo Fisher Scientific, MA, USA). DNA extracts were preserved at −20 °C for PCR.

The bacterial amplicon library was obtained using the 338F (5′-ACTCCTACGGGAGGCA GCA-3′) (Wu et al., 2016) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) (Wu et al., 2016) primers, which target the V3–V4 region of the 16S rRNA gene. Personalbio Technology (Shanghai, China) generated 300 bp paired-end reads from the qualified libraries using an Illumina MiSeq platform (Illumina, CA, USA) following the manufacturer’s instructions. The raw reads were deposited to the NCBI Sequence Read Archive under accession number PRJNA662686.

Sequencing data analysis, OTU production and annotation

Split reads were sorted into each sample by their unique barcodes using QIIME V.1.8.0 (Bokulich et al., 2013; Caporaso et al., 2010). The sequencing barcodes and primer sequences were cut off from original paired-end reads using the FLASH V1.2.7 (Magoč & Salzberg, 2011). The USEARCH software (v5.2.236, http://www.drive5.com/usearch/) was used to detect and remove chimera sequences.

The effective tags with ≥97% similarity were assigned to the same operational taxonomic units (OTUs) using the Uparse software (v8.1.1861). The sequence with the highest occurrence frequency in each OTU was considered the representative sequence for further annotation (Edgar, 2013). Each representative sequence was annotated using the UCLUST method (Edgar, 2010) on the Greengenes database (Release 13.8, http://greengenes.secondgenome.com/) (DeSantis et al., 2006). The annotation levels were: kingdom, phylum, class, order, family, genus and species to determine the community composition of each sample. Rare OTUs (singletons to tripletons), potentially originating from artificial sequences, were removed. The read counts were normalized considering the smallest read number per sample.

Statistical analysis

All statistical analyses were calculated using the QIIME software (Version 1.8.0) and displayed using the R software (Version 2.15.3). All data were analyzed using a one-way analysis of variance (ANOVA) (P value < 0.05) with the SPSS19.0 software (SPSS Inc., Chicago, USA). The differences between groups were compared using Fisher’s least significant difference test.

Alpha diversity was applied to analyze the complexity of species diversity for each community through four indices, including abundance-based coverage (ACE) estimator, Chao 1, Shannon diversity and Inverse Simpson indices. The ACE estimator and Chao 1 index determined the community richness, while the Shannon and Inverse Simpson indices estimated the community diversity. All samples were randomly resampled at the lowest sequencing depth level of 90% to correct the diversity difference caused by sequencing depth. Then, QIIME software is used to calculate alpha diversity.

The principal component analysis (PCA) determined the beta diversity to reduce the dimension of the original variables. Hierarchical clustering was performed using the Unweighted Pair-Group Method with Arithmetic Means (UPGMA) clustering to interpret the distance matrix using average linkage.

The Spearman’s rank correlation analysis performed using the Mothur software (Schloss et al., 2009). The Cytoscape tool (http://www.cytoscape.org/) displayed the first 50 correlated groups. The PICRUSt tool (http://huttenhower.sph.harvard.edu/galaxy/tool_runner?tool_id=PICRUSt_normalize) determined the enriched pathways following the Kyoto Encyclopedia of Genes and Genomes (KEGG) of the microbiota. The variation analysis was performed following the functional abundance of the samples.

The SPSS 19.0 software determined the Spearman’s rank correlation analysis between the microbial community composition and soil physicochemical parameters. Meanwhile, the Redundancy analysis (RDA) was performed using the Canoco 5.0 tool to identify the relationships between the bacterial community structures and soil physicochemical parameters.

Sequencing statistic

A total of 660,115 raw tags were generated from 16S rRNA genes of the V3–V4 region. After trimming and filtering, 430,843 high-quality reads were obtained, with a 415 bp average reading length. The OTU clustering (97% similarity) generated 33,858 OTUs from the nine samples, with 99.97% classified under phylum. The rarefaction curves indicated that the sampling sufficiently captured the overall composition of each sample (Fig. S3).

Taxonomic composition of bacterial assemblages in rhizosphere soil

The taxonomy of 16S rRNA gene amplicons and the relative abundance of the top twenty taxa estimated the bacterial assemblages and variability in American ginseng rhizosphere soils from three different shade sheds (Fig. S4). At the phyla level, Proteobacteria was the most prominent taxon (30.5%) of the total bacterial community, followed by Chloroflexi (17.2%), Acidobacteria (15.4%), Actinobacteria (15.3%), and Gemmatimonadetes (5.8%) (Fig. S4A). The dominant classes with >1% relative abundance across all samples included Alphaproteobacteria, Acidobacteriia, Gammaproteobacteria, KD4-96, Gemmatimonadetes, Bacteroidia, Actinobacteria, Thermoleophilia, Deltaproteobacteria, Subgroup_6, Saccharimonadia, Ktedonobacteria, Anaerolineae, Acidimicrobiia, Clostridia, AD3, Verrucomicrobiae, Planctomycetacia, Blastocatellia (Subgroup_4), and Parcubacteria (Fig. S4B). Bacterial variation increased from phylum to genus, following analysis of the ten most-abundant classes at the order, family and genus levels (Figs. S4C–S4E).

Alpha and beta diversity analysis of samples

The ACE, Chao 1, Shannon diversity, and Inverse Simpson indices quantified the alpha diversity and per sample bacterial diversities. The ACE and Chao 1 in PP were significantly lower (P < 0.05) than GP and CTP (Fig. 2), demonstrating that shading affected the rhizosphere bacterial richness. The Shannon diversity and Inverse Simpson indices showed no significant differences in bacterial diversity among the three shade sheds.

At the OTU level, PC1 explained 60.77% and PC2 13.07% of the total bacterial variation, clearly distinguishing the community structure of soil bacteria under the three shade sheds (Fig. 3A). All the samples collected from the rhizosphere were distinguishable based on the UPGMA-generated cluster trees (Fig. 3B).

Variations in the bacterial community of American ginseng rhizosphere soil under three different shade sheds

The least discriminant analysis (LDA) effect size taxonomic cladogram compared the phylum to the genus in the three sample communities to expose the bacterial variation in the American ginseng rhizosphere soil under three different shade sheds. The relative abundances of bacterial phyla were different across the three shade shed rhizosphere soils (Fig. 4). The relative abundances of Proteobacteria, Nitrospirae and WS2 in GP were significantly higher than PP and CTP (LSD, p < 0.05). Dependentiae, Elusimcrobia, FCPU426 and WPS-2 were significantly abundant in the pp than GP and CTP (LSD, p < 0.05). At a finer taxonomic level, the three shade sheds were significantly different in some bacterial genera. The abundance values of the 50 most abundant genera from the three shade sheds are shown in Fig. S5.

Correlation and KEGG functional analysis in bacterial communities

The relationships between abundance values of the various bacterial genera were determined following the abundance changes of different genera in the soil samples. The genera relationship network indicated a complex functional collaboration within the microbiota (Fig. S6). The KEGG functional enrichment analysis of bacterial microbiota from rhizosphere soils under the three different shade sheds determined their bacterial functional diversities. The abundance of the top 50 functional groups was analyzed and displayed on a heat-map (Fig. 5). A Venn diagram (Fig. 6) shows the number of OTUs of the common functional soil bacteria in the three shade sheds, indicating the functional divergence of bacterial microbiota in the rhizosphere of American ginseng roots.

Relationship between microbial communities and soil chemical properties

The Spearman’s rank correlation coefficient determined the relationship between soil microbial diversity and soil chemical properties (Fig. 7). The Chao 1 and Shannon diversity indices positively and significantly correlated with the soil pH (P < 0.05). However, Chao 1 and ACE indices negatively correlated with soil TOM and AP contents (P < 0.01), while the Inverse Simpson index negatively correlated with HN (P < 0.05).

RDA analysis was performed to observe the potential correlations between the community structures at genus levels and soil chemical properties in the GP, PP and CTP. The RDA analysis results showed that the first two RDA axes (RDA1 and RDA2) account for 60.30% and 12.67% of the variance, respectively (Fig. 8). Soil pH (F = 9.2, P = 0.016) was the most important contributor to bacterial community variation (Table S7).