Diversity and structure of the root-associated bacterial microbiomes of four mangrove tree species, revealed by high-throughput sequencing

Plant materials and sampling

The sampling site was located at Dongzhaigang National Natural Reserve (110°32′110° 37′E; 19°95′20°1′N), Haikou City, Hainan Province, China (Fig. 1), which harbors mangrove communities with the highest abundance in China. Four different mangrove tree species, A. ilicifolius, B. gymnorrhiza, C. inerme, and L. racemosa were collected on May 3, 2021. For each mangrove species, five biological replicates, each comprising five pooled subsamples were collected, as well as their surrounding soils. Compartments surrounding the plant root system are known as rhizocompartments, which comprise the endosphere (i.e., all inner root tissues), the rhizosphere soil (i.e., the soil immediately surrounding the root) and the non-rhizosphere soil (i.e., the soil not immediately surrounding the root). These distinct compartments collectively influence various physiological and ecological processes that shape the overall functioning of plants. To disclose the bacterial microbiome composition across three rhizocompartments of four mangrove species, the samples we collected from each species were also divided into three compartments: endosphere (R), rhizosphere soil (S), and non-rhizosphere soil (N). They were processed within 12 h after sampling. The soil around the root not adhering directly to the roots was defined as the non-rhizosphere soil fraction. The soil (about one mm thickness) is attached to the root by washing with sterile Phosphate Buffered Saline (PBS) solution was collected as a rhizosphere fraction. To obtain the endosphere biomass, root systems were washed by constant shaking in TE buffer containing 0.1% Triton X-100, followed by washing consecutively in 75% ethanol and 2.5% Sodium hypochlorite to further rinse the root surfaces, and then washed five times (3 min each) with sterile distilled water.

DNA extraction, PCR amplification, and high-throughput sequencing

The root was homogenized in advance by bead beating for 1 min before the DNA isolation (Mini Beadbeater, Jingxin, and Shanghai). Root, rhizosphere, and non-rhizosphere soil (approximately 0.5 g) with five biological replicates used for DNA isolation. The total DNA for each sample was then isolated using the protocol of the MoBio PowerSoil^® DNA isolation kit (Mo Bio Laboratories Inc, USA, catalog No. 12888-100) with slight modifications. Briefly, the kit mainly utilized a bind/wash/elute workflow. The test sample was treated with lysis buffer (Solution C1). Following a brief centrifugation at 10,000×g for 30 s, the resulting supernatant was carefully collected and transferred to a pristine two mL Collection Tube. To effectively eliminate impurities, repeated additions of Solution C2 were employed, with subsequent supernatant collection achieved through centrifugation. To ensure a robust binding of the DNA to the Spin Filter, Solution C4, characterized by high salinity, was introduced and thoroughly mixed with the supernatant. Subsequently, the resulting solution, comprising approximately 675 µL, was loaded onto the Spin Filter. This was accompanied by centrifugation at 10,000 x g for 1 min, with the process repeated three times to enhance the binding efficacy. Following this binding step, the Spin Filter underwent successive washes using Solution C5 (500 µL). Ultimately, the elution of the DNA of interest was achieved through two rounds of centrifugation at 10,000 x g for 30 s each, utilizing solution C6 (50 µL). The quality of DNA was determined by OD(260/280) and OD(260/230) ratios through a NanoDrop spectrophotometer.

The 16S rRNA genes (V3 −V4 region) were amplified using primer pairs (forward primer, 5′-ACTCCTACGGGAGGCAGCA-3′; reverse primer, 5 ′- GGACTACHVGGG TWTCTAAT-3′) as described previously (Wang, Chen & Zhang, 2017). The polymerase chain reaction (PCR) was carried out in a reaction volume of 50 µL, containing 5 µL 10 ×Pfu Buffer, 2 µL 2.5 mM dNTP Mixture, 2 µL each primer set (10 µM), 10 µL 5xGC Enhancer 1 µL Pfu (2.5 U/ µL), 50 ng template DNA. PCR was performed according to a previous study (Zhuang et al., 2020). Before amplicon sequencing, PCR products of all samples were detected by running agarose gel and quantified through a NanoDrop spectrophotometer. The high-quality 16S rRNA amplicons with the same concentrations from each sample were pooled and then sent for sequencing using the Illumina Novaseq 6000 platform by generating PE reads in targeted regions (Biomarker, Beijing, China).

