Spatial and Species Variations of Bacterial Community Structure and Putative Function in Seagrass Rhizosphere Sediment

Seagrasses are an important part of the coral reef ecosystem, and their rhizosphere microbes are of great ecological importance. However, variations in diversity, composition, and potential functions of bacterial communities in the seagrass rhizosphere of coral reef ecosystems remain unclear. This study employed the high-throughput sequencing based on 16S rDNA gene sequences and functional annotation of prokaryotic taxa (FAPROTAX) analysis to investigate these variations based on seagrass species and sampling locations, respectively. Results demonstrated that the seagrass rhizosphere microbial community was mainly dominated by phylum Proteobacteria (33.47%), Bacteroidetes (23.33%), and Planctomycetes (12.47%), while functional groups were mainly composed of sulfate respiration (14.09%), respiration of sulfur compounds (14.24%), aerobic chemoheterotrophy (20.87%), and chemoheterotrophy (26.85%). Significant differences were evident in alpha diversity, taxonomical composition and putative functional groups based on seagrass species and sampling locations. Moreover, the core microbial community of all investigated samples was identified, accounting for 63.22% of all obtained sequences. Network analysis indicated that most microbes had a positive correlation (82.41%), and two module hubs (phylum Proteobacteria) were investigated. Furthermore, a significant positive correlation was found between the OTUs numbers obtained and the functional groups assigned for seagrass rhizosphere microbial communities (p < 0.01). Our result would facilitate future investigation of the function of seagrass rhizosphere microbes.


Introduction
Coral reef plays important roles in the marine ecosystem for their high biodiversity, productivity, and economic and ecological services to millions of people [1]. Their biogenic carbonate structures are accumulated over many years by calcifying organisms, including corals and algae [2]. However, climate change (e.g., ocean acidification and global warming) and human activities both exert direct and indirect effects on the coral reef locally and globally, leading to a wide decline in the coral reef ecosystem [3,4]. Seagrass meadows spatial scales, were selected for seagrass rhizosphere microbe investigation ( Figure 1). Seagrass forms dense meadows in all investigated coral reef ecosystems. Daya Bay, a shallow semi-closed bay (22.54 • E, 114.45 • N) with an area of~600 km 2 , is under the typical subtropical marine climate. Numerous coral colonies distribute sporadically within the bay around the inner islands, such as Dalajia and Xiaolajia Islands. Seagrass Halophila ovalis (H) is the only species forming seagrass meadow in DY Bay. For the Luhuitou fringing reef of Sanya Bay (18.2 • E, 109.47 • N), seagrass Thalassia hemprichii (T) appeared following coral reef degradation, and it is the only seagrass species detected in this fringing reef. There are four seagrass species (Cymodocea nodos (C), T. hemprichii, H. ovalis (H), and Syringodium isoetifolium (S)) growing in one seagrass meadow in the Xisha Islands (XS) (16.84 • E, 112.34 • N). Besides, the Nansha Islands (NS) are in the southern South China Sea and under the effect of a typical tropical climate. Two seagrass species, H. ovalis and T. hemprichii, were collected from Nanxun Reef and Chigua Reef, respectively. The two species are the dominant seagrass species of where they are growing. The seagrass rhizosphere samples were named by combing the sampling locations and the seagrass species, e.g., DYH indicates the sample collected from Daya Bay and seagrass species is H. ovalis. Four investigated coral reef ecosystems are all under the effect of the East Asian monsoon (Northeast in the winter, and Southwest in the summer). The average atmospheric temperature in Sanya Bay was 30.74 • C, with warm summers (34.75 • C) and cold winters (27.20 • C) [24]. The temperature range and mean value for Daya Bay were 14.4-32.4 • C and the annual mean temperature was 22.4 • C, respectively [25]. The temperature range for Yongxing Island was from 15.5 to 37 • C, and the annual average temperature was 26.5 • C [26]. For the Nansha Islands, the temperature only slightly varies throughout the year, with an annual average of 28.1 • C [27].

