Spatial and temporal dynamics of the bacterial community under experimental warming in field-grown wheat

Experimental site and design

The soil warming experiment was established at the Luancheng Agro-Ecosystem Experimental Station of the Chinese Academy of Sciences on the North China Plain, Hebei, China (37°53′N, 114°41′E) in 2008. The field was cultivated with local winter wheat (Triticum aestivum L.) cultivar ‘Shixin 828’. The soil was classified as a sandy loam with a pH of 8.1, 15.1 g kg⁻¹ organic matter, and 1.1 g kg⁻¹ total N in the 0–20 cm soil layer (Liu et al., 2016). The control and warming treatments were set up in a randomized block design, each replicated three times with an individual plot size of 4 m × 4 m. The warmed plots were heated with three infrared heaters (1000 W, size of 2 m × 0.02 m) that were installed at the center of the plot 2 m above the ground, which were distributed equally in the 2 m long area, and the radiation area was 2 m × 2 m. In the control plots, “dummy” heaters were installed with no power to imitate shading effects (Fig. S1). The daily average temperature of the topsoil in the warming plots was 1.5 °C higher than that of the control plots. The nitrogen fertilizer was urea, half of which was applied before sowing in October and the other half in April of the following year. All phosphate fertilizers were applied at 65 P₂O₅ kg hm⁻².

Plant and soil sampling

Plant root-zone (root, rhizosphere and bulk soil) samples were collected in November (Feekes growth stage 2–3, tillering stage), March(Feekes stage 6–7, jointing stage) and May (Feekes stage 11, ripening stage) during the wheat growing season. At each growth stage, root samples and rhizosphere soil of the control and warming treatments were randomly taken with three replicates using the method described previously (Chen et al., 2019). After gently shaking the roots to remove loosely bound soil clumps, the rhizosphere soil was carefully brushed out of the roots (Clemensson-Lindell & Persson, 1992). The roots were washed with sterile distilled water and used for endosphere bacterial community analyses. The sample collection method did not discriminate between microbial communities at the root surface and in the endosphere; therefore, we considered the root fraction as the ‘root microbiome’ (Hu et al., 2018). To sample bulk soil, three soil cores (0–20 cm soil depth) were randomly taken from each plot to form a composite sample. In total, three replicate plots from the warming treatment and three from the control were sampled at each growth stage. Roots, rhizosphere and bulk soil samples were stored at −80 °C before DNA extraction and at 4 °C for measurement of dissolved organic carbon (DOC) and enzyme activities.

Soil properties and dehydrogenase activities

DOC is defined as the dissolved part of organic carbon in the soil. DOC in the rhizosphere was extracted with sterile distilled water and measured by an elemental analyzer (vario TOC; Elementar Analysensysteme GmbH, Langenselbold, Germany). Soil pH was measured using a pH electrode (1:5, soil:water) (Rayment & Higginson, 1992), soil temperature and soil moisture were monitored continuously using T-type thermal couples and time domain reflectometry (TDR 100 system, USA), respectively, and data were recorded every hour using a data logger (CR10X, Campbell, USA) in control and warmed treatment plots.

The dehydrogenase enzyme activity is considered a marker of microbial activities in the soil and was assessed for the rhizosphere soils in the current study. Triplicate samples collected at the tillering, jointing and ripening stages from both the control and warming treatments were used for enzyme activity determination. The soil was passed through a two mm sieve and stored at 4 °C, and activity was measured within one week of sampling. In brief, 6 g of fresh soil and 60 mg of CaCO₃ were incubated with 1 ml of 3% triphenyl tetrazolium chloride (TTC) and 2.5 ml of deionized water at 37 °C for 24 h. Triphenylformazan (TPF), a product from the reduction of TTC, was extracted with methanol in a 100 ml volumetric flask, and the color intensity was measured at a wavelength of 485 nm (Tabatabai, 1994).

DNA extraction, PCR amplification and sequencing

The DNA was extracted from 0.5 g of fresh root powder that was obtained by grinding with liquid nitrogen. Total DNA of the rhizosphere and bulk soil was extracted using an E.Z.N.A. Soil DNA Kit (Omega Biotek, Inc., Norcross, GA, USA) following the manufacturer’s instructions. The concentration and quality of extracted DNA were determined using a NanoDrop spectrophotometer (NanoDrop-2000c Technologies, Inc., Wilmington, DE, USA), and extracted DNA was stored at −20 °C until further use.

The bacterial 16S rRNA gene (V3–V4 region; approximately 460 bp) was amplified with primers 341F (5′-CCTACGGGNGGCWGCAG-3′) and 785R (5′-GACTACHVGGGTATCTAATCC-3′) (Yasir et al., 2015). Overhanging bases were added to the primers to connect the Illumina sequencing adapters and dual-index barcodes in a second round of PCR. Each PCR was performed in a 25 µl mixture containing 12.5 µl of PCR Premix Ex Taq™ (Takara Biotech, Dalian, China), 1 µl of each primer (10 µM), and 1 µl of DNA template (approximately 20 ng DNA). The PCR conditions were as follows: 95 °C for 3 min; 25 cycles of 30 s at 95 °C, 30 s at 55 °C and 30 s at 72 °C, with a final extension at 72 °C for 10 min. The PCR products were visually examined on agarose gels and then purified with AMPure XP beads (Beckman Coulter, Inc., Brea, CA, USA) following the manufacturer’s protocol. Subsequently, eight-cycle PCR was carried out to add Illumina sequencing adapters and dual-index barcodes to each sample, and then the PCR product was purified using AMPure beads. The libraries were then normalized according to the Nextera XT (Illumina) protocol, and samples were sequenced on a MiSeq PE300 platform (GENEWIZ, Suzhou, China).

Sequence processing and analysis