Amplicon-based microbiome data analysis

Valid amplicon sequences were obtained from the verified libraries using the Illumina Novaseq 6000 system. These sequences were then transformed into raw reads through a base-calling algorithm, and assembly was performed using FLASH (Version 1.2.11). To generate high-quality amplicon sequences, a series of tools including Trimmomatic (Version 0.33), Cutadapt (version 1.9.1), Usearch (Version 10.0), and UCHIME (Version 4.2) were employed. During the filtering process with Trimmomatic, sequences were trimmed using a 50-bp moving window approach, while maintaining a quality threshold Q-score of 20. Cutadapt was utilized to identify and remove primer sequences, employing a maximum mismatch ratio of 0.2 and a minimum coverage of 80%. For the assembly of paired-end (PE) reads, Usearch was utilized, requiring a minimum overlap length of 10 bp, a minimum similarity within the overlapping region of 90%, and a maximum accepted mismatch of 5 bp. Subsequently, UCHIME was utilized to eliminate chimeric sequences, using a similar threshold of greater than 80% compared to the query sequence. The resulting high-quality reads, obtained through the aforementioned steps, were used for subsequent analyses. The operational taxonomic units (OTU) as the final valid reads were generated by 97.0% sequence similarity with Usearch software (Version 10.0) and the conservative threshold for OTU filtration is 0.005%. The annotation of feature sequences was performed using a bayesian classifier compared with the SILVA reference sequences (Release 132, http://www.arb-silva.de) (Quast et al., 2013). The sequences matching the reference sequences annotated as “Chloroplast” and “Mitochondria” were removed from the datasets.

Bioinformatics analysis was carried out using the cloud platform BMKCloud (http://www.biocloud.net). QIIME2 (https://qiime2.org/) was used for evaluating the abundance of different species in samples, and the histogram based on frequency distributions at different taxonomic levels was performed using an R package. Subsequently, the α-diversity index of samples was also obtained by the QIIME software, using Student’s t-test to analyze the bacterial diversity and richness of each mangrove species. The ACE and Chao1 indices were used for estimating the bacterial abundance, while the Simpson and Shannon indices were been applied to estimate the bacterial diversity. Principal coordinates analysis (PCoA) was performed to test the distinction in β-diversity through the https://en.wikipedia.org/wiki/UniFrac (Binary Jaccard and Unweighted Unifrac) algorithms to measure the distance (OTU level). Non-metric multidimensional scaling (NMDS) obtained from the Binary-Jaccard distance was applied to investigate the difference between different groups. Biomarkers (significantly different OTUs) between different groups were filtered out by LEfSe (Line Discriminant Analysis (LDA) Effect Size) (Quast et al., 2013). The taxonomic ranks were from phylum to species and the LDA threshold was set to 4.0. A network diagram is used to reveal the correlation. Spearman rank correlation analysis was applied to evaluate the associations and variations in abundant species (—r—>0.9; p < 0.05) and then used for building the network. To predict potential functions based on the 16S rRNA gene amplicon data, the Functional Annotation of Prokaryotic Taxa (FAPROTAX Version 1.2.6) database (http://www.loucalab.com/archive/FAPROTAX) was utilized (Louca, Parfrey & Doebeli, 2016). The FAPROTAX database was used to predict possible functions based on the data of 16S rRNA gene amplicons. FAPROTAX comes with a versatile script (collapse_table.py) for converting taxonomic microbial community profiles (e.g., in the form of an OTU table) into putative functional profiles based on the published and verified culturable bacteria literature. The complete database for FAPROTAX includes over 7600 functional annotations covering over 4600 taxa, and over 80 functions, being freely available at http://www.loucalab.com/archive/FAPROTAX. It maps prokaryotic taxa (e.g., genera or species) to metabolic or other ecologically relevant functions (e.g., nitrification, denitrification or fermentation).

General analysis of the sequencing data

The sequencing of the bacterial 16S rRNA amplicons of four distinct mangrove species, A. ilicifolius, B. gymnorrhiza, C. inerme, and L. racemosa yielded 4,314,420 PE reads in total from 60 samples. After quality filtering and read merging, 3,244,252 high-quality sequences in total were obtained (Table S1). With a threshold of 97% sequence identity and deleting low-abundance OTUs, 2636 OTUs were detected (Table S1). Among them, 2,119 bacterial OTUs were in common among the four distinct mangrove species in all samples (Fig. 2A), and 1,608 bacterial OTUs were shared among the three different rhizocompartments (Fig. 2B). Nine unique OTUs of C. inerme were found and the number was much higher than the other three mangrove species (Fig. 2A), this may be related to the fact that the sampling site of C. inerme was on the shore, and the soil was dry, while the soil of other three mangrove species was muddy. Meanwhile, more unique OTUs were observed in the endosphere than that in the other two compartments (Fig. 2B). This is probably because the endophytic microorganisms in the root are affected not only by the soil environment but also by physiological conditions, growth stages, and metabolites of the plant. The sample-based Shannon curves and OTU species accumulation curves of bacterial OTUs nearly achieved the saturation state (Fig. S1), and Good’s coverage obtained for all the test samples was over 99 percent (Table S2), manifesting that the sequencing depth was enough to cover most bacterial ranks in all the tested samples.

Composition of bacterial communities

We compared the composition analyses of the root-associated compartments in four different mangrove species, based on the 16S rRNA amplicon dataset. The total OTUs in all samples were classified and designated into 34 phyla, 94 classes, 232 orders, 385 families, and 646 genera.

The abundant bacteria at the taxonomic ranks (phylum, class, and genus) among three different compartments of four different mangrove species are shown in Figs. 3A–3C. At the phylum level, the dominant bacteria in all compartments of all samples (endosphere, rhizosphere soil, and non-rhizosphere soil) was Proteobacteria, accounting for 62.60%, 47.99%, and 45.53%, respectively, while the second dominant bacteria in each compartment was different. As the core compartment, the endosphere was distinctive. The second most common phylum was Actinobacteria (20.57%) in the root endosphere compartment, while it was Chloroflexi in the other two compartments (Fig. 3A). At the class level, Gammaproteobacteria was the most abundant which was 58.14% in endosphere, 23.04% in rhizosphere soil, and 19.21% in non-rhizosphere soil, respectively. Actinobacteria (18.64%) was the second predominant class in the endosphere compartment. Alphaproteobacteria was the second dominant class in other compartments except in roots, the relative abundances were 12.85% in the rhizosphere and 16.67% in the non-rhizosphere respectively (Fig. 3B). At the genus classification level, we found significant distinctions across three root-associated rhizocompartments. The root endosphere had a significantly higher proportion of Marinomonas, Vibrio, and Acidothermus than the other two rhizocompartments, whereas uncultured members of Gammaproteobacteria, Gaiellales, and Rhodobacteraceae were almost depleted in the endosphere (Fig. 3C).

The abundant bacteria at the classification levels (phylum, class, and genus) in four different mangrove species are shown in Figs. 3D–3F. The phylum-level similarity of the bacterial community was found in different mangrove species, although the proportion was different. C. inerme had a larger proportion of Actinobacteria, Acidobactera, and Gemmatimonadetes, while C. inerme had a smaller proportion of Proteobacteria. However, in B. gymnorrhiza, the bacteria with the lower abundance belonged to Actinobacteria. At the class level, the bacterial composition of C. inerme was quite different from the other three species, with a smaller proportion of Gammaproteobacteria and Deltaproteobacteria. In the top 15 abundant classes, most of the taxa had a similar ratio. Notable genus-level differences among the four mangrove species were reported in this study. In A. ilicifolius and B. gymnorrhiza, the predominant genus was Vibrio, while the dominant bacteria in L. racemosa was Marinomonas. In C. inerme, an uncultured terrestrial bacterium was dominant, which was consistent with the root’s surroundings. The mangrove tree C. inerme grows in drier, higher land, and is submerged by sea water for less time, so its root-associated microorganisms are characteristic of land plants.

At the phylum level, Actinobacteria was the second common phylum in the root endosphere compartments of four different mangrove species except for B. gymnorrhiza, (Fig. 4A). B. gymnorrhiza trees selected in our study were soaked in seawater and mud all year round, resulting in an oxygen shortage of the roots, thus affecting the growth and reproduction of aerobic actinomycetes. However, at the genus level, the bacterial composition of the same compartment in different mangrove species is quite different from each other (Fig. 4B).

In each compartment of each mangrove species, the microbial composition and structure are similar to that of all samples with some differences. At the phylum level, in all compartments of C. inerme and L. racemosa, Actinobacteria accounts for a greater proportion, whereas it is depleted in A. ilicifolius and B. gymnorrhiza (Figs. 4C–4D). In the four different mangrove plants, the microbial composition in the roots of A. ilicifolius, B. gymnorrhiza, and L. racemosa was significantly different from that in the other two compartments. In C. inerme, the structure of the bacterial community between the root and other compartments are the most similar, which may be because the C. inerme trees we selected grow in a nutrient-poor shore environment, which is less affected by seawater, and thus the connection between the root and the rhizosphere was more closely. At the genus level, the differences between the three compartments were more obvious. Although there are similar groups at the phylum level, the predominant bacteria are quite different at the genus level, which may be relevant to the host and environmental differences between them.

Comparisons of bacterial communities among distinct root-associated compartments

The richness and diversity of microbial communities among different root-associated compartments in four mangrove species were estimated by α-diversity indices, including ACE, Chao1, Shannon, and Simpson (Table S2). In terms of α-diversity, there was no significant difference between rhizosphere and non-rhizosphere microbial diversity whether ACE, Chao1, Shannon, or Simpson indices were selected. However, the diversity between endosphere and the other two compartments were remarkably different (P < 0.05) (Figs. 5A–5D). In general, the diversity in the root endosphere compartment was much lower than that in the other two compartments. This phenomenon was consistent with previous findings, illustrating that the abundance of microbial communities in the rhizosphere was higher compared to the endosphere under a high-nutrient soil environment (Orozco-Mosqueda et al., 2018).

α-diversity between rhizospheric compartments illustrated a decreasing gradient in bacterial diversity from the rhizosphere to the endosphere independent of mangrove species (Figs. 5A–5D and Table S2). Bacterial communities of the endosphere had the lowest α-diversity, while that of the rhizosphere had the highest α-diversity, although α-diversity indices were similar between the two exterior rhizocompartments. As shown in Fig. S2 α-diversity between each mangrove species revealed no obvious difference, although the structure of bacterial communities was different.

The microbial community structure comparisons among different root-associated compartments in four mangrove species were performed by PCoA using the Binary Jaccard and Unweighted Unifrac algorithms (Figs. 6A–6B). In this type of analysis, the distance shows the similarity of bacterial populations. Therefore, samples with high similarity are inclined to cluster together. As shown in Figs. 6A–6B, bacterial communities of the two exterior rhizocompartments were located closer and were more similar to each other than those of the endosphere. Remarkable differences were obtained between the root endosphere and the other two compartments along the PC1 axis (P value <0.05), so the β-diversity of bacteria in the endosphere compartment exhibits a uniqueness from the other two compartments. In addition, the pattern of separation is in accord with a gradient of bacterial populations from the interior of the root into the exterior of the root. The β-diversity of bacteria among four mangrove species exhibits no significant difference (Figs. 6C–6D).

As shown in Fig. 7A, the UPGMA clustering tree and histogram of species composition based on UniFrac distance showed that the root endosphere of four mangrove species was clustered together, while the samples from the other two compartments were clustered together separately, but there were some overlaps between these two compartments samples. These results were consistent with PCoA (Figs. 6A–6B).

Non-metric multidimensional scaling analysis (NMDS) using the Binary-Jaccard algorithm was carried out to compare the bacterial community structure in different compartments or different mangrove species (Figs. 7B–7C). NMDS results in Fig. 7B well revealed the difference in bacterial community structure in root compartments (R² = 0.36, P = 0.001). The bacterial communities in the endosphere compartment differed significantly as compared to the other two compartments, suggesting the selective entry and exit of bacteria communities on the surface of the root, while the differences in bacterial communities among different mangrove species were not significant (R² = 0.16, P = 0.002).

Differential analysis between different compartments and different mangrove species

LEfSe analysis was applied to seek the biomarkers with statistically significant differences between the groups. LEfse analysis sought out significantly different genera abundant in three root-associated compartments or four distinct mangrove species (P < 0.05, LDA score>4.0). The bacteria with high relative abundances among root-associated compartments exhibited remarkable differences. A histogram of the LDA value distribution showed that 36 taxa were enriched including seven genera (Vibrio, Marinomonas, an uncultured bacterium of the family Enterobacteriaceae, Escherichia Shigella, Shewanella, Acidothermus, and an uncultured bacterium of the family Spongiibacteraceae) in the endosphere compartment, 20 taxa were enriched including one genus (Sulfurimonas) in the rhizosphere compartment, while forty-one taxa were enriched including nine genera (an uncultured bacterium of the family Rhodobacteraceae, an uncultured bacterium of S0134_terrestrial_group, an uncultured bacterium of Gitt_GS_136, an uncultured bacterium of the class Gammaproteobacteria, an uncultured bacterium of the order Actinomarinales, an uncultured bacterium of the class Alphaproteobacteria, an uncultured bacterium of the family Xanthobacteraceae, an uncultured bacterium of the order Gaiellales, and Subgroup 10) in the non-rhizosphere compartment (Fig. 8A and Fig. S3). The result indicated these genera enriched in each compartment could be used as biomarkers. In the four mangrove species, two biomarkers were identified in the samples of A. ilicifolius, eighteen biomarkers were found in the samples of B. gymnorrhiza, one biomarker was identified in the samples of C. inerme, while twelve biomarkers were found in the samples of L. racemosa (Fig. S3).

The ternary plot uses an equilateral triangle to describe the relationship among the ratio of attributes of three rhizocompartment. In our study, a ternary plot was used to compare the richness of the top five phyla in three different root-associated compartments (endosphere, rhizosphere, and non-rhizosphere). As shown in Fig. 8B, the top five phyla are widely distributed in all three compartments, but the relative abundance was quite different. The endosphere compartment (R) had the highest abundance of Proteobacteria and Actinobacteria compared to the other two compartments whereas the distribution of the top five phyla between the rhizosphere (S) and non-rhizosphere (N) compartments were more similar. These results were consistent with the composition analysis of microbial communities discussed above.

Correlation network analysis of bacterial communities between mangrove species

The co-occurrence networks of the endosphere (R) and rhizosphere (S) compartments were analyzed by the Pearson algorithm. As shown in Fig. 9A, the endosphere compartment showed a completely different network structure compared to the rhizosphere compartment. Among the top ten most abundant genera, the network in the endosphere compartment consists of 46 edges. An uncultured bacterium of the family Moraxellaceae, an uncultured bacterium of the order Alteromonadales, an uncultured bacterium of the family Enterobacteriaceae, and Acidothermus were found to be the hub genera (≥6 edges per node) in the network in which an uncultured bacterium of the family Enterobacteriaceae and Acidothermus were the biomarker taxa. Vibrio, the most dominant genus in the samples, and also one of the biomarker taxa, had a negative relationship with Acidothermus and Mycobacterium but had a positive relationship with Shewanella. In the rhizosphere compartment, the network structure with more edges was more complicated compared to the endosphere compartment. As the most dominant genus and the only biomarker genus, Sulfurimonas had the highest degree of connections (13 edges). It had a positive relationship with seven genera and a negative relationship with six genera (Fig. 9B). Interactions between microbes have a great impact on the growth and production of metabolites, therefore, further validation is needed in future studies.

16S functional genes prediction of bacterial communities of root-associated compartments

In our study, functional prediction of the top 25 biological functions of the bacterial populations was performed using the FAPROTAX dataset (Figs. 10A–10D). Among different mangrove species, nitrate reduction in the endosphere was much stronger than that in the other two compartments. This result was consistent with the differences in the bacterial community since the genera Vibrio, Marinomonas, and unclassified Enterobacteriaceae were involved in the function of nitrate reduction (Fig. 10A). In the root endosphere of A. ilicifolius and L. racemosa, nitrogen and nitrate respiration had a higher proportion than that of other compartments, which was consistent with the corresponding functions of genera Marinomonas and unclassified Enterobacteriaceae (Fig. 10B). Similarly, in the non-rhizosphere and rhizosphere compartments, the function of photoautotrophy was generally higher than that of the endosphere compartment, while the function of photoautotrophy with a great proportion in all the compartments of C. inerme, which was also related to the microbial composition, as these compartments were rich in a high percentage of related bacteria, i.e., Chloroflexi (Figs. 10B–10C).