The Yellow River Delta is the second largest estuary delta in China, which is one of the areas with the highest biodiversity in the world^[1]. The succession process occurs from ocean to terrestrial land in this region, along with a decrease in the soil salt content and an increase in the organic matter content. At present, few studies of the microbial characteristic associated with different successional process in the wetlands have been reported. This estuarine delta is naturally characterized by extensive coverage of primary salinization, and different types of halophytes grow along the natural saline gradient. Saline-alkali stress is one of the main abiotic factors that greatly restricts the survival, growth, and productivity of plants by affecting their physiological and development processes^[1-2]. However, soil microorganisms can promote the growth and yield of the plants by increasing their tolerance to salinity and nutrient availability^[3-5]. In turn, halophytes could improve the microenvironments of the soil, which promotes the survival of rhizosphere microorganisms^[6].

Soil microorganism play critical roles in maintaining the ecological balance and have been regarded as sensitive indicators of the soil environments^[7-8]. Microorganisms play important roles in biogeochemical cycling in the rhizosphere via nitrogen fixation, phosphorus uptake, and nutrient cycling^[9-10]. The composition and abundance of soil bacteria are often related to soil chemistry and plant species^[11-14]. Many studies have demonstrated that soil pH and salinity are the main factors that affect the composition of microbial communities^[15-18]. In a global survey of bacterial communities from terrestrial environments, salinity has been reported to be the major environmental determinant of microbial community composition^[19-20]. The structure of soil bacterial communities has also been shown to be associated with differences in soil pH^[21-22]. In addition, plant species can influence the soil characteristics and microbial communities through the excretion of root exudates and oxygen into the rhizosphere^[23-25]. However, the combined effects of plant species and soil properties on soil microorganisms are poorly understood.

The Yellow River Delta forms a natural saline-alkali gradient from the coastal to terrestrial land. Suaeda glauca (SG), Glycine soja and Phragmites australis (GP) are three typical plant communities along the succession region in the Yellow River Delta since 1990s. In this study, we compared the microbial community compositions and functional gene contents of rhizosphere microbial communities of Suaeda glauca (SG) and Glycine soja-Phragmites australis (GP). We compared halophyte rhizosphere microbial communities with the bulk soils of barren wetland in the coastal area of the Yellow River Delta. Therefore, the aim of this study was to investigate the combined effects of different vegetation covers and soil factors on both microbial communities and functional genes in Yellow River Delta rhizosphere soils.

1 Materials and Methods 1.1 Site description and sampling

Three sampling sites were located in the Yellow River Delta (119.1°E, 37.7°N) in Dongying City, Shandong Province, China. The rhizosphere soils were collected as described previously with a few modifications^[26]. A nested five-point sampling method was used. At each sampling site, surface sediments (0-5 cm) were sampled at five points, and the distance between each point was 5 m. Three subsamples were collected in each sampling site. After that, all the subsamples of each site were manually mixed to homogeneity and immediately transported to the laboratory. The soil samples were divided into two parts. One was stored at 4 ℃ to analyze its soil properties, while the other was stored at -80 ℃ until DNA extraction.

1.2 Soil physiochemical properties

The soil pH was measured as previously described using a pH meter (Mettler Toledo Instruments, Shanghai, China)^[27]. The electrical conductivity (EC) was measured using a conducting meter (DDS-11A, Shanghai, China). The total organic carbon (TC), total nitrogen (TN), and total phosphorus (TP) were measured as previously described^[28].

1.3 Metagenome sequencing and assembly