Sample Collection and Physic-Chemical Property Measurement
Samples of seagrass rhizosphere sediment were collected in triplicate for each sample, and 24 samples in total were obtained. Sediment collection and physic-chemical property measurements were performed, according to Ling et al. [23]. Environmental parameters, e.g., pH, salinity, and DO, were simultaneously measured in situ. The overlying seawater column's temperature and salinity above the seagrass meadows were measured by the YSI 6600V2 water quality sonde TM (YSI, Yellow Springs, OH, USA). A portable pH/DO meter (Thermo Fisher Scientific, Inc., Beverly, MA, USA) measured the dissolved oxygen (DO) concentrations and pH values. The sampling collection was performed at the low tide of about 50~100 cm. Ten (for T. hemprichii, Cymodocea nodos and Syringodium isoetifolium) or thirty (for H. ovalis, small plant size) seagrass shoots containing leaves, rhizomes, and roots were collected manually with great caution to avoid the damages, and distance between them was about one meter. After the seagrass plants were separated from the substrate, the loosely attached sediment was gently shaken from the roots manually with approximately 1 mm of sediment still attached to the roots. Samples for microbial and physicochemical analysis were immediately put in an icebox (4 • C) in the dark and then transported to the laboratory or the cruise ship within three hours. Samples collected in Daya Bay were preprocessed in Daya Bay Marine Biology Research Station, Sanya Bay in Tropical Marine Biological Research Station in Hainan, the Xisha Islands in Xisha Deep Sea Marine Environment Observation and Research Station, and Nansha island on the survey vessel. The seagrass roots were placed in a 50 mL sterile conical tube containing 25 mL phosphate-buffered saline (PBS) solution and vortexed for 30 s. The obtained turbid solution was filtered through a 100-µm nylon mesh, and then the filtrate was centrifuged to obtain the seagrass rhizosphere sediment [28,29], then the samples were stored at −80 • C in the lab or lipid nitrogen until DNA extraction. Inorganic nutrients in seawater, including ammonium, nitrate, nitrite, and phosphate, and total phosphorus (TP) in sediment were measured in triplicate according to the standard oceanography methods (General Administration of Quality Supervision, Inspection, and Quarantine of the People's Republic of China, 2002). Total carbon (TC), total organic carbon (TOC), and total nitrogen (TN) were measured by CHNS Vario E1 III elemental analyzer (Hanau, Germany), according to Sun et al. (2020) [30].

DNA Extraction, PCR, Library Preparation, and Illumina MiSeq Sequencing
Microbial DNA from approximately 1 g of rhizosphere sample was extracted with EZNA ® Soil DNA kit (Omega Bio-Tek, Norcross, GA, USA), quantified by Qubit 2.0 fluorometer (Invitrogen, Carlsbad, CA, USA), and checked on 1% agarose gel. The 16S rDNA sequencing library preparations and sequencing were conducted at GENEWIZ, Inc. (Suzhou, China). Thirty nanograms of total DNA was used to amplify the target amplicons with a MetaVx TM library preparation kit (Genewiz, Inc., South Plainfield, NJ, USA). The forward primers sequence "CCTACGGRRBGCASCAGKVRVGAAT" and reverse primers "GGACTACNVGGGTWTCTAATCC" were used to amplify the V3-V4 region of 16S rDNA [31]. Two primer pairs were used to amplify more microbes, and then more variations of the microbes were identified in the samples. The PCR mixture contained 20 ng DNA template, 2.5 µL 10× TransStart Buffer, 0.5 µL 2.5 U/µL TransStart Taq, 2.0 µL dNTPs (2.5 mM each), an d2.5 µL 1× primers mix. Distilled H 2 O was used as the template for negative controls. The PCR protocol consisted of an initial 3-min denaturation at 94 • C, followed by 14 cycles of denaturing at 94 • C for 50 s, annealing at 57 • C for 90 s, extension at 72 • C for 10 s, and completed with a final extension at 72 • C for 5 min. The libraries for downstream NGS sequencing were constructed according to Li et al. (2017) [31], and then were loaded and sequenced on an Illumina MiSeq instrument according to the manufacturer's instructions (Illumina, San Diego, CA, USA). Sequencing was performed using a 2 × 300/250 paired-end (PE) configuration. All the raw sequences obtained from this study were deposited in the NCBI sequence read archive (SRA) under accession number PRJNA497291.

