The rhizosphere and root selections intensify fungi-bacteria interaction in abiotic stress-resistant plants

Plant growth and sampling

Eleven plant species (Table S1), including Casuarina equisetifolia, Calophyllum inophyllum, Fagraea ceilanica, Guettarda speciosa, Heritiera littoralis, Hernandia sonora, Melaleuca bracteata, Pongamia pinnata, Portulaca pilosa, Ruellia brittoniana, and Scaevola sericea, natively grow and thrive under the environment of high soil salinity and alkalinity, high temperature, strong light, and frequent drought in Hainan, China (Tong et al., 2013; Liu et al., 2015; Ren et al., 2017; Li et al., 2018). Their seeds were collected in Hainan, and sown in pots (24 × 26 cm, diameter × height) filled with field soils near a greenhouse of Wenchang Lingye Gardening Limited Company, Hainan, China. For each plant species, 24 seeds were selected to germinate, and four seeds were sown in a pot. After germination, only one healthy seedling was left in the pot, and kept growing in the same greenhouse from February to August, 2017. All the plants were harvested on August 14, 2017. The soils collected before planting were designated as the bulk soils. The plant roots were dug out from the pot soils, and then shaken to remove loosely attached soils. Only the soils, closely adhered to the roots, were mechanically brushed and kept as the rhizosphere soils. The fine roots, with diameter smaller than 2 mm, were washed in 1 × TE buffer (added with 0.1% Triton X-100) for 30 s, 75% ethanol for 15 s, 2% bleach for 15 s, and rinsed three times with sterile water (Agler et al., 2016). All samples were put into sterile plastic bags and stored under −80 °C.

DNA extraction and amplicon sequencing

Microbial genomic DNA was extracted from 250 mg soils or 100 mg ground roots using a Powersoil^® Kit (MoBio Laboratories Inc., Carlsbad, CA, USA). Briefly, the sample was lysed via chemical and mechanical homogenization using lysis buffer from the kit and a standard bench-top vortex. The crude lysate was collected and cleaned up to remove PCR inhibitors. Then, the purified lysate was mixed with an equal volume of DNA binding solution, passed and washed twice through a silica spin filter membrane. The silica-bound DNA was eluted using a 10 mM Tris elution buffer and kept for downstream applications. The partial nucleotide sequences of fungal nuclear ribosomal internal transcribed spacer (ITS rDNA) and bacterial 16S rRNA were amplified by polymerase chain reaction (PCR) with fungal primer set ITS5-1737F (5′-GGAAGTAAAAGTCGTAACAAGG-3′)/ITS2-2043R (5′-GCTGCGTTCTTCATCGATGC-3′) and bacteria primer set 515F (5′-GTGCCAGCMGCCGCGGTAA-3′)/806R (5′-GGACTACHVGGGTWTCTAAT-3′), respectively (Caporaso et al., 2011). PCR reaction was carried out using Phusion High-Fidelity PCR Master Mix (New England Biolabs Inc., Ipswich, MA, USA). The PCR cycle was carried out with an initial 94 °C for 5 min, followed by 35 cycles of 94 °C for 45 s, 56 °C for 30 s, 72 °C for 30 s, and a final extension at 72 °C for 10 min. For each sample, the PCR was conducted in triplicate and then pooled to make one PCR mixture. The harvested DNA was verified on 2% agarose electrophoresis gel and purified with Qiagen Gel Extraction Kit (Qiagen, Germany). Sequencing libraries were generated using TruSeq^® DNA PCR-Free Library Preparation Kit (Illumina, San Diego, CA, USA), and assessed on the Qubit@ 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) and Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). The amplicon sequencing was carried out on the Illumina HiSeq2500 platform (Illumina, San Diego, CA, USA) by Novogene Co., LTD (Novogene, Beijing, China).

Bioinformatics analysis