Illumina HiSeq-2000 sequencing of metagenomic DNA from the three samples generated 239489174 raw reads with a total of 24.19 billion bps. Sequences were cleaned and assembled using SeqPrep (https://github.com/jstjohn/SeqPrep), Sickle (https://github.com/najoshi/sickle) and SOAP^[29]. Open reading frames (ORFs) were predicted using MetaGene software^[30]. The ORFs with a length of ≥100 bp were selected and translated into amino acid sequences using the NCBI database. The predicted genes (with 95% sequence identity and 90% coverage) were clustered by CD-HIT software^[31], and the longest sequences from each cluster were selected as representative sequences to build a non-redundant gene catalog. The reads were then mapped to the representative sequences with 95% identity using SOAP aligner^[32], and the gene abundance in each sample was evaluated.

Representative sequences of the non-redundant gene catalog were aligned to the NR database using BLASTp for taxonomic annotations. The reads were also used for the taxonomic assessment by screening for SSU rRNAs using Metaxa2^[33]. The non-redundant gene catalog was aligned against the Non-supervised Orthologous Groups (eggNOG) database using BLASTp for Clusters of Orthologous Groups (COG) annotation with an E-value cut-off of 1e^-5[34]. The non-redundant gene catalog was then aligned against the Kyoto Encyclopedia of Genes and Genomes (KEGG) using BLASTp against the KEGG database with an E-value cutoff of 1e^-5[35].

1.4 Accession number

The raw sequencing data have been deposited into the NCBI Sequence Read Archive (SRA) database. The SRP accession number is SRP253127.

2 Results and Discussion 2.1 Soil properties

Soil samples were taken from three sample sites, where the vegetation was Suaeda glauca (SG), Glycine soja-Phragmites australis (GP) and barren wetland (BW). The chemical composition of the soils is summarized in Table 1. All the soil samples were highly saline and had a high pH. The pH values are alkaline and decline along from the land to coast. The electrical conductivity (EC) values varied from 7.42 to 8.98. A previous study showed that the main component of soil salinity is sodium chloride, which exhibits a significant linear correlation with the EC in Yellow River Delta^[36]. The total nitrogen and total carbon were highest in GP, which could have resulted from the improvement of soil nutrients by the legume Glycine soja. It has been estimated that soil microorganisms require a C:N ratio of approximately 25:1 to reach 40% of their growth efficiency^[37]. The C:N ratio varied from 14.98 to 20.43 in the three samples, indicating that the microbial communities facilitated a high degree of organic nitrogen mineralization of the soil organic matter.

2.2 Microbial taxonomic analysis at the phylum level

We first explored the microbial community composition by analyzing the reads using the ribosomal small subunit (SSU)^[38]. The results demonstrated that the niche contains very high bacterial abundance (90%-93%), with a significantly lower proportion of Eukaryota (1.90%-2.88%) and Archaea (0.43%-3.56%) (Figure 1-A). The populations of Eukaryota and Archaea are found in low amounts in the three soil samples, which is consistent with the results of other environments^[39-40]. The taxonomic analysis of SSUs of the bacterial fraction revealed that Proteobacteria was the most prevalent phylum represented in the metagenome, while other phyla, such as Acidobacteria, Chloroflexi, Actinobacteria, Bacteroidetes, Planctomycetes, and Gemmatimonadetes, were present as relatively minor contributors to the total bacterial abundance (Figure 1-A).

We further analyzed the microbial community structure of the three samples by comparing the relative reads number. The results showed that there were 10 phyla in common that comprised a total of 90.32%, 85.05% and 85.95% in SG, GP, and BW, respectively. As shown in Figure 1-B, the dominant phyla in the three sample soils were Proteobacteria, Firmicutes, Actinobacteria and Chloroflexi. Proteobacteria was the most abundant phylum in all samples, accounting for a total of 71.5%, 61.4% and 55.5% of SG, GP and BW, respectively. Many studies have shown that Proteobacteria predominate in different soil types, which could indicate an involvement in the biogeochemical process in ecosystems^[41-44]. The differences in distribution are probably related to the organic matter content in the soil, since the relative abundance of Proteobacteria can be influenced by nutrients availability^[45].

The next most abundant phylum varied among the three samples, Bacteriodetes (5.87%) in SG, Actinobacteria (5.88%) in GP, and Gemmatimonadetes (7.89%) in BW, respectively. These results are consistent with those of previous studies, in which Bacteroidetes, Actinobacteria and Gemmatimonadetes were reported to be prevalent in the saline-alkaline sediment environment^[46-48]. Bacteriodetes and Actinobacteria have been demonstrated to thrive in plant soils, while Gemmatimonadetes were enriched in barren soils in the Yellow River Delta^[49]. The relative abundance of Bacteriodetes was decreased from 5.87% (SG) to 0.6% (GP), whereas that of Acidobacteria was increased from 0.50% (SG) to 5.88% (GP) and that of Nitrospirae was increased from 0.32% (SG) to 2.49% (GP). Since Acidobacteria serve as decomposers, they participate in the cycling of organic matter derived from plants, fungi and insects^[50]. In this study, the relative abundance of Acidobacteria was higher in the rhizospheres of GP, which could be related to the higher contents of root exudates in that region.

2.3 Microbial taxonomic analysis at the genus level

The richness of the bacterial community at the general level is illustrated in Figure 2. Twenty-four main genera (each with a relative abundance of > 1%) were detected in the three soil samples, and the main genera represented 15.83%-22.88%. The composition of bacterial communities differed substantially across different samples. Marinobacter (3.94%), Streptococcus (3.92%), Alcanivorax (2.83%), Gemmatimonas (1.73%), Pseudomonas (1.25%), Steroidobacter (1.09%), Thioalkalivibrio (1.07%), and Fluviicola (1.03%) were the dominant genera in SG. Sinorhizobium (10.82%), Nitrospira (2.32%), Pseudolabrys (2.25%), Pyrinomonas (1.87%), Bradyrhizobium (1.68%), Streptococcus (1.49%), Lysobacter (1.40%), and Rhodoplanes (1.06%) were the dominant genera in GP. Gemmatimonas (6.81%), Kiloniella (1.51%), Streptococcus (1.38%), Alcanivorax (1.38%), Fodinicurvata (1.37%), Rhodovibrio (1.25%), Azospirillum (1.14%), and Thioalkalivibrio (1.00%) were the dominant genera in BW. Among the top 24 genera, the dominant genera were composed of Sinorhizobium, Gemmatimonas, Marinobacter and Streptococcus, many species of which have been reported to have beneficial effects for C or N cycling in the soil^[48]. In contrast with those in the BW and SG, the main microbes in the communities of GP at later succession stages included the genera of Sinorhizobium (10.82%), Nitrospira (2.32%), and Bradyrhizobium (1.68%). These species have been demonstrated to play important roles in nitrogen cycle^[51]. The legume Glycine soja appears to facilitate the expansion of the Sinorhizobium species in GP, since this microorganism is a nitrogen-fixing mutualistic symbiont of the legume. The relative lower population of Bradyrhizobium could be attributed to the influence of saline soil. Previous studies have identified that Sinorhizobium as the dominant species in alkaline-saline soil, whereas Bradyrhizobium was found to be dominant in neutral to acidic soil^[52-54]. Nitrospira plays a role in nitrification, which is also exist in YRD agriculture soils in previous studies^[55-56]. In concert, these results suggest that different vegetation covers and saline soils have an important impact on the bacterial communities.

2.4 Distribution of putative metabolic genes of microbial communities

To better understand the functional potential of the microorganisms represented in the three sampling sites, the metabolic functions of the microbial communities in SG, GP and BW were investigated by aligning to the eggNOG databases. The results demonstrated that the number of total predicted functional genes of SG and GP was 25.2% and 71.2% lower than those of the BW, respectively (Figure 3-A). In addition to the function of unknown genes, the main functional gene classes in the three samples were amino acid metabolism, energy production and conversion, and inorganic ion transport and metabolism genes (Figure 3-A).

Furthermore, the predicted genes were aligned to the KO (Kyoto Encyclopedia of Genes and Genomes Orthology) database. The results demonstrated that the total number of predicted functional genes of SG and GP was 19.8% and 68.0% lower than that of the BW, respectively (Figure 3-B). Many genes were attributed to three categories, namely amino acid transport and metabolism, carbohydrate metabolism, and energy metabolism (Figure 3-B). Moreover, we compared the functional genes at KEGG level 3 using a heatmap analysis. As shown in Figure 3-C, a higher number of sequences in quorum sensing and ABC transporters was found in GP compared with those of SG and BW. The quorum sensing genes were potentially implicated in regulating the legume-rhizobia nodulation symbiosis process^[57]. The ABC transporters genes could be related to the high EC values of the soil sample in GP, since the ABC transporter proteins have been demonstrated to play a role in osmotolerance^{[46, 58]}. In summary, our data show that the vegetation types and soil characteristics play important roles in the structure of bacterial communities.

We further compared functional genes that participate in the nitrogen cycling in the three samples. As shown in Figure 4, the ATP-dependent glutamine synthetase-encoding gene glnA had the highest abundance in each sample. The next abundant genes were gltB and gltD, which encode glutamate-oxoglutarate amidotransferase and can cooperate with glnA in indirect ammonia fixation. The abundance of N fixation gene (nifH) in SG and GP was higher than that of BW. And more dissimilatory N reduction gene (nasA, nrfA, and nirADK) was detected in SG than that of GP and BW. For the ammonification process, the abundance of gdh of SG was lower than those of GP and BW. This observation could be related to C:N ratio in the soil, as the ammonification process can be influenced by nitrogen availability^[59]. The mere presence of gene does not necessarily mean the fuction, and metatransciptomic analysis and cultivation of targeted functional guilds are warranted for further study.