16S rRNA gene amplicon-based metagenomic analysis of bacterial communities in the rhizospheres of selected mangrove species from Mida Creek and Gazi Bay, Kenya

Prokaryotic communities play key roles in biogeochemical transformation and cycling of nutrients in the productive mangrove ecosystem. In this study, the vertical distribution of rhizosphere bacteria was evaluated by profiling the bacterial diversity and community structure in the rhizospheres of four mangrove species (Sonneratia alba, Rhizophora mucronata, Ceriops tagal and Avicennia marina) from Mida Creek and Gazi Bay, Kenya, using DNA-metabarcoding. Alpha diversity was not significantly different between sites, but, significantly higher in the rhizospheres of S. alba and R. mucronata in Gazi Bay than in Mida Creek. Chemical parameters of the mangrove sediments significantly correlated inversely with alpha diversity metrics. The bacterial community structure was significantly differentiated by geographical location, mangrove species and sampling depth, however, differences in mangrove species and sediment chemical parameters explained more the variation in bacterial community structure. Proteobacteria (mainly Deltaproteobacteria and Gammaproteobacteria) was the dominant phylum while the families Desulfobacteraceae, Pirellulaceae and Syntrophobacteraceae were dominant in both study sites and across all mangrove species. Constrained redundancy analysis indicated that calcium, potassium, magnesium, electrical conductivity, pH, nitrogen, sodium, carbon and salinity contributed significantly to the species–environment relationship. Predicted functional profiling using PICRUSt2 revealed that pathways for sulfur and carbon metabolism were significantly enriched in Gazi Bay than Mida Creek. Overall, the results indicate that bacterial community composition and their potential function are influenced by mangrove species and a fluctuating influx of nutrients in the mangrove ecosystems of Gazi Bay and Mida Creek.

The mangrove ecosystems provide important ecological and economic functions, including protecting coastlines from storm damage and erosion, degrading environmental contaminants, and providing nursery habitats for numerous aquatic organisms [4]. In addition, they form an enormous food web, providing a myriad of microorganisms with nutrients [5]. The regular tidal variations, pH, temperature, salinity, light, rainfall and nutrient availability create a dynamic environment within the mangrove ecosystems [6], which provides a niche for a wide range of organisms. Thus, the complex mangrove ecosystems have generated increased interest among microbial ecologists [7]. Complex interactions of the microbes in plant ecosystems maintain the harmony of different biogeochemical processes and sustain the nutritional status and ecological balance [8]. Special environmental conditions such as high salinity, high sulfur and low oxygen, drives the proliferation of unique microorganisms that contribute to sulfate reduction, methane cycling and ammonia oxidation within mangrove sediments [9]. Therefore, it is important to understand the bacterial structure and species composition underlining mangrove sediments, particularly within the rhizospheres of specific mangrove species.
A recent study by Wu et al. [10] on the vertical distributions of rhizosphere bacteria from three mangrove species in Beilun Estuary, South China showed that bacterial communities from mangrove rhizosphere sediments were dominated by Proteobacteria (mostly Delta-proteobacteria and Gamma-proteobacteria) and to some extent Chloroflexi, Bacteroidetes, Planctomycetes and Acidobacteria. Furthermore, their study revealed that plant species and depth played a role in shaping the microbial communities, however, the influence of deterministic factors (environmental parameters) such as pollution and mangrove encroachment were not considered. Considering the important ecological functions of mangrove microbiomes and the fact that this unique ecosystem lies at the interface between land and ocean, investigating the deterministic factors that shape microbial communities in line with increasing climate change is crucial.
Previous studies on mangroves in Kenya have concentrated on floristic composition and distribution of mangrove species, economic utilization and regeneration strategies of the principal species [11]. However, data on microbial community diversity is limited due to low attention in exploring the mangrove habitats for microbial diversity [12,13], an information gap that this study has also addressed. The few microbial diversity studies that have been conducted were mainly culture-dependent [12]. Considering that about 98% of microbes are uncultivable, the 'omic' technologies are emerging as potential tools to ensure holistic insight into environmental systems [14] and this has led to the increasing application of high-throughput sequencing in studying microbial communities [10,15]. In this study, the vertical profiles of bacteria communities in the rhizospheres of four mangrove species namely, Sonneratia alba, Rhizophora mucronata, Ceriops tagal and Avicennia marina, were determined and compared. Since we anticipated some possible contamination at Mida Creek attributed to humaninduced activities such as firewood harvesting, pollution from plastics and faeces, pollution from oil spills, overharvesting for building materials and encroachment for settlements [16,17], we obtained our samples from two distinct mangrove ecosystems; viz, a pristine site in Gazi Bay, located in the South Coast of Kenya and the polluted site in Mida Creek, located in the North Coast of Kenya. We successfully applied for the first-time in-depth analysis of the bacterial communities within the rhizospheres of four mangrove species from Kenya, using high-throughput sequencing of the V3-V4 regions of the 16S rRNA gene on Illumina MiSeq platform.