After sequencing, the raw reads were filtered by Trimmomatic version 0.39 (Bolger, Lohse & Usadel, 2014) and cutadapt v1.9.1 (Martin, 2011) to remove low-quality reads (with more than 10% unknown nucleotides, and less than 80% Q-value >20 bases), short sequences (<200 bp), adapters, primers, and poly bases. The left high-quality reads were assembled by FLASH (v1.2.7, http://ccb.jhu.edu/software/FLASH/) with a minimal overlap of 16 bp and mismatch of 5 bp. De-noise and removal of chimeric sequences were processed by dada2 (Callahan et al., 2016) in QIIME2 2020.6 (Bolyen et al., 2019). The amplicon sequence variants (ASVs) were clustered by dada2, and annotated by bayesian classifier using UNITE (https://unite.ut.ee/) and SILVA138 (http://www.arb-silva.de/) as fungal and bacterial reference database, respectively.

Network analysis

The microbial interaction network was constructed using SParse InversE Covariance Estimation for Ecological Association Inference (SPIEC-EASI) R package v1.0.7 (Kurtz et al., 2015). For microbial network inference, ASVs with an average relative abundance lower than 0.02% were filtered out and only top fungal and bacterial ASVs were used. All networks were constructed following the same workflow. The specific construction workflow could be found in an R document (https://www.rdocumentation.org/packages/SpiecEasi/versions/1.0.7), the parameters were set as default. The network layout was shown, and also the network descriptive parameters was calculated, in R package igraph v1.3.2 (https://igraph.org/) based on the Fruchterman-Reingold algorithm.

Statistical analysis

The ASVs table for all samples was generated in QIIME2 2020.6 (Bolyen et al., 2019), and then normalized to the sample with the least sequence number (39,262 for ITS and 37,161 for 16S rRNA sequences; Fig. S1), the alpha diversity indices (Observed species, Shannon, Simpson, ACE, Chao1, Good’s coverage, and Phylogenetic diversity) of each sample and the relative abundances of each taxon were subsequently calculated in the same software. All statistical analyses were conducted in the R software version 4.0.0 (R Core Team, 2022). The estimation of fold change of fungal and bacteria ASVs was performed in DESeq2 (Love, Huber & Anders, 2014). The correlation analysis was carried out in the package Hmisc (https://hbiostat.org/R/Hmisc/). The data presented by heatmap were carried out and plotted in the package pheatmap (Kolde, 2019). The ANOVA followed by Tukey post-hoc pairwise test was used for comparisons among the bulk soils, rhizosphere soils, and roots. The P values were adjusted using false discovery rate (fdr) method, P values smaller than 0.05, except for the correlation analysis, were accepted with significance.

The rhizosphere and root selections on fungal and bacterial communities

A total of 97 samples, including 17 samples of bulk soils, 42 of rhizosphere soils, and 38 of plant roots, were successfully sequenced with both qualified ITS and 16S rRNA sequences (Table S1; Figs. S1A and S1B, respectively). Accordingly, 7,142,267 ITS sequences (mean 73,632 ± standard deviation 12,779 per sample) and 6,937,220 16S rRNA sequences (71,518 ± 11,483) were harvested, respectively. After normalization (at 39,262 for ITS and 37,161 for 16S rRNA sequences) and deletion of singletons, the sequences were clustered into 15,626 (989 ± 448) fungal and 23,423 (999 ± 578) bacteria ASVs, respectively. By compartment, 2,978 (19.1% of all fungal ASVs, Fig. 1A) fungal and 9,762 (41.7% of all bacteria ASVs, Fig. 1B) bacteria ASVs were found in the bulk soils, rhizosphere soils, and roots. In addition, 6,268 (2,978 + 3,290, 40.1%) fungal and 15,942 (9,762 + 6,180, 68.1%) bacteria ASVs were found in bulk soils and rhizosphere soils, only 1, 110 (1,033 + 77, 7.1%) fungal and 1,449 (1,291 + 158, 6.2%) bacteria ASVs found in bulk soils were not found in rhizosphere soils; 4,461 (2,978 + 1,483, 28.5%) fungal and 11,445 (9,762 + 1,683, 48.9%) bacteria ASVs were found in rhizosphere soils and roots, 9,212 (5,922 + 3,290, 59%) fungal and 10,178 (3,998 + 6,180, 43.4%) bacteria ASVs found in rhizosphere soils were not found in roots.

The mean numbers of both the fungal and bacteria ASVs decreased significantly from bulk soils to rhizosphere soils (1,411 ± 271 to 1,220 ± 286, P = 0.03 for fungi; 1,617 ± 323 to 1,277 ± 367, P < 0.01 for bacteria), and further from rhizosphere soils to roots (to 546 ± 234, P < 0.01; to 416 ± 249, P < 0.01). The fungi to bacteria ratios (f:b) of both the observed species and shannon diversity index were not significantly different between bulk soils and rhizosphere soils, but significantly improved from rhizosphere soils to roots (1.0 ± 0.4 to 1.6 ± 0.8, P < 0.01 for f:b of the observed species, Fig. 1C; 0.6 ± 0.1 to 0.8 ± 0.3, P < 0.01 for f:b of the shannon diversity index, Fig. 1D).

The rhizosphere and root selections on fungal and bacterial ASVs

After the dropout of low abundant ASVs (the mean relative abundance <0.01%), 2,249 fungal and 4,523 bacteria ASVs were used for the comparisons between bulk soils and rhizosphere soils, and between rhizosphere soils and roots. Following the criteria (|Log₂FoldChange| ≥ 2 & adjusted P value < 0.01), 56 fungal and two bacteria ASVs were depleted, while 17 fungal and nine bacteria ASVs were enriched in rhizosphere soils in comparison to bulk soils (Figs. 2A and 2B). In addition, 290 fungal and 691 bacteria ASVs were depleted, while 75 fungal and 12 bacteria ASVs were enriched in roots in comparison to rhizosphere soils (Figs. 2C and 2D). Not only the relative abundance of the fungal and bacteria ASVs, the correlation between the fungal and bacteria ASVs was also affected by rhizosphere and root selections. Following the criterion (adjusted P value < 0.01), the number of significant correlations between the top 50 fungal and the top 50 bacteria ASVs increased from 99 in bulk soils to 110 in rhizosphere soils, and further increased to 142 in roots (Fig. 3).

The rhizosphere and root selections on microbial interaction networks

The interaction networks were constructed for the communities of fungi (n = 92, Fig. 4A), bacteria (n = 83, Fig. 4B), and combined fungi and bacteria (n = 175, Fig. 4C) in bulk soils (left), rhizosphere soils (middle), and roots (right), respectively. For all the three communities, the complexity of the networks increased after rhizosphere and root selections (Fig. 4). In specific, the density of the fungal networks increased from 0.017 in bulk soils to 0.028 in rhizosphere soils, and further increased to 0.035 in roots; the density of the bacteria networks increased from 0.021 to 0.023, and further to 0.025; the density of the combined network increased from 0.021 to 0.026, and further to 0.028 (Fig. 5A). Similarly, the number of edges of the fungal networks increased from 72 in bulk soils to 119 in rhizosphere soils, and further increased to 145 in roots (Fig. 5B); the number of edges of the bacteria networks increased from 71 to 79, and further to 84 (Fig. 5C); the number of edges of the combined networks increased from 319 to 396, and further to 430 (Fig. 5D).

As the edge was represented by microbe-microbe interaction, all the edges were divided into either positive (orange) or negative (light blue) microbial interactions. All the microbial networks were dominated by positive interactions (Fig. 4), the mean ratio of positive interactions was 74.2% ± 9%, and the mean ratio of negative interactions was 25.8% ± 9% (Figs. 5B–5D). After the rhizosphere and root selections, both the number of positive (from 61 in the bulk soils to 102 in the rhizosphere soils, and further to 112 in the roots) and negative (from 11 to 17, and further to 33) interactions increased in the fungal networks (Fig. 5B). However, only the numbers of negative interactions increased consistently in the bacteria (from 16 to 26, and further to 36, Fig. 5C) and the combined (from 76 to 99, and further to 140, Fig. 5D) networks. On the contrary, the numbers of positive interaction varied differently in the bacteria (from 55 to 53, and further to 48, Fig. 5C) and the combined (from 243 to 297, and further to 290, Fig. 5D) networks after the rhizosphere and root selections. At last, the variation of the positive and negative interactions (Fig. 5D) was similar to the variation of the intra-kingdom and inter-kingdom interactions (Fig. 5E) in the combined communities. After the rhizosphere and root selections, the number of intra-kingdom interaction increased from 231 to 304, and further to 307; the number of inter-kingdom fungi-bacteria interaction increased from 88 to 92, and further to 123.