The raw sequences were processed mainly with QIIME2 (2020.11) (Bolyen et al., 2019). The 16S rRNA gene sequences were quality filtered and denoised using DADA2, followed by the creation of amplicon sequence variants (ASVs) using the Deblur tool (Callahan et al., 2016). The taxonomic identities of the ASVs were obtained using the QIIME2 feature-classifier plugin (sklearn method) against the SILVA v.138 database (Quast et al., 2013). The alpha and beta diversity of the microbial communities were calculated within QIIME2 based on the standardized ASV table. PICRUSt2 (https://github.com/picrust/picrust2) was utilized to predict the functional potential of the microbial community (Douglas et al., 2020). Sequencing data were deposited into the European Nucleotide Archive under accession number PRJEB37653.

Statistical analyses

All statistical analyses were carried out with SPSS 20.0 (IBM, Chicago, USA) and R v4.0.3 (Team RC, 2014). Student’s t test was used to compare the means of the control treatment and warming treatment at the P < 0.05 level. The R packages “ggplot2” and “pheatmap” were used to draw the point plots, bar plots, and heatmaps of bacterial diversity and community composition. Principal component analysis (PCA) was performed on the bacterial community at the genus level using the packages “vegan” and the results were visualized with “ggplot2” to determine the effects of warming and wheat developmental stage on the community structure in the root, rhizosphere and bulk soil.

Soil characteristics and microbial activity

Warming was simulated in the wheat field plots by using infrared heaters. The soil moisture was significantly lower in the warmed plots than in the control plots (Fig. S2). The DOC increased under warming treatments (Fig. 1A). DOC increased with wheat development and was found to be higher at the ripening stage. Dehydrogenase activity (i.e., microbial activity) was higher (P < 0.05) under warming at all developmental stages (Fig. 1B).

Bacterial diversity and richness

High-throughput sequencing was carried out on the root, rhizosphere and bulk soil samples collected at the tillering, jointing and ripening stages of wheat grown under control and warming conditions. The effects of warming, root proximity and plant growth stages on bacterial diversity (Shannon and Simpson) and richness (Chao1) indexes were determined (Fig. 2). In the root compartment, bacterial diversity(Shannon index) generally decreased in the warmed plots compared to the control plots, and the Shannon and Simpson indexes were significantly higher at the ripening stage than at the tillering and jointing stages. In contrast, bacterial diversity and richness increased with warming, and the Shannon and Chao1 indexes were generally higher at the tillering and ripening stages in the rhizosphere soil. Alpha diversity also decreased in the bulk soil at the early wheat developmental stage but was not statistically significant.

Soil warming along with wheat development stage considerably influenced the root bacterial community when compared to the rhizosphere and bulk soil (Fig. 3, Fig. S3). The PCA showed greater separation for communities between the control and warmed plots in the root and rhizosphere compared with the bulk soil. In roots, PC1 and PC2 explained 50.0% and 21.3% of the variability in the bacterial community of samples from all wheat development stages. The separation of the bacterial community between the control and warmed plots decreased with the wheat development stage in root samples. In the rhizosphere, PC1 and PC2 accounted for 51.2% and 23.7% of the variance in the data, respectively. In bulk soil, there was no significant separation between treatments and the development stage.

Bacterial community responses to warming

The bacterial community responded differently to warming at spatial (root, rhizosphere and bulk soil) and temporal (tillering, jointing and ripening stages) scales. For example, warming considerably affected the bacterial community composition in roots compared with the rhizosphere and bulk soil at an early wheat development stage (Figs. 3 and 4A). Proteobacteria (sum of the Alphaproteobacteria, Betaproteobacteria, Deltaproteobacteria and Gammaproteobacteria) was the most abundant phylum across the control and warming treatments and was relatively higher in the roots (47–55%) than in the rhizosphere (38–43%) and bulk soils (24–29%). In the root, rhizosphere and bulk soil compartments, warming increased the abundance of Alphaproteobacteria at the tillering stage and decreased Betaproteobacteria at the tillering and jointing stages. Cyanobacteria were substantially more abundant in the roots (6–26%) than in the rhizosphere (0.2–0.6%) and bulk soils (0.3–1.9%), while Bacteroidetes was more abundant in the rhizosphere (13–19%) than in the roots (3.3–11%) and bulk soils(3.3–4.1%). The relative abundances of Acidobacteria and Planctomyces increased considerably with distance from the roots. The relative abundance in response to warming was calculated (Fig. 4B). The relative abundances of Actinobacteria increased in response to warming and responded more significantly in the rhizosphere and bulk soil. In the warming treatments, Acidobacteria and Alphaproteobacteria increased at the tillering stage, while Gammaproteobacteria increased in the root and rhizosphere at the tillering and jointing stages (Fig. 4B).

In the root compartment, the relative abundance of the order Rickettsiales was considerably higher than that in the rhizosphere and bulk soil compartments and decreased markedly from the tillering to ripening stages (Fig. S4). An increasing trend was observed for Rhizobiales from tillering to ripening in root and rhizosphere samples, while the abundance did not change significantly under warming. In the roots and rhizosphere at the tillering, jointing and ripening stages, the relative abundances of Actinomycetales, Bacillales and Pseudomonadales increased in response to warming, while the relative abundances of Sphingobacteriales and Burkholderiales decreased at the tillering and jointing stages in response to warming. Moreover, the abundance of Actinomycetales was higher in roots, while Sphingobacteriales and Burkholderiales were more abundant in the rhizosphere. In the roots and rhizosphere, the relative abundance of Rhizobiales increased with wheat development and was 2-3-fold higher at the ripening stage than at the tillering and jointing stages. The abundance of Rickettsiales was significantly higher in the roots (12–28%) than in the rhizosphere (0.2−0.4%) and bulk soil (0.06−0.31%).

Warming considerably influenced the bacterial genera at all wheat growth stages, and we analyzed the top 22 dominant bacterial genera (average relative abundance greater than 0.2%) (Fig. 5). In the roots and rhizosphere, the relative abundances of Pseudomonas, Promicromonospora, Saccharothrix, Bacillus and Arthrobacter increased after warming at the tillering and jointing stages. Moreover, the abundance of Pseudomonas was dramatically higher in the root and rhizosphere soil than in the bulk soil, and the abundances of Devosia and Streptomyces decreased with distance from the roots.

LEfSe analysis was performed to further identify differential species between the warming and control treatments in both root and rhizosphere bacterial communities. The differential species that met the linear discriminant analysis (LDA) significance threshold greater than 2.0 are shown in Fig. 6. The results confirmed that in the roots and rhizosphere, the phylum Actinobacteria, the order actinomycetes and the genera Saccharothrix, Arthrobacter, Promicromonospora, Glycomyces and Cellulosimicrobium increased significantly under warming conditions, and the relative abundance of the order Pseudomonadales and genus Pseudomonas also increased after warming.

Functional prediction of bacterial communities

Differences in the function of the bacterial communities in the warming and control treatments of different root zone compartments and growth stages were assessed using PICRUSt2 with the Kyoto Encyclopedia of Genes and Genomics (KEGG) database. Six types of biological metabolic categories at KEGG level 1 were obtained (Fig. S5) with metabolism as the primary category (46.7%–50.0%). The functional profiles of energy metabolism and carbohydrate metabolism at KEGG level 3 were then predicted and plotted as a heatmap for comparison (Fig. S6). We found that for the carbohydrate metabolism, warming elevated the functional categories including galactose metabolism, ascorbate and aldarate metabolism, pentose and glucuronate interconversions, pentose phosphate pathway, and pyruvate metabolism in the rhizosphere.