Ethical statement
The National Commission for Science, Technology and Innovation of Kenya (NACOSTI) approved this research study, National Environmental Management Authority of Kenya (NEMA) provided the access permit (for field sampling), Kenya Wild life Services (KWS) and Kenya Plant Health Inspectorate Services (KEPHIS) provided permits that facilitated the shipment of samples to Laval University, Canada. The field studies neither involved endangered nor protected species.

Study site
Two study sites (Mida Creek and Gazi Bay) in Kenya were investigated. Mida Creek is located in Kilifi County (03˚34 0 S, 039˚96 0 E), about 88 km North of Mombasa and approximately 25 km south of Malindi town in a planigraphic area of 32 km 2 [18]. It is characterized by a hot and humid tropical climate with an average annual temperature of 27˚C. Humidity is high throughout the year, up to 90% relative humidity during the rainy season [18]. In addition, the creek is affected by anthropogenic degradation such as overharvesting of mangroves for firewood, timber and fish traps, pollution from plastics, faeces and oil spills, clearing of mangrove and conversion to other land uses such as aquaculture, urban development and tourism [16,17].
Gazi Bay is located in Kwale County (04˚44 0 S 039˚51 0 E), approximately 55 Km from Mombasa, South Coast of Kenya. The Bay is sheltered from strong waves by the presence of the Chale peninsula to the East and a fringing coral reef to the South. The climate is hot and humid. The average annual temperature and humidity are about 28˚C and upto 95%, respectively [18]. The mangrove forests in both sites display a similar zonation and mangrove species contribution among the four dominant species A. marina (30%), R. Mucronata (25%), S. alba (15%) and C. tagal (10%), which contribute about 80% of the total mangrove formation in both forests. Sonneratia alba (about 6-10 m tall) forms the outermost zone (seaward side) towards the open water followed by pure stands of Rhizophora mucronata (about 8-12 m tall) or mixed stands of Rhizophora mucronata and Bruguiera gymnorrhiza (about 10-20 m tall) and in turn these stands are followed by pure stands of Ceriops tagal (about 3-5 m tall) and Avicennia marina (about 12-18 m tall) along the creek [19][20][21]. Rhizophora mucronata have well-developed prop roots that accumulate large stocks of debris, perhaps contributing to some accretion that supports the extensive tidal flats seen in the area [21]. in a dry iced box before they were transported and stored at -20˚C before further analyses was performed at the Laval University, Canada.

Nutrient analysis of sediment samples
Nutrient analyses of soil samples for nitrogen, carbon, phosphorus, potassium, calcium, magnesium and sodium were conducted according to standard methods [23]. Determination of pH and electrical conductivity was done using the calcium chloride method at a ratio of 1:2 using a digital pH meter [Corning pH meter 140] [24] and electrical conductivity meter [Conductivity meter type CDM 2d radiometer Copenagen] [25], respectively.