Statistical Analysis
The resultant amplicon sequencing, bioinformatics, was mainly performed in R package EasyAmplicon v1.0 according to Liu et al. [32]. The raw amplicon paired-end reads were grouped based on their barcode sequences and then were merged to obtain amplicon sequences. The joined sequences were then quality-filtered by "usearch10 -fastq_filter -fastq_maxee_rate 0.01", and the barcode and primers sequences were removed during this process. Criteria for the quality filtering were as following: sequence length >200 bp, ambiguous bases, mean quality score >20. All representative sequences were grouped into operational taxonomic units (OTUs) using the clustering program USE-ARCH (https://www.drive5.com/usearch/ (accessed on 17 August 2020)) against the Silva 138 Reference database [33]. The singletons and chimeric sequences were also removed. The obtained representative sequences of each OTU were used to construct the phylogenetic tree by FastTree after aligned by sequence alignment program Multiple Alignment using Fast Fourier Transform (MAFFT) 7.0 [34][35][36].
All OTUs identified as belonging to the chloroplast, mitochondria, Archaea, and Eukaryote were removed from the data set before downstream analyses. The remaining sequences were rarefied before calculating alpha diversity statistics analysis, and all statistical analyses were conducted in the open source statistical language R (R Core Team, 2018, version 3.4.4). The alpha diversity (observed OTU richness, the Shannon index, Simpson index, and phylogenetic diversity (PD) were calculated in R ("vegan"). Venn diagram plot shared OTUs, and Unique OTUs were generated by the R package "VennDiagram" [37].
Non-parametric tests (999 permutations) of the multiple-response permutation procedure (MRPP) were performed based on weighted UniFrac distance, Bray-Curtis distance, Euclidean distance, and Sorensen distance to measure the dissimilarity of microbial communities, respectively [38,39]. Gene functional annotation based on all obtained sequencing was performed against the FAPROTAX dataset; the dataset and associated software are available at http://www.loucalab.com/archive/FAPROTAX/ (accessed on 17 August 2020) [40]. FAPROTAX dataset is a manually constructed database that maps prokaryotic taxa (e.g., genera or species) to putative functions based on the literature on cultured representatives [40]. The core microbial community in this study was defined as those present in all samples, including all three replicates [11], and their network analysis was further analyzed based on a random matrix theory (RMT)-based method according to Deng [41] and Luo [42]. Network analysis of the core microbial community was processed using by molecular Ecological Network Analysis Pipeline (http://ieg4.rccc.ou.edu/MENA/ (accessed on 17 August 2020)) and visualized in an open-source software platform Cytoscape (v.3.6.0).

Sampling Locations and Physicochemical Properties of the Seagrass Rhizosphere
The sampling location and all the measured physicochemical parameters of investigated samples are shown in Table 1. All four seagrass species across four coral reef ecosystems were studied. There were variations in the physicochemical parameters for different samples. For instance, pH values varied between 8.14 and 8.24 ± 0.03. The salinity range was from 24.70 ± 0.3 to 34.57 ± 0.08 practical salinity unit (PSU), with the lowest value detected at SYB. The highest concentration of nitrate was recorded in the region DYB as 0.153 ± 0.004 mg/L, while SYB had the highest concentration of ammonium (0.121 ± 0.002 mg/L). The range of the DO was from 5.82 ± 0.05 to 9.64 ± 0.59 mg/L. The highest TN, TC, and TOC concentration existed in the sediment of XSC.

Taxonomy, Phylogenetic Diversity, and Composition of Rhizosphere Bacterial Communities of the Seagrass
Total sequences assigned to all samples were 1,601,000. After removing chimeric and singleton, the remaining sequences were 1,332,470 belonging to 2423 OTUs, and the minimum read of all samples was 33,987. Then, sixteen OTUs related to chloroplast/mitochondria were removed. The resampling depth was 33,980 merged sequences, with 2405 OTUs left after resampling. All OTUs were grouped into 33 phyla and 290 genera. The rarefaction curves showed that the sequencing depth was relatively enough to cover bacterial diversity as all the rarefaction curves almost reached saturation plateaus ( Figure S1), and the minimum coverage of the samples was 98.96% (Table S1). The phylum Proteobacteria accounted for 33.47% of bacterial communities, followed by phylum Bacteroidetes (23.33%) and Planctomycetes (12.47%). At the genus level, the dominant genera were Bacillus (phylum Firmicutes), and the range of its relative abundance for each sample was from 4.61% to 24.04%. The highest abundance was detected in sample XSH (the Xisha Islands), while the lowest sample was DYH (Daya Bay). The other dominant genera were unclassified Desulfobulbaceae (class Deltaproteobacteria), unclassified Rhodobacteraceae (class Alphaproteobacteria), unclassified Bacillaceae 2 (phylum Firmicutes), Lactococcus (phylum Firmicutes), and unclassified Desulfobacteraceae (class Deltaproteobacteria and Desulfopila (class Deltaproteobacteria) ( Figure 2).
The alpha diversity, including PD, Chao1, Richness (Observed OTUs), Shannon, and Simpson index of all samples, was calculated, and is shown in Figure S2. The highest Faith's PD value was detected in the samples of NST with a value of 80.92 ± 7.95, while the lowest Faith's PD was 59.3 ± 5.13 for the samples of DYH. Meanwhile, the highest value of Chao1, Richness, Shannon, and Simpson index also existed in the samples of NST. Based on seagrass species (C. nodosa, T. hemprichii, H. ovalis and S. isoetifolium) collected from XS, the alpha index of PD, Richness, Shannon, and Simpson demonstrated significant differences with the p-value below 0.05 (Table 2). For different sampling locations, seagrass H. ovalis from Daya Bay, the Xisha Islands, and the Nansha Islands showed significant variations in Richness, Shannon, and Simpson (p < 0.05), while for seagrass T. hemprichii from Sanya Bay, the Xisha Islands and the Nansha Islands exhibited significant in Richness and Shannon (p < 0.05). The beta diversity analyzed by MRPP based on bray-cutis dissimilarity, Euclidean distance, and Sorensen distance demonstrated that there were significant variations (p < 0.05) for four seagrass species from the same location, while no significant variations were detected for seagrass T. hemprichii and H. ovalis from different sampling locations (p > 0.05) ( Table 3).

Potential Functional Roles of Microbial Played in Seagrass Rhizosphere
The FAPROTAX database comprises 7820 members (4724 unique members) belonging to 90 groups. The results of this study showed that there were 920 assignment records to groups, and 427 out of 2405 records (17.75%) were assigned to at least one group. In total, 36 functional groups were represented (i.e., associated with at least one record). As illustrated in Figure S3, the top four dominant functional groups were sulfate respiration (129 records), respiration of sulfur compounds (131 records), aerobic chemoheterotrophy (192 records), and chemoheterotrophy (247 records), respectively. The ranges of relative abundance of these four groups in all the samples were 11.75-21.44%, 11.75-21.48%, 7.72-11.34%, and 3.29-8.20%. Further analysis revealed that the sulfate respiration group mainly consists of nine genera, including Desulfatiglans, Desulfobulbus, Desulfocarbo, Desulfomonile, Desulfopila, Desulfosarcina, Desulfovibrio, unclassified Desulfobacteraceae, and unclassified Desulfobulbacea, while for the respiration of sulfur compounds, 11 genera had participated in this process. The diversity of microbes that participated in the procedure of chemoheterotrophy and aerobic chemoheterotrophy is very high. Ninety and sixty-eight genera of microbes were identified, respectively. Spearman's correlation was employed to test the relationship between alpha diversity (richness) and functional groups obtained, and results showed that there was a positive correlation (r = 0.74, p < 0.01) (Figure 3).

Venn Diagram Analysis of the Variations in Taxonomy Species and Functional Groups
Based on the seagrass species, the OTUs shared by four seagrass species collected in XS were 1451, and each species had its unique OTUs (Figure 4). Among them, sample XSH had the highest unique OTUs with 47, followed by XSC. Meanwhile, seagrass H. ovalis shared 1362 OTUs with the same species from three sampling locations based on the sampling location. Moreover, for seagrass T. hemprichii, the shared OTUs were 1434, with sample NST possessing 192 unique OTUs. As for the functional structure, from the seagrass species perspective, the functional groups they shared were 31, and no unique functional groups were detected. XST had 33 functional groups and followed by NST having 31 functional groups. Moreover, for seagrass H. ovalis, samples from three sampling locations shared all their functional groups (31), while for seagrass T. hemprichii, they shared 28 functional groups. Furthermore, for seagrass T. hemprichii, the lowest number of functional groups is 28 detected in seagrass SYT, and the highest value is 33 from samples XST (Figure 4). Variations in seagrass rhizosphere microbial communities at taxonomical and functional levels were analyzed based on species and locations, respectively. The top abundant 50 genera were included for further taxonomical structure analysis (Tables S2-S4), and all detected functional groups (36) were included for analysis (Tables S5-S7). Most of the investigated genera demonstrated significant differences between different species based on seagrass species (Table S2). Several genera for species' comparison between two seagrass species from the Xisha Islands, such as Desulfopila, unclassified Bacteroidales, and Eudoraea, which showed no significant differences (p > 0.05). For site-based analysis, many of those genera exhibited substantial variations among the sampling sites, such as genus unclassified Desulfobulbaceae and unclassified Chloroflexi. In contrast, Desulfopila, Oceanobacillus, unclassified Syntrophobacterales, Desulfobulbus of seagrass T. hemprichii (p < 0.05), and Desulfosarcina of seagrass H. ovalis showed no significant differences (p > 0.05).
All investigated samples shared many of the functional groups, but significant differences were also detected in the samples at both species-based and location-based levels. For instance, methanogenesis's functional groups, by reducing methyl compounds with H 2 , Hydrogenotrophic methanogenesis, and methanogenesis, could be found in all the seagrass species samples collected from the Xisha Islands. However, none of the analyzed genera and detected functional groups showed both species differences and location differences.

Core Microbial Community in Seagrass Microbial Rhizosphere
The co-occurrence network method was used to explore the interaction between the rhizosphere microbes and to identify the keystone species. In all, 308 of 2405 OTUs were identified as core OTUs shared by all samples. They accounted for 12.81% of all obtained OTUs, and their relative abundance was 61.83% ( Figure 5A). The core OTUs belonged to 14 phyla and 89 genera, and the predominant phyla were Proteobacteria (24.37%), Firmicutes (21.03%), and Bacteroidetes (3.37%). Afterward, 197 OTUs were selected for network analysis, and the correlation network was generated with a coefficient cutoff of 0.760, as determined by the RMT-based algorithm. There was a total of 773 edges (136 negative correlations and 637 positive correlations), and most of the correlation was positive (82.41%) ( Figure 5B). In all, 19 modules were constructed, and the biggest module was Module 1, consisting of 58 OTUs, followed by Module 3 with 44 OTUs. The OTU 114 and OTU 1807 (phylum Proteobacteria) were identified as module hubs (OTU highly connected in the own module) in Figure S4. Modules with more than five OTUs were included for correlation analysis of module eigengenes and environmental factors. The modules' responses to the environment were different, and Figure 5C showed that the measured physiochemical factors were significantly correlated with module eigengenes of Module 2, 3, and 5. TOC, ammonium, and nitrate were negatively correlated with Module 2 while positively correlated with Module 3. However, Module 1 and 4 were not significantly correlated with the environmental parameters.

Variations in the Taxonomical, Phylogenetical Diversity and Composition of Bacterial Communities
Significant variations in PD and taxonomical diversity of bacterial communities could be detected based on seagrass species within a coral reef ecosystem, but only significant variations in taxonomical diversity for the same seagrass species from different sampling locations (Table 2). Moreover, significant taxonomical and phylogenetic variations only existed among different seagrass species collected from the Xisha Islands. Therefore, the coral reef ecosystem's seagrass species may be one important factor in shaping the rhizosphere bacterial communities.
We also found significant differences in the taxonomy composition of rhizosphere bacterial communities at the genus level based on different seagrass sampling locations. This was partly consistent with the investigation result of Cúcio et al. (2016) [16], the result of which demonstrated that significant differences were detected for the same seagrass species from different sampling locations, but no significant differences existed between the rhizobiomes of different seagrass species from the same sampling location. The reason for this phenomenon may be that different seagrass species were included in each study. Three different seagrass species, namely Z. marina, Z. noltii, and Cymodocea nodos, were studied for Cúcio et al. (2016) [16], while four seagrass species (C. nodos, T. hemprichii, H. ovalis, and S. isoetifolium) were examined in our investigation. Another reason for this discrepancy may be the different growth habits. The seagrass habitats for their study was in the intertidal regions, while all the seagrasses in this study were collected from the coral reef ecosystem [16]. Moreover, previous studies have highlighted the importance of temperature in constructing the rhizosphere bacterial community anwhich exhibited seasonal variations [43,44]. Therefore, there may also be seasonal variations in the seagrass rhizosphere bacterial community. More investigation on the temporal scale in the future needs to be performed.
Proteobacteria (class alpha-, beta-, delta-, gamma-, and epsilon-proteobacteria) and the Firmicutes were the two most predominant phyla across the four coral reef ecosystems. Besides, class Deltaproteobacteria accounted for over 20% of all investigated bacterial communities. Cúcio et al. (2016) also reported that the phylum Proteobacteria was the most dominant in the rhizomes of seagrass Z. marina, Z. noltii, and Cymodocea nodosa, with the proportion ranging from 65% to 68%. The existence of plants played a crucial role in shaping the microbial community in the rhizosphere of seagrasses as the seagrass rhizosphere bacterial community composition was quite different from that of the surrounding water and bulk sediment [16]. Besides, seagrass (Z. marina) colonization increased the abundance of the nitrogen fixation bacteria and other bacteria involved in benthic carbon and sulfur cycling [45].
Moreover, some OTUs were peculiar to one coral reef ecosystem, and each coral reef had its own individual OTUs in our study. For instance, OTU1109 was affiliated to class Phycisphaerae SHA-43 belonging to the phylum Planctomycetes and could only be discovered at XS. It may play an important role in the nitrogen cycle by participating in the anammox process, which was assumed as a predominant source of N 2 production in anoxic marine environments [46][47][48]. Moreover, bacteria from the family Rhodothermaceae (phylum Bacteroidetes) were specially retrieved from Sanya Bay, and microorganisms from this family were usually isolated from the extreme environments and exhibited extreme thermophilic or halophilic characteristics [49,50].

The Functional Structure of Microbial Communities in Seagrass Rhizosphere
Seagrass holobionts have been reported to play essential roles in the cycle of sulfur, nitrogen, and carbon, at both microbial structural and functional levels [11,13], and Ugarelli et al. (2019) [51] reported that the seagrass plant and its microbiome were highly interlinked in the cycle of sulfur, nitrogen, and carbon. Likewise, FAPROTAX analysis revealed that many microbes in the seagrass rhizosphere of coral reef ecosystems participated in these processes ( Figure S3).
Previous studies showed that increased sulfide concentration in the sediment caused by the activity of sulfate-reducing prokaryotes was one of the main reasons for seagrass death all over the world [11,13]. However, seagrass could oxygenate their roots [41] and lose radial oxygen in the rhizosphere of young roots to lower the concentration of sulfide to protect themselves [52]. What is more, sulfur-oxidizing bacteria in this ecosystem may also alleviate the sulfide stress for seagrass by oxidation of sulfide [53]. A higher abundance of genes was found to participate in the process of sulfur oxidation than sulfate reduction in the rhizosphere of the seagrass Z. marina [14]. We also found a high percentage of sulfate respiration (129 records) and respiration of sulfur compounds (131 records) in the FAPROTAX analysis result. This may indicate that microbes in the seagrass rhizosphere also play an important role in the sulfur-related cycle in the coral reef ecosystem.
Bioavailable nitrogen is crucial to all living organisms, but it is still a limiting nutrient globally [54]. The nitrogen enters the ecosystem from the air in the form of ammonia by the microbial nitrogen fixation, which is an essential link of the nitrogen cycle due to nitrogen usually acting as the limiting factor for productivity in the oligotrophic seagrass meadow and coral reef ecosystems [55]. Welsh et al. (2002) found that the microbes capable of sulfate-reducing are the significant component of the diazotrophs in many seagrass ecosystems [56]. Besides, the microbes involved in the nitrogen cycle, as revealed by FAPROTAX analysis in this study, mainly involved in the process of nitrification, aerobic ammonia oxidation, nitrate reduction, nitrate respiration, and nitrogen respiration.
Furthermore, the microbes conducted of nitrification activity mainly came from the genus Nitrosopumilus, Nitrososphaera, and Nitrospira. Nitrification is a process of oxidizing ammonia via nitrite to nitrate, which was assumed as a two-step process catalyzed by chemolithoautotrophic microorganisms before 2015 [57,58]. Daims et al. (2015) have reported that a completely nitrifying bacterium from the genus Nitrospira, which was present in diverse environments, and those findings confirmed that completely nitrifying Nitrospira played important roles in the nitrogen cycle-related microbial functional groups [58]. Although the ammonia available concentrations in most ocean waters are low, this is suitable for the living of comammox organisms. However, no comammox gene has been found in ocean waters until now. To explore microbes capable of comammox, a future research hotspot for environmental microbiologists is underway [54].
The diversity of carbon metabolism found in this study was very high, such as aerobic chemoheterotrophy, chemoheterotrophy, fermentation, aromatic compound degradation, photoautotrophy, methanogenesis, and methylotrophy ( Figure S3). Many microbes of the phylum Planctomycetes were involved in the process of aerobic chemoheterotrophy. Like the genus Blastopirellula, a dominant chemoorganotrophic genus in the Black Sea sediments, are chemoheterotrophic [59,60], and their the major carbon and energy sources are carbohydrates [59]. Eight OTUs were detected in the methylotrophy from the genus Methanomassiliicoccus, unclassified Methylophilaceae, and Methylophaga, which accounted for 0.87% of all detected functional groups. Moreover, the putative methylotrophic bacteria, such as Methylotenera and Methylophaga, were more abundant in healthy seagrasses and could be used as indicators of seagrass health root microbiomes [61]. Besides, the microbes involved in sulfur-cycling, including sulfide-oxidizing (e.g., Candidatus Thiodiazotropha and Candidatus Electrothrix) and sulfate-reducing (e.g., SEEP-SRB1, Desulfomonile, and Desulfonema), were more abundant in stressed seagrass [61]. Hence, there is a need to investigate the relationship between the composition and functions of rhizosphere microbes and seagrass health.

The Core Microbial Community in Seagrass Rhizosphere across the Four Coral Reef Ecosystems
Identification of the core microbial community may provide the cues for understanding the key players in sustaining the growth and health of the seagrass, regardless of the seagrass species and locations. The taxonomy of the predominant core microbial community in this study was Desulfobulbaceae (phylum Proteobacteria), Bacillaceae 1 (phylum Firmicutes), Rhodobacteraceae (phylum Proteobacteria), and Streptococcaceae (phylum Firmicutes). While for seagrass Z. marina, Z. noltii, and Cymodocea nodosa [16], the core seagrass rhizobiome consisted of 0.2% of all OTUs, about 12.81% of all obtained OTUs were identified as core OTUs for the sample investigated in this study. The core microbial communities of different seagrass species and distributing locations may have different community composition and species specialty. The effect of the different environmental factors in the different sampling sites could one of the reasons contributing to this phenomenon [51].
Most of the bacterial core community in the seagrasses rhizosphere was involved in the sulfur cycle [16]. The FAPROTAX analysis from our investigation of the sulfur-related biogeochemical cycle showed that sulfate respiration, dark oxidation of sulfur compounds, and sulfite respiration were dominant in all putative functional groups. For instance, the FAPROTAX analysis showed that genus Desulfocarbo, Desulfatiglans, Desulfosarcina, Desulfobulbus, Desulfopila, Desulfovibrio, Desulfomonile, unclassified Desulfobulbaceae, Desulfovibrio, Desulfomonile, and unclassified Desulfobacteraceae participated in the sulfate respiration process, and 18.76% of the core community was affiliated to above-related genera. Network analysis demonstrated that two module hubs were OTU114 and OTU1807, and they were affiliated to family Rhodobacterales and class Gammaproteobacteria, respectively. Family Rhodobacterales was predicted to play an important role in the process of aerobic chemoheterotrophy and chemoheterotrophy, thus they were involved in the carbon cycle of the seagrass ecosystem. Compared to the neighboring bulk sediment, the bacterial production in the seagrass rhizosphere exhibited a diel pattern that the production rates could be two times higher during daytime than at night [62,63]. In the oligotrophic seagrass meadow, DOC excreted from the seagrass rhizosphere accounted for a large proportion of carbon source for sediment bacteria [64].
In this study, the microbes' carbon cycle was mainly involved in chemoheterotrophy, aerobic chemoheterotrophy, phototrophy, and fermentation.

Conclusions
This investigation mainly focused on the taxonomy structure and functional variations, and core microbial community of seagrass rhizosphere across four coral reef ecosystems. The predominant microbial phyla were Proteobacteria, Bacteroidetes, and Planctomycetes, and the four major functional groups were sulfate respiration, respiration of sulfur compounds, aerobic chemoheterotrophy, and chemoheterotrophy in this study. In the aspect of alpha diversity, significant differences existed based on seagrass species and sampling locations, while no such variations were detected for beta diversity between different sampling locations. The investigated seagrass rhizosphere microbial communities demonstrated significant community composition variations at the genus level, and functional groups also differed among different seagrass species (four seagrass species in the Xisha Islands) and the same seagrass species (seagrass H. ovalis and T. hemprichii) from different locations. The core microbial community of all microbial communities was identified, and most of the microbes had a positive correlation (82.41%). Two module hubs, OTU114 and OTU1807, affiliated to family Rhodobacterales and class Gammaproteobacteria, respectively, were identified as the keystone species. In addition, TOC, ammonium, and nitrate were the significant environmental factors correlated with the core microbial community structure. This study will provide new insight into the seagrass rhizosphere microbiome of coral reef ecosystems and will contribute to more effective management for ecological conservation and restoration policymaking for coral reef ecosystems, such as supplying the scientific data for directed isolation of functional bacteria belonging to the core microbe or species-specific seagrass rhizosphere, thus allowing for characterization of the isolated strain's functional ability and measurement of their seagrass-promoting traits by bacteria inoculant.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/life11080852/s1, Figure S1: Rarefaction curves of OUT diversity in each sample; Figure S2. Phylogenetic and taxonomic alpha diversity of seagrass rhizosphere microbial communities (A: phylogenetic diversity; B: richness; C: Shannon; D: Simpson); Figure S3. The bubble plot of functional groups in all rhizosphere microbial communities of seagrass; Figure S4. Zi-Pi plot of network nodes of core microbial community shared by all samples; Table S1. The sequences characteristics of the all obtained seagrass rhizosphere samples; Table S2. The significance of differences in taxonomy structure of microbial community between different seagrass species (A) and seagrass T. hemprichii (B) and H. ovalis (C) from different sampling sites on the genus level; Table S3. The significance of differences in taxonomy structure of microbial community between seagrass H. ovalis from different sampling sites on the genus level; Table S4. The significance of differences in taxonomy structure of microbial community between seagrass T. hemprichii from different sampling sites on the genus level; Table S5. The significance of differences in functional groups of microbial community between different seagrass species on the genus level; Table S6. The significance of differences in functional groups of microbial community between seagrass T.hemprichii from different sampling sites on the genus level; Table S7. The significance of differences in functional groups of microbial community between different seagrass H.ovalis from different sampling sites on the genus level.  Data Availability Statement: All the raw sequences obtained from this study were deposited in the NCBI sequence read archive (SRA) under accession number PRJNA497291.