Total community DNA extraction, PCR protocol and Illumina MiSeq sequencing
The excised roots containing rhizosphere sediment were placed in a 50 mL sterile falcon tube containing autoclaved phosphate buffer and shaken for 2 min to release the rhizosphere sediment from the surface of the roots. Total genomic DNA was extracted directly from 0.25 g of the sediment using Power Soil DNA isolation kit (DNeasy PowerSoil Kit, Qiagen, Germany) in accordance with the manufacturer's protocol. The extracted DNA was quantified using NanoDrop 1000 spectrophotometer (Thermo Scientific, USA). For each PCR reaction, 20 ng of total DNA was used as a template. Briefly, amplification of the bacterial 16S rRNA gene fragment (V3-V4 region) was performed using the sequence specific regions described in [26]. The following generic oligonucleotide sequences were used for amplification: Bakt_341F-long AATGATACGGCGACCACCGAGATCTACAC[index1]TCGTCGGCAGCGTCAGATGTGTATAA GAGACAGCCTACGGGNGGCWGCAG and Bakt_805R-long CAAGCAGAAGACGGCATACGAGAT [index2]GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC. The primers used in this work contain Illumina specific sequences protected by intellectual property (Oligonucleotide sequences © 2007-2013 Illumina, Inc., all rights reserved). The PCR was carried out in a total volume of 50 μL that contains 1X Q5 buffer (NEB), 0.25 μM of each primer, 200 μM of each dNTPs, 1 U of Q5 High-Fidelity DNA polymerase (NEB) and 1 μL of template DNA. The PCR started with an initial denaturation at 98˚C for 30s followed by 10 cycles of denaturation at 98˚C for 10s, annealing at 55˚C for 10s, extension at 72˚C for 30s. An additional 25 cycles were as follows; denaturation at 98˚C for 10s, annealing at 65˚C for 10s, extension at 72˚C for 30s, followed by a final extension step at 72˚C for 2 min. The PCR amplicons was purified using the Axygen PCR cleanup kit (Axygen). Quality of the purified PCR product were checked on a DNA7500 BioAnalyzer chip (Agilent) and quantified using a Nanodrop 1000 (Thermo Fisher Scientific). Barcoded Amplicons were pooled in equimolar concentration and sequenced on the Illumina Miseq (paired-end 300 bases with two index reads) at the Institute for Integrative System Biology (IBIS)/ (Laval University's Genomic Analysis Platform, Quebec, Canada). Sequences from this study have been deposited into the SRA under the BioProject ID PRJNA644929.
Predictive functional analysis of bacterial communities was done using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States 2 (PICRUSt2). To predict functions, representative sequences were first aligned to HMMER [34]. The alignment was then placed into a reference tree using EPA-NG [35] and gappa [36]. Multiple 16S rRNA gene copies were then normalized and gene families predicted using Castor-a hidden state prediction tool [37]. The predicted gene families were subsequently collapsed into MetaCyc pathways using MinPath [38].
Statistical analyses. Unless otherwise stated, statistical analyses were performed using R v3.6.1 [33]. A three-factor (sites, mangrove species and depth differences) test of differences in α-diversity and physicochemical parameters was done by the non-parametric Kruskal-Wallis H test using the agricolae [39] R package. Post hoc test for mean separations was based on Fisher's least significant difference. Comparison of microbial community differences (β-diversity) in the two study sites, among mangrove species and between sampling depths were based on Bray-Curtis dissimilarities and principal coordinates analysis (PCoA) using ampvis2 and ape [40] R packages. Permutational multivariate analysis of variance (PERMANOVA) was used to test the significance of the differences in multivariate space using the vegan package. Furthermore, a post-hoc test on significant PERMANOVA (q � 0.05) was performed using the "adonis_pairwise" function in the metagMisc R package. The influence of physicochemical parameters on community-level differentiation (community structure) was determined by constrained redundancy analysis (RDA). RDA was performed on Hellinger-transformed bacterial and environmental data using a combination of ampvis2 and vegan R package functions. The significance of the constraining variables was determined based on permutation test using the vegan function 'anova.cca'. Further investigation on the contribution of environmental factors on the microbial community differentiation was achieved by variance partitioning (based on chi-square) using the vegan R package.
Differential abundance (Benjamini-Hochberg false discovery rate (FDR)-adjusted, p�0.05) of OTUs in both sites of study and among mangrove species was determined using aldex2 [41] R package. Linear Discriminant Analysis Effect Size (LEfSe) [42] was used to determine the presence of biomarker PICRUSt2 predicted KEGG orthologs (enzymes) by applying the Kruskal-Wallis alpha significance threshold of � 0.01 and an LDA (linear discriminant analysis) score of 3.0. The GraPhlan software [43] was subsequently used to visualize the detected biomarkers. Pathways significant differentiation was determined by also applying the Kruskal-Wallis alpha significance threshold of � 0.05 and FDR-adjusted Fisher's least significant difference post hoc analysis using the agricolae R package.

Sequencing statistics, diversity measures and richness estimates
A total of 3,038,718 sequence reads were obtained from the 64 mangrove rhizosphere samples, out of which 189,254 high-quality 16S rRNA sequence reads were clustered into 4,295 OTUs. The rarefication of counts to 1,126 reads per-sample was sufficient for explaining differences in bacterial diversity (S1 Fig). Approximately 56% of the OTUs were shared between Gazi Bay and Mida Creek. Also, mangrove plant species comparison revealed that A. marina and C. tagal had the most OTUs (1341 and 1332, respectively). Sonneratia alba comprised 978 OTUs while R. mucronata had the least number of OTUs (901). Comparison of shared OTUs in the two study sites revealed that 24.3% of the core OTUs (OTUs found in � 20% of the samples from both sites) were shared between the rhizospheres of A. marina in Gazi Bay and Mida Creek (S2A Fig). For S. alba, 22.6% of the core OTUs were shared between the two sites while 9.8% and 4.6% of the core OTUs were respectively shared in the rhizosphere of R. mucronata and C. tagal for Gazi Bay and Mida Creek. Overall, 12.8% of the core OTUs were shared between the two sites (S2B Fig).
Bacterial diversity in the rhizosphere of the mangrove plant species revealed significant differences on comparison of the two sites investigated (Table 1). Shannon-Wiener, Observed and Chao1 diversity indices were all significantly higher (Fisher's LSD, p�0.01) in the rhizosphere of S. alba in Giza bay compared to Mida Creek. Other rhizosphere species comparison in Giza bay and Mida Creek were not significantly different apart from Shannon-Wiener diversity for A. marina, which was found to be significantly higher (p � 0.05) in Mida Creek (Table 1) Spearman's rank correlation of α-diversity measures to physicochemical parameters revealed significant correlations. 24, p = 0.05) and K (r = -0.34, p = 0.006) all significantly correlated negatively with Observed species, while salinity, pH, Na, EC and phosphorus did not correlate significantly. For Shannon-Wiener measure of diversity, all the physicochemical parameters showed significant negative correlation except for pH, which did not correlate significantly (S1 Table). The correlation results of Chao1 to physicochemical properties was similar to the Observed as both matrices measures specie richness. For Chao1, Ca (r = -0.31, p = 0.01), C (r = -0.37, p = 0.002), Mg (r = -0.38, p = 0.001) and K (r = -0.33, p = 0.007) were all negatively correlated to species richness, while salinity, pH, Na, EC and phosphorus did not correlate significantly. Overall, only pH failed to correlate with any of the α-diversity matrices.

Bacterial composition and differences across sites and mangrove species
The bacterial composition of Gazi Bay and Mida Creek was found to be highly diverse, comprising a total of 25 bacterial phyla with at least 1% abundance in the samples. Proteobacteria (41%) was the most abundant phylum in both sites. This was followed by Planctomycetes

PLOS ONE
Desulfobulbaceae, Desulfarculaceae, Ectothiorhodospiraceae and Thiotrichaceae were differentially enriched (p � 0.05). Differentially abundant families in the rhizosphere of S. alba were Flavobacteriaceae, Moduliflexaceae, Spirochaetaceae, Kineosporiaceae and Chthoniobacteraceae while Desulfobacteraceae was enriched (p�0.05) in the rhizospheres of both R. mucronata and S. alba. Overall, there were no significant differences (p > 0.05) in bacterial families between these two-mangrove species.
Several species were not classified (35.06%) at the genus level using the SILVA reference database. Among the top 10 classified bacterial genera across all sites and plant species were Sulfurimonas > Sva0081 sediment group > Pir4 lineage > Arcobacter > Woeseia > Blastopirellula > Spirochaeta > wb1-A12 > Psychrilyobacter > Desulfatitalea, according to sequence of abundance ( Fig 1B). Investigation of the bacterial genera for differential abundance between mangrove species in Gazi Bay and those of Mida Creek revealed significant differences. Nitrospira  Table).

Bacterial community differentiation based on site, depth and mangrove species
Beta-diversity based on Bray-Curtis distance revealed a significant clustering of samples within the ordination space. Samples clustered mainly based on mangrove species, and much less, based on the site differences and on depth of sampling (Fig 2). The differences observed in multivariate space for mangrove species (PERMANOVA R 2 = 16%, p = 0.001), site of sampling (PERMANOVA R 2 = 4.5%, p = 0.001) and depth of sampling (PERMANOVA R 2 = 4.0%, p = 0.02) were all significant. Pair-wise post hoc analysis of PERMANOVA revealed that the bacterial community composition and structure are significantly different between the mangrove plant species (FDR-adjusted p = 0.001) as well as between the sites (p = 0.001) and sampling depth (p = 0.017). The bacterial community composition and structural differences in the mangrove species rhizospheres is presented in details in S3 and S4 Tables.

Site and mangrove species differences in physicochemical parameters
Site-based comparison of the different mangrove trees were mostly significant, apart from salinity and electrical conductivity that were not influenced by geographical distance for S. alba and A. marina (Table 2). pH and calcium were significantly higher (Kruskal-Wallis, p � 0.05) in all mangrove plant species of Mida Creek compared to Giza bay. Apart from the rhizosphere of R. mucronata, all other mangrove species rhizosphere in Giza bay had significantly lower (p < 0.05) physicochemical properties compared to the species in Mida Creek. The physicochemical parameters that were higher in the rhizosphere of R. mucronata in Gazi Bay included potassium, sodium, phosphorus, total carbon, nitrogen, salinity and electrical conductivity.

Influence of environmental factors on bacterial composition and structure
Redundancy analysis led to the explanation of 91% of the observed variation in β-diversity out of which 18% was explained by the constrained relation of site differences, mangrove species differences and depth of sampling, while 72% was explained by the unconstrained environmental variables (Fig 3A). The total explained variation of the microbial community composition of the first two axes (RDA1 and RDA2) reached 51%. The constrained variation resulting from the differences in site, species and depth were 6.5% and 3.7% for RDA1 and RDA2, respectively. The subjection of RDA to statistical test revealed that the observed differentiation was significant (p = 0.001). On the overall, the significance of environmental terms fitted into the RDA model indicated that calcium, potassium, magnesium, electrical conductivity, pH and salinity with significance values of p = 0.001 and carbon with significance of p = 0.002 contributed the most to the species-environment relationship ( Table 3). Other contributing physicochemical factors were nitrogen (p = 0.007) and sodium (p = 0.01), while phosphorus (p = 0.07) was the only physicochemical factor without a significant influence on the bacterial community. Though pH influenced the bacterial community across sites, its effect also separated the samples

PLOS ONE
from the different sites (Gazi Bay and Mida Creek) along the y-axis (Fig 3A). Desulfobacteraceae and Syntrophobacteraceae were projected in the positive direction of RDA1 and thus observed to be positively correlated with carbon, nitrogen, potassium and magnesium, while Ardenticatenales appeared not to be strongly influenced by physicochemical factors.
To further disentangle the respective contribution of plant species, chemical parameters, depth and geographical variation in explaining the bacterial community structure, variance partitioning was performed (Fig 3B). The sediment chemical parameters had the highest influence on the bacterial community compared to mangrove species and site/depth. The sediment chemical parameters alone explained for 6% and in combination with other factors for a total of 13% of the bacterial community pattern followed by mangrove species (alone: 4% while in combination with other factors: 10%) and site/depth (alone: 1% while in combination with other factors: 4%). Eighty-two percent of the differences in bacterial community structure was not explained by a combination of all the measured variables, indicating that other underlying factors not investigated are also influencing the community differences. Spearman's rank correlation of physicochemical parameters to the top bacterial genera revealed significant associations (S5 Table). Pir4 lineage, Desulfobacter, Pleurocapsa PCC.7319, Rhodopirellula, Vibrio, Desulfatiglans and Marinobacterium correlated the most to physicochemical terms.

Predictive functional analysis of bacterial communities
A total of 8819 PICRUSt2 predicted KEGG orthologs (enzymes) were collapsed into 406 Meta-Cyc pathways. Special focus was given to enzymes associated with energy metabolism (nitrogen metabolism, carbon fixation, methane metabolism, sulfur metabolism and photosynthesis) and the biosynthesis of important secondary metabolites. One hundred and   KEGG orthologs associated with carbon metabolism in prokaryotes were differentially enriched in the rhizosphere of R. mucronata in Gazi Bay. Desulfatiglans, Vibrio, Desulfovibrio, Desulfatitalea and Woeseia were the potential phylotypes contributing to carbon metabolism. Investigation of genes for methane metabolism also revealed that R. mucronata associated microbiomes were enriched with KOs for methane metabolism in Gazi Bay. Potential contributing phylotypes included Desulfobacteraceae (Desulfosarcina), Desulfarculaceae (Desulfatiglans) uncultured Dehalococcoidia, Moduliflexaceae, uncultured Crenarchaeota, Syntrophobacteraceae, Thalassospiraceae (Thalassospira), Pirellulaceae, Psychromonadaceae (Psychromonas) and Bogoriellaceae (Georgenia).
Pathway differentiation according to site and mangrove species also indicated that they were significant differences in the abundance of predicted pathways for sulfur, carbon, nitrogen and methane metabolism (Fig 5). It was observed that pathways for both sulfur and carbon metabolism were significantly higher (Fisher's LSD, p<0.05) in Gazi Bay than Mida Creek while other difference in predicted pathways were mangrove species-dependent. Meanwhile, for all mangrove species, depth-based comparison within site and across sites indicated that there were no significant (Fisher's LSD, p>0.05) differences.

Discussion
Mangrove ecosystems are widespread in estuarine and coastal regions of the tropics, providing protection and nutrients to several microbial communities [5]. The dynamic environment within mangrove ecosystems provides an excellent niche for a wide range of organisms. In this study, we determined the bacterial communities of rhizospheric sediments from four mangrove species (Sonneratia alba, Rhizophora mucronata, Ceriops tagal and Avicennia marina), which are common along the Coastline of Kenya, using the Illunima MiSeq sequencing technology. Also, the influence of mangrove species, sampling depth, geographical location and other physicochemical factors, in shaping bacterial communities in the mangrove ecosystem was determined.
An important finding drawn from the comprehensive study of the mangroves in Gazi Bay and Mida Creek is that the different mangrove species and chemical properties had a greater influence on the rhizosphere microbiome diversity compared to sampling depth. For instance, despite the differences in geographical location, A. marina in both Gazi Bay and Mida Creek had the most bacterial diversity and richness (Table 1). This finding is consistent with the study by Wu et al. [10], where they observed that differences in mangrove species (Bruguiera gymnorrhiza, Kandelia candel and Aegiceras cornicu-latum) had a strong pattern on bacterial α-diversity. Besides, it has been demonstrated that plant roots are capable of imposing a selective force on rhizospheric microbiomes in both mangrove swamps [44] and other terrestrial environments [45], and that this phenomenon is usually plant species dependent. As earlier stated, another important variable influencing α-diversity were the sediment chemical parameters, and this finding is consistent with previous studies of bacterial communities in mangrove ecosystems [46,47]. In this study, most of the investigated physicochemical parameters correlated inversely with both species diversity and richness matrices (S1 Table); suggesting that an influx of nutrients from anthropogenic sources reduces bacterial diversity and the influence of mangrove species on mangrove sediment microbiomes. Notable examples were the significantly higher diversity and species richness among some mangrove species in Gazi Bay (considered pristine) compared to Mida Creek where nutrient inputs were observed to be significantly higher. Also, the low percentage of shared OTUs between same mangrove species in the two study sights (S2 Fig) highlights the importance of deterministic factors in shaping mangrove sediment microbiomes. Overall, this pattern of bacterial response indicates that in Gazi Bay and Mida Creek, the bacterial alpha diversity pattern is determined by an interaction of both mangrove species and physicochemical factors than with sampling depth.
Though differences in sampling depth did not have a strong pattern on α-diversity, investigation of bacterial community composition revealed it significantly contributed in shaping the

PLOS ONE
bacterial community (Fig 2). The clustering of samples in the PCoA plot according to mangrove species, geographical location and depth of sampling were all significant and respectively explained 16%, 4.5% and 4% of the total variation in bacterial community composition (S3 Table). In an exploratory study of the prokaryotic communities in mangrove sediments of southern China, Zhang et al. [47] indicated that there were significant interactions between sampling site and plant species that shaped the bacterial community. Similarly, Wu et al. [10] demonstrated that differences in mangrove plant species imposed a strong selective pattern on the rhizosphere microbial communities compared to depth, which did not show a clear clustering pattern, while further evidence on the influence of geographical location on mangrove microbial community composition has been presented in recent studies [48,49].
RDA revealed a strong bacterial species-environmental variable relationship that influenced the community composition. Also, variable partitioning revealed the importance of the inter-relationship among environmental properties and between bacterial species (Fig 3B). In this study, Calcium, potassium, magnesium, electrical conductivity, pH and salinity influenced the bacteria community composition the most (Table 3). These environmental variables have been reported as factors that influence bacterial community composition and structure [15,47,50]. Zhang et al. [47] reported that total organic carbon and mean annual precipitation were the main environmental variables influencing the bacterial community composition of mangrove sediments in Southeastern China.
Although evidence abounds on the effect of sediment chemical parameters, mangrove species and sampling depth on bacterial communities, it is important to state that diversity and assembly patterns for this study also followed known intertidal zonation patterns. For example, alpha diversity reduced in the seawards tidal zone compared to the upper intertidal zones. Also, regardless of site differences, samples drawn from the lower intertidal zones appeared to cluster together in the PCoA ordination space. This observation is consistent with previous studies [51,52] and thus highlights the importance of mangrove forest intertidal zonation on bacterial alpha diversity.
The phylum Proteobacteria was the most abundant in the rhizosphere of all mangrove species in both sites. Also, a total of 25 phyla had relative abundance above 1% (Fig 1A). This observed high number of bacterial phylotypes is one of the trademark features used in distinguishing the microbial communities of mangrove sediments from other biomes [47]. The dominant phyla from this study have been consistently detected in similar studies that evaluated mangrove rhizosphere bacterial communities [48,[53][54][55]. The high abundance of Proteobacteria can be traced to the dominant physicochemical parameters in mangrove sediments, particularly the nitrogen deficiency and the usually rich sulfur and carbon content. Dominant proteobacteria classes such as Delta-proteobacteria are reported to be metabolically versatile, and have a strong affinity to saline environments [56]. Besides, they actively participate in organic matter decomposition, ammonia oxidation and sulfate reduction [47]. The phylum Planctomycetes has been reported to be a component of microbial core in the mangrove ecosystems; taking part in methane and sulphur metabolism [57]. Also, Cyanobacteria, Chloroflexi and Planctomycetes are reportedly involved in nitrogen cycling in mangrove ecosystems [58]. Related studies by Yamada et al. [59] and Ghosh and Bhadury [60] also demonstrates additional importance of Chloroflexi in the decomposition of organic matter in mangrove ecosystems.
The most abundant genus, Sulfurimonas (Fig 1B) has been detected and recovered mostly from sulfur rich marine environments [7,52]. Members of the genus Sulfuriomonas are versatile hydrogen and sulfur oxidizing chemolitotrophs [61,62] and have been detected abundantly in both high [50] and low [52] tidal sediments as well as in hydrothermal vents [61]. Other abundant genera including Sva0081 sediment group (Desulfobacteraceae), Pir4 lineage (Pirellulaceae), Arcobacter, Woeseia, Blastopirellula, Spirochaeta, wb1-A12 (Methylomirabilaceae), Psychrilyobacter and Desulfatitalea are also associated with sulfur rich environments and some tend to co-occur. For instance, Sva0081, Woeseia and Spirochaeta were recently detected among differentially enriched bacterial genera in Tanmen harbor, east coast of Hainan Island, China [63]. Sva0081, Pir4 lineage, Spirochaeta and Desulfatitalea are established sulfur oxidizers [63][64][65] and their abundance in the samples can be linked to the rich sulfur and carbon content of the study sites. Pir4 lineage and Spirochaeta have also been reported to contribute to the decomposition of organic matter in marine environments [65,66]. Besides participating is sulfur oxidation, Desulfatitalea have been reported to also possess the nrfA gene for dissimilatory nitrate reduction to ammonium [64], while both Arcobacter and Desulfatitalea are suspected to play a role in the degradation of microplastics in marine environments [67]. The versatility of these organisms suggests the reason for their high abundance in the sites of study, and underlines their potential importance in pollution mitigation and biogeochemical cycling in mangrove sediments. The differentially abundant phyla and genera detected in the same mangrove species from different sites can be described as species sensitive to changing environmental factors. These sensitive bacterial species can be of importance in evaluating the impact of anthropogenic factors or climate change on mangrove microbiomes and their corresponding ecological function. Overall, the bacterial composition of Mida Creek and Gazi Bay suggests a bacterial interaction that leads to a coupled cycling of sulfur, carbon and nitrogen.
The prediction of bacterial functional profiles focused on key enzymes and pathways associated with energy metabolism and the biosynthesis of important secondary metabolites (Fig 4). The results demonstrated the functional potential of the mangrove rhizospheres bacterial communities. Oxyphotobacteria, which was predicted to contribute to photosynthesis are established primary producers in marine environments and are known to also play a role in nitrogen fixation [68]. Also, the functional prediction suggests sulfate oxidizing bacteria were involved in several biogeochemical cycles in the mangrove ecosystem including sulfur, nitrogen, methane and carbon metabolism. Key bacterial family like Syntrophobacteraceae are known to possess genes (Apr/Dsr) for dissimilatory sulfate reduction to H 2 S and are usually found to be among dominant families in sulfate rich sediments [69]. Blastopirellula have been implicated in anaerobic oxidation of ammonium [70]. Also, Marinobacterium is an established nitrogen fixing bacteria [71] while the family Desulfobacteraceae are versatile sulfate reducers [72] capable of coupled nitrogen metabolism [54,73] and hydrocarbon degradation in marine environments [72,74].
Diverse members of bacterial communities contributed differentially to carbon fixation among mangrove species. Also, predicted pathways for both nitrogen, sulfur and carbon metabolism was significantly higher in Gazi Bay compared to Mida Creek (Fig 5). This finding suggests that the inverse association of bacterial diversity and richness with the influx of nutrients, may have impacted the functional potential of the bacterial communities in Mida Creek. Further insights on genes associated with carbon metabolism revealed differential enrichment in R. mucronata rhizosphere. The preferential abundance of genes for carbon metabolism in mangrove sediments has been reported previously [54,75]. Also, Zhao et at. [75] reported that total organic carbon in marine environment significantly affects the distribution of genes for carbon metabolism. For this present study, we consider that the characteristic feature of the attachment of debris to the long roots of R. mucronata may have influenced the bacterial community response to carbon metabolism. Overall, the site and mangrove species-based differences in predicted enzymes for energy metabolism is more likely to be a function of differences in co-interacting physicochemical parameters and critical mangrove species requirements. Though the functional prediction gives a first overview of the potential contribution of microbial communities to energy metabolism and biosynthesis of biomolecules within mangrove rhizospheres, it is recommended that a multi-omics (metagenomics, metatranscriptomics, metabolomics and proteomics) approach, which eliminates known biases associated with PICRUSt2 prediction be applied in future explorations.

Conclusion
This study has generated baseline data of bacterial diversity associated with rhizospheres of specific mangrove species from Kenya. The study revealed that mangrove species and interacting physicochemical terms imposes a strong pattern on the bacterial community structure and diversity of Mida Creek and Gazi Bay. The influx of nutrients into the mangrove environment inversely correlated with bacterial diversity and richness. PICRUSt2-based inference of the potential functional impact of such negative association suggested that critical pathways for sulfur and carbon metabolism could be potentially affected. However, the wide distribution of predicted genes associated with energy metabolism and the synthesis of important secondary metabolites suggest that the investigated mangroves are potentially rich hubs for the recovery of novel bacterial strains / products with wider biotechnological applications.