Effects of tidal influence on the structure and function of prokaryotic communities in the sediments of a pristine Brazilian mangrove

Mangrove forests are ecosystems that constitute a large portion of the world's coastline and span tidal zones below, between, and above the waterline, and the ecosystem as a whole is defined by the health of these tidal microhabitats. However, we are only beginning to understand tidal-zone microbial biodiversity and the role of these microbiomes in nutrient cycling. While extensive research has characterized microbiomes in pristine vs. anthropogenically impacted mangroves, these have, largely, overlooked differences in tidal microhabitats (sublittoral, intertidal, and supralittoral). Unfortunately, the small number of studies that have sought to characterize mangrove tidal zones have occurred in impacted biomes, making interpretation of the results difficult. Here, we characterized prokaryotic populations and their involvement in nutrient cycling across the tidal zones of a pristine mangrove within a Brazilian Environmental Protection Area of the Atlantic Forest. We hypothesized that the tidal zones in pristine mangroves are distinct microhabitats, which we defined as distinct regions that present spatial variations in the water regime and other environmental factors, and as such, these are composed of different prokaryotic communities with distinct functional profiles. Samples were collected in triplicate from zones below, between, and above the tidal waterline. Using 16S ribosomal RNA (rRNA) gene amplicon sequencing, we found distinct prokaryotic communities with significantly diverse nutrient-cycling functions, as well as specific taxa with varying contributions to functional abundances between zones. Where previous research from anthropogenically impacted mangroves found the intertidal zone to have high prokaryotic diversity and be functionally enriched in nitrogen cycling, we find that the intertidal zone from pristine mangroves has the lowest diversity and no functional enrichment, relative to the other tidal zones. The main bacterial phyla in all samples were Firmicutes, Proteobacteria, and Chloroflexi while the main archaeal phyla were Crenarchaeota and Thaumarchaeota. Our results differ slightly from other studies where Proteobacteria is the main phyla in mangrove sediments and Firmicutes makes up only a small percentage of the communities. Salinity and organic matter were the most relevant environmental factors influencing these communities. Bacillaceae was the most abundant family at each tidal zone and showed potential to drive a large proportion of the cycling of carbon, nitrogen, phosphorus, and sulfur. Our findings suggest that some aspects of mangrove tidal zonation may be compromised by human activity, especially in the intertidal zone.

1 Introduction

Mangrove ecosystems exist at the interface of land and sea and constitute a large portion of tropical and subtropical coastlines, spanning 118 countries and approximately 137 000 km² (Giri et al., 2011). Within these ecosystems, the recycling of elements occurs via tightly coupled exchanges between plants, microorganisms, and microbes, with Bacteria and Archaea playing important roles in the biogeochemical cycling of carbon, nitrogen, and phosphate (Imchen et al., 2017; Lin et al., 2019; Alongi, 1988; Reis et al., 2017; Holguin et al., 2001). Mangroves span tidal zones that are characterized by periodic tidal flooding, such that environmental conditions such as salinity vary greatly across small spatiotemporal scales. We refer to the differences in microhabitats and biodiversity in mangrove sediments as zonation, and we consider the mangroves to be divided into three main tidal zones, where the sublittoral zone is always immersed, while the intertidal zone is exposed to daily variations in water content and the supralittoral zone is the region that presents mangrove vegetation but normally is above the sea level. With each tidal cycle the levels of nutrients, oxygen, and salinity fluctuate, resulting in frequent anaerobic conditions and a wide range of redox potentials in the sediments (Holguin et al., 2001; Andreote et al., 2012; Marcos et al., 2018; Lin et al., 2019). While there has been little research on microbial diversity within mangrove tidal zones, research in non-mangrove tidal ecosystems has found distinct microbiomes with significantly different biogeochemistries (Musat et al., 2006; Boehm et al., 2014; Azzoni et al., 2015). Understanding the interrelationship of mangrove tidal zones, their microbial populations, and their role in nutrient cycling is a key component in understanding mangrove adaptation and conservation (Coldren et al., 2019; Saintilan et al., 2020; Shiau and Chiu, 2020; Allard et al., 2020).

Numerous studies have sought to characterize mangrove microbiomes in both pristine and anthropogenically impacted areas (see Supplement file 2). This has led to an increased understanding of the sensitivity of mangrove microbiomes to perturbations, such as pollution (Peixoto et al., 2011; Andreote et al., 2012; Nogueira et al., 2015), environmental degradation (Marcial Gomes et al., 2008; Tong et al., 2019; Cabral et al., 2018), and rising sea levels (Ceccon et al., 2019; Yun et al., 2017; Saintilan et al., 2020). Recently, additional work has sought to understand the diversity of microbial populations across microhabitats defined by distinct biotic and abiotic characteristics (Luis et al., 2019; Lv et al., 2016; Zhang et al., 2018; Rocha et al., 2016). While the particular effects of these perturbations vary from site to site, broadly these result in significantly altered population structures with substantially different functional profiles. This is important to consider since research in mangrove tidal-zone microhabitats (Zhang et al., 2018; Jiang et al., 2013; Zhou et al., 2017) has so far been carried out in anthropogenically impacted sites, thus confounding the interpretation of these findings.

In Brazil, mangrove forests are primarily composed of the genera Rhizophora, Avicennia, Laguncularia, and Conocarpus (Pupin and Nahas, 2014) and are part of the Atlantic Forest, one of the most biodiverse biomes on the planet. Unfortunately, the Atlantic Forest and associated ecosystems are highly threatened by anthropogenic disturbances, resulting in a severe decline from the forest's approximate pre-colonial area (Brooks et al., 1999; Ditt et al., 2013). However, in the southern part of Bahia State, Brazil, a significant fragment of the Atlantic Forest remains preserved within the Environmental Protection Area of Pratigi (MMA, Ministério do Meio Ambiente, 2004). Recent studies on the area have shown that preservation efforts have been generally effective, resulting in a high environmental quality relative to most mangrove forests, both in Brazil and globally (Ditt et al., 2013; Lopes, 2011; Mascarenhas et al., 2019). This preserved area constitutes an increasingly rare site for the understanding of the ecology of unimpacted mangrove forests.

In this study we aimed to characterize the prokaryotic microbiota from sediments of three tidal zones in a pristine mangrove located in the Serinhaém estuary, within the Environmental Protection Area of Pratigi, using 16S ribosomal RNA (rRNA) gene amplicon sequencing. Considering that mangrove zonation is driven, primarily, by tidal variation, we hypothesized that sediments of different tidal zones would differ significantly in the diversity, composition, and functioning of prokaryotic communities as a result of differences in the physicochemical properties affected by tides. We also hypothesized that the intertidal zone, being a highly dynamic environment, would have the highest diversity as reported previously (Zhang et al., 2018; Brose and Hillebrand, 2016). Several important ecological consequences follow on from high levels of biodiversity, such as increased functional redundancy and increased resilience to environmental perturbations (Girvan et al., 2005). However, unlike the findings from impacted mangroves (Zhang et al., 2018; Jiang et al., 2013) we found the intertidal zone to be the least diverse, suggesting that it may, in fact, be more sensitive to anthropogenic disruption. Our study provides insight into the role of microbes in the functioning of mangrove forests and establishes a baseline for monitoring the health of this important ecosystem, since this information is still scarce for pristine mangrove sites.

This work was conducted before a massive oil spill occurred off the coastline of Brazil in August 2019, impacting hundreds of miles of coastline including the Serinhaém estuary, where this research was conducted. This work therefore serves as a baseline measure of the prokaryotic communities of the tidal zones of what was a pristine mangrove forest. We hope that this will spur subsequent research into the effects that anthropogenic interferences have on mangrove ecosystems.

Figure 1 Map and schematic of sediment sampling sites. Here we show the locations of the sampling sites relative to Brazil (a) and Bahia (b). A picture shows the relation of the three sampling sites within each zone (1 sublittoral, 2 intertidal, 3 supralittoral; c). A schematic shows the topographic and tidal relation of each sampling site (d).

2 Materials and methods

2.1 Study area

The Serinhaém estuary is located in the Low South Region of Bahia State, Brazil (Fig. 1), between the coordinates 13^∘35^′ and 14^∘10^′ S, 39^∘40^′ and 38^∘50^′ W. The estuary is within the Environmental Protection Area of Pratigi, an Atlantic Forest region with a total area of 85 686 ha, enclosing a 32 km long portion of the Juliana River and emptying directly into Camamu Bay (Corrêa-Gomes et al., 2005).

2.2 Sampling and DNA extraction

For clarity, here we refer to the location of a sample as a “site” and the collection of sites within a tidal zone as a “zone”. Samples were collected from three tidal zones (centered around 13^∘42^′59.0^′′ S, 39^∘01^′35.9^′′ W) in the Serinhaém estuary in July 2018 during the morning low-tide period. No sites exhibited signs of anthropogenic disturbance or pollution. The three collection zones were chosen based on tidal influence: sublittoral, intertidal, and supralittoral regions (Fig. 1). We evaluated tidal zones based on water level at morning low tide; additional factors were also considered to differentiate between the intertidal and supralittoral zones, such as the presence of surface soil moisture, proximity to restinga vegetation, and the direction of a local guide. From each tidal zone, three samples of superficial sediments (top 10 cm of the surface layer) were collected with a cylindrical sediment core sampler. Precautions were taken to avoid the disruption of rhizospheres associated with vegetation, which could otherwise lead to contamination of the soil by rhizosphere-specific microbes. Sample sites were located at a minimum of 15 m from each other in a triangle. Plant and other organic material was subsequently removed from core samples.

Physical–chemical parameters were measured using a multiparameter monitoring system (YSI model 85, Columbus). Metal concentration analysis had been previously performed by our laboratory (da Silva Pereira, 2016) and found no significant difference in metal concentrations relative to background within the Serinhaém estuary. For each sample an aliquot was taken for DNA extraction, while the remainder was used for measuring organic matter content using “loss on ignition” (Nelson and Sommers, 1996). The total genomic DNA was extracted from 0.25 g of sediment using the PowerSoil DNA Isolation Kit (Qiagen, Carlsbad, CA, USA).

2.3 Library preparation and sequencing

After DNA extraction, we amplified the V4 region of the prokaryotic 16S rRNA gene using PCR with the primer pair 515F-Y (Parada et al., 2016) and 806R-XT (Caporaso et al., 2011) using KAPA HiFi for 25 cycles (see Supplement methods, file 1). Illumina sequencing adapters and dual-index barcodes were added to the amplicon target using Nextera XT according to the manufacturer's directions (Illumina, San Diego, CA, USA). DNA sequencing was performed using the Illumina MiSeq platform, V2 kit (300 cycles).

2.4 Data analysis

2.4.1 Sequence trimming, denoising, and amplicon sequence variant (ASV) clustering

Trimmomatic (Bolger et al., 2014) was used to filter and trim demultiplexed sequences. Paired reads were combined using QIIME 2 for reads with an overlap of greater than 9, while paired reads with an overlap of greater than 6 were combined using a custom script (see Supplement methods, file 1). Reads were denoised using DADA2 (Callahan et al., 2016) in QIIME 2 (Bolyen et al., 2019) and then clustered into amplicon sequence variants (ASVs). Alpha rarefaction was calculated using QIIME 2 (Fig. S1 in the Supplement). We performed a variety of alpha-diversity (Figs. S2–S4 in the Supplement), beta-diversity (Fig. S5 in the Supplement), and hierarchical correlation clustering (Fig. S6 in the Supplement) tests using QIIME 2.

Figure 2 Normalized taxonomic abundances from each tidal zone. Taxa identified with an abundance higher than 1 % are shown as stacked bar plots; the horizontal axis is the phylum level while the stacked bars are the class level. Taxa abundance has been normalized by downsampling sequence depth.

2.4.2 Taxonomic assignment, community visualization, and environmental tests

Taxonomic assignment was performed using QIIME 2's naive Bayes scikit-learn classifier (Bokulich et al., 2018) trained using the 16S rRNA gene sequences in the SILVA database (Silva SSU 132) (McDonald et al., 2012) (see Supplement methods, file 1). Additional QIIME 2 visualizations for (1) ASV abundance (Fig. S7 in the Supplement), (2) proportional representation between zones (Supplement file 3), and (3) taxonomy (Supplement file 4) are available as a Supplement files. Phylogenetic reconstruction was carried out in QIIME 2 using the feature classifier trained on 16S rRNA gene sequences. All groups were required to be present within at least two samples with a minimum of three reads each. In one event a taxonomic family was manually relabeled from “Uncultured eukaryote” to “Uncultured Bathyarchaeia” (Figs. 2 and S8 in the Supplement).

The vegan R package (Dixon, 2003; Oksanen et al., 2015; R Core Team, 2019) was used to test correlations between community structure and environmental variables. Distances were calculated using metaMDS (distance used was Bray–Curtis), and environmental variables were fit using envfit (Table S2 in the Supplement). Additionally, we also tried a constrained ordination analysis using vegan's capscale method. However, the linear model fit of the complete data set using this approach did not identify a statistically significant correlation (Fig. S9 in the Supplement).

2.4.3 Functional analysis using PICRUSt2

Functional analysis was performed using PICRUSt2 (Douglas et al., 2019; Barbera et al., 2019; Czech et al., 2020; Louca and Doebeli, 2018; Ye and Doak, 2009) with default settings. As we used 16S rRNA gene amplicon sequencing a crucial limitation must be considered in evaluating our results (Sun et al., 2020), which is that all KEGG ortholog (KO) abundances (and subsequent pathway analysis) are derived from the inferred abundance of genes from the closest reference taxa that matches our supplied taxa generated by QIIME 2. As such, there is the potential for disagreement between the actual organism and the reference we rely on for inference. To address this we use PICRUSt2's nearest sequenced taxon index (NSTI) as a heuristic cutoff in the evaluation of taxon level functional analysis (Sect. 3.3, Table S1 in the Supplement). Although we are using taxonomic families in this study, we will continue to use the NSTI cutoff of 0.15, despite it being derived for species level comparisons. Where relevant families with median NSTI scores above 1 SD are labeled with an asterisk (*), families with an NSTI outside of this range were not included in this analysis. A zone-specific analysis of significant differential abundances was performed using the general linear model method of the ALDEx2 package (see Supplement methods, files 1 and 6). A heatmap of KOs with significant differential abundance between zones was then generated (Fig. 6).

2.4.4 Zone-specific taxonomic and functional enrichment

In order to identify which species were significantly different in abundance in each zone, we performed taxa enrichment analysis (Fig. 4), using a custom Python (Van Rossum and Drake, 2009) script (see Supplement methods, files 1 and 6). To determine which taxa were driving the differences in functional abundance, we calculated taxa-specific KO enrichment (Fig. 4) using the relative functional abundances from PICRUSt2 (see Supplement methods, files 1 and 6). KOs were limited to the KEGG metabolic pathways (ko00030, ko00195, ko00680, ko00710, ko00720, ko00910, ko00920). KOs were deemed enriched in a taxa if that taxa accounted for at least 10 % of the total relative functional abundance for that KO. A taxa was deemed to be enriched for a metabolic pathway (Fig. 4) if it was enriched for three KOs within that pathway. A more detailed analysis of metabolic pathway modules was performed using KEGG Mapper reconstruction to identify instances of taxa enrichment that included complete pathway modules (Fig. 5, Supplement file 7).

2.5 Accessibility

The entire computational workflow is available on GitHub: https://github.com/pspealman/COSantana_2021, last access: 25 March 2021.

The data have been deposited as PRJNA608697 in the NCBI BioProject database: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA608697, last access: 25 March 2021.

3 Results

3.1 Structural composition of prokaryotic communities

Sequence clustering yielded a total of 2030 ASVs, of which 94.3 % were assigned to Bacteria, 5.3 % were assigned to Archaea, and 0.4 % were not assigned to any prokaryotic kingdom (Fig. S7). All sites combined resulted in 36 phyla, 158 orders, 79 classes, 164 families, 135 genera, and 95 species. Approximately 89 % of all the sequences belonged to 6 phyla: Firmicutes (33.7 %), Proteobacteria (32 %), Chloroflexi (7 %), Planctomycetes (5.9 %), Actinobacteria (5.7 %), and Crenarchaeota (4.4 %). Crenarchaeota, Thaumarchaeota, and Asgardaeota were the representatives of Archaea groups with more than 1 % abundance in the samples (Fig. 2). Taxonomic profiles are broadly similar across the tidal zones at the family level (Fig. S10 in the Supplement). However, several families break from this trend (Supplement file 7), such as the two families of the order Alteromonadales (Fig. S10 – 2. Marinobacteraceae and Pseudoalteromonadaceae) are only abundant in the supralittoral zone. From the class Deltaproteobacteria (Fig. S10 – 9. Desulfarculales, 10. Desulfobacterales, 14. Myxococcales, 23. Syntrophobacterales), the majority of orders and families are either more abundant or only present at the sublittoral zone as well.

Figure 3 Correlation between environmental variables and prokaryotic communities. Dots represent taxonomic abundances per site as plotted by vegan metaMDS using nonmetric multidimensional scaling (NMDS). Arrow length is a representative of the predictor strength of environmental variable vectors, with red arrows having statistical significance as calculated by envfit (p value < = 0.05).

Alpha-diversity analysis showed a trend of the intertidal zone having the lowest prokaryotic diversity (Shannon's index – 3.8), while the sublittoral zone had the highest (Shannon's index – 4.8) and the supralittoral region had intermediary diversity (Shannon's index – 4.2) (Fig. S2). However, these differences between diversity indices were not statistically significant between the tidal zones, suggesting that the communities are similar in terms of biodiversity. Similarly, the beta-diversity analysis revealed no statistical differences in the taxonomic structure of prokaryotic populations between zones (Fig. S5).

3.2 The influence of environmental variables

We also aimed to determine if population structures across mangrove tidal zones correlated with the abiotic environmental variables of salinity, water content, organic matter, and temperature from each tidal zone (Table S2). The results revealed a significant correlation between the prokaryotic populations with salinity and organic matter (Fig. 3, Table S3 in the Supplement). Specifically, increased organic matter was positively correlated with sublittoral population abundances, while increased salinity was positively correlated with supralittoral population abundances. Neither water content nor temperature measures reflect a significant difference in community structure between zones.

Figure 4 Zone-specific measures of taxa showing significant enrichment in sediments of different tidal zones. Here we show prokaryotic families that have significantly different abundances between zones and whose mean effect size exceeded 10 % in at least one zone. To be labeled with a metabolic pathway, taxa were required to have at least 10 % of three KOs in that pathway in any zone. Families with a median NSTI within 1 SD of 0.15 are labeled with *.

Figure 5 Zone-specific enrichment of pathway modules. Here, we show the impact that the differential abundances of major families (Fig. 4) have on potential functional profiles in each zone at the level of pathway modules. To be included here each module must be completely enriched within a single taxon in a single zone, such that the taxa account for at least 10 % of the total potential functional abundance for each KO involved in the pathway. Families with a median NSTI within 1 SD of 0.15 are labeled with *.

3.3 Tidal-zone-specific taxa enrichment and their contributions to metabolic functional potentials

To determine which specific taxa were significantly different in abundance between zones, we performed a zone-specific enrichment test for each taxa (Fig. 4). At the taxonomic level of family, taxa that were a substantial component of the population of at least one zone (at least 10 % of the population) and whose abundance in another zone was significantly different were defined as zone-specific enriched taxa. Additionally, these taxa were also evaluated for their role in metabolic processes; those with enrichment of metabolism-associated KOs (accounting for more than 10 % of the total of a given KO for at least three KOs of a single metabolic pathway) were then labeled for that pathway. We also performed pathway analysis wherein a pathway module would be defined as enriched if the complete set of KOs for that pathway was enriched (Fig. 5).

We found eight families with statistically different population abundances between tidal zones (Bacillaceae, Flavobacteriaceae, Micrococcaceae, Clostridiaceae, Peptostreptococcaceae, Desulfobacteraceae, Rhodobacteraceae, and Vibrionaceae); four of these families (Flavobacteriaceae, Micrococcaceae, Clostridiaceae, and Vibrionaceae) had their lowest population abundances in the intertidal zone, and only one family (Bacillaceae) had its highest abundance in this zone. No archaea satisfied the statistical criteria necessary for inclusion. This dominance by Bacillaceae and reduced abundance of other families is consistent with our previous observations of the intertidal zone having a low diversity. Furthermore, these differences in population abundance can also affect the functional profile of these zones. Each family also contributed substantially (≥ 10 %) to the total potential functional abundance of at least one of the metabolic pathways in one zone (Fig. 4). Of those eight families with the exception of Peptostreptococcaceae, we found each contributed significantly to carbon metabolism, which includes photosynthesis (Ph, two families), carbon fixation in photosynthetic (Co, six families) or heterotrophic (Cp, seven families) organisms, and methane metabolism (M, seven families). All families contributed to sulfur metabolism-associated KOs, and these were mostly associated with assimilatory sulfate reduction (M00176) and dissimilatory sulfate metabolism (M00596). Seven families contributed substantially to nitrogen metabolism-associated KOs. Only six families contributed importantly to the general phosphorus metabolism-associated KOs. Taken together, this suggests that zone-specific taxa enrichment may also contribute to differential metabolic activities in these zones.

We can also analyze KO enrichment at a higher resolution by considering taxa enrichment of metabolic modules (Fig. 5). We find greater taxa-enriched KOs in the sublittoral zone, where 24 out of 29 (82.8 %) taxa-enriched pathways are found, vs. 21 (72.4 %) and 22 (75.9 %), in the intertidal and supralittoral zones, respectively. Notably, Bacillaceae is the most enriched family in each zone and, as such, has the potential to drive a large proportion of functional activity, being enriched in 17 metabolic modules in eight pathways. These include modules from the carbohydrate metabolism pathways (citrate cycle and pentose phosphate), as well as from sulfur and nitrogen metabolism pathways. Notably, some of these are enriched in only one zone, such as sulfur metabolism and the citrate cycle, suggesting that these roles are occupied by other families in these tidal zones.

Figure 6 Zone-specific measures of metabolism-associated KEGG orthologs. The heatmaps show the metabolism-associated KOs that have significant differences in functional abundance between zones for carbon, phosphorus, sulfur, and nitrogen pathways. Notably, the majority of KOs (45/54, 83 %) have their highest relative functional abundance in the sublittoral zone, with the intertidal (5/54, 9.3 %) and supralittoral (4/54 7.4 %) zones as near-equal minorities.

3.4 Zone-specific differences in metabolism-associated KEGG orthologs

To determine if there were significant differences in potential metabolic activity between tidal zones, we calculated the potential functional abundance of metabolic KOs for each zone. KOs with significantly different functional abundances between zones are shown in Fig. 6. Overall, the functional profiles showed that the majority of metabolic pathways for the transformations of carbon, nitrogen, phosphorus, and sulfur are more abundant in the sediments of the sublittoral zone. Taken together, the results suggest a positive correlation between biodiversity and potential metabolic functional abundances.

4 Discussion

In this work, we extended the study of preserved mangrove areas to characterize the variance in prokaryotic populations within microhabitats defined by the influence of the tidal regime. Previous research has shown that mangrove forests exhibit zonation driven by different biotic and abiotic factors which influence microbial community structure. However, most of this research has been conducted in anthropogenically impacted areas (Pupin and Nahas, 2014; Marcial Gomes et al., 2008; Alzubaidy et al., 2016; Rocha et al., 2016; Ceccon et al., 2019; El-Tarabily, 2002; Imchen et al., 2017; Lin et al., 2019; Zhang et al., 2018; Zhou et al., 2017), thus confounding the makeup of the microbial population structures. Conversely, studies that sought to identify differences induced by pollution and urbanization on mangroves found broad differences in prokaryotic populations in impacted areas compared to preserved mangrove areas (Pupin and Nahas, 2014; Nogueira et al., 2015) but did not study the population differences of distinct microhabitats within mangroves.

Our work shows that the different tidal zones of the Serinhaém estuary's mangrove have similar community structures with distinct functional profiles. We observed a prevalence of Firmicutes in the sediments of all tidal zones, accounting for 34 % of all reads. This differs from previous studies of mangroves, both around the world and in Brazil (Alzubaidy et al., 2016; Andreote et al., 2012; Mendes and Tsai, 2014; Nogueira et al., 2015; Ceccon et al., 2019; Imchen et al., 2017; Cheung et al., 2018; Zhou et al., 2017), which rarely found the relative frequency of Firmicutes to be above 10 % (Figs. S11–S13 in the Supplement). Instead, these studies found Proteobacteria to be the dominant phyla and, while Proteobacteria was also abundant in our samples (32 % of all reads), it never exceeded the abundance of Firmicutes. Notably, Firmicutes were mainly represented by the family Bacillaceae which in our potential functional analysis was an important driver of metabolic cycles. We also compared our observations to coastal marine sediments (Supplemental file 2), and here too, we found Proteobacteria to be the dominant taxa at anthropogenically impacted sites. However, in one study performed at a marine sanctuary there was also the observation of high-abundance Firmicutes (18 %–31 %).

Previous work on contaminated mangroves in Brazil has suggested that Gammaproteobacteria is a bioindicator of an anthropogenically impacted mangrove (Andreote et al., 2012a). Comparing normalized reads across our samples for all Firmicutes and Gammaproteobacteria (Table S4 in the Supplement), we found that the two taxa are significantly negatively correlated (Pearson's R = −0.84; p value = 0.004), consistent with Firmicutes being a bioindicator of a pristine environment. However, Firmicutes were also found at elevated levels at several anthropogenically impacted sites (Torres et al., 2019; Haldar and Nazareth, 2018; Tiralerdpanich et al., 2018). While it could be that their presence is due to the pristine nature of the area, it may also be due to the local climate (Nogueira et al., 2015), an unmeasured environmental factor (Tong et al., 2019), or the presence of another organism such as a non-prokaryote or prokaryote which is difficult to detect using 16S rRNA gene amplicon sequencing (Eloe-Fadrosh et al., 2016; Zhang et al., 2021).

The results show that both salinity and organic matter were significantly correlated with communities in different tidal zones, confirming the hypothesis that changes in physicochemical parameters lead to changes in the prokaryotic communities inhabiting these sediments. It has been previously observed that tidal variations can cause nutrient washout or erosion (Behera et al., 2019; Zhang et al., 2020), which in turn can change the bioavailability patterns for organic matter and essential nutrients along the distinct tidal zones of a mangrove habitat. Salinity levels are also influenced by the tidal variations, and it has been previously determined that the increased salinity as well as the prevalence of humic-like substances can interfere with the accessibility of the microbial community to the sediment organic matter (Wang et al., 2010; Ceccon et al., 2019).

It is important to note that other environmental variables such as granulometry, vegetation, and pollutant distributions may also impact mangrove sediment communities, as observed previously (Peixoto et al., 2011; Colares and Melo, 2013; Rocha et al., 2016). Furthermore, it is also possible that they are influenced by additional biological factors, such as fungi and other eukaryotic microbes (Simões et al., 2015), and plant rhizome contamination (Bennett and Klironomos, 2019; Miller et al., 2019), as observed in previous work (Rocha et al., 2016; Zhang et al., 2017). Thus, our understanding of prokaryotic community structures will be greatly increased if complemented with rhizome and eukaryotic population information.

It is also important to notice that the use of 16S rRNA gene amplicons potentially impairs the identification of some important prokaryotic groups, especially from the domain Archaea, which in some environments can reach 10 % of the total microbial community (Eloe-Fadrosh et al., 2016). Recent research has emphasized how the use of 16S rRNA gene amplicon sequencing can lead to failure in detecting some archaeal populations which can play important roles in nutrient cycling (Zhang et al., 2021). Taken together, however, this work suggests that a variety of environmental factors, potentially driven by tidal cycles, can generate niche variations that influence the structure and function of communities.

Although the differences in diversity indices were not statistically significant, the data show a clear trend of differentiation between zones that suggests that further investigation, with a larger number of samples, would be sufficient to resolve differences with statistical significance. The highest biodiversity was found in the sublittoral mangrove sediments, while the intertidal zone had the lowest biodiversity. The finding of taxonomic groups that were more abundant and/or exclusive to specific zones is also suggestive of the importance of the physical–chemical parameters in shaping the microbial communities. The prevalence of the groups of Deltaproteobacteria in the sublittoral zone, for example, could be associated with the prevalent anaerobic conditions of such sediments, which could result in the selection of specific prokaryotic groups such as sulfate-reducing bacteria (Andreote et al., 2012), as seen with the increased abundance of Desulfobacteraceae, Desulfobulbaceae, and Syntrophaceae.

Our data suggest that the intertidal sediments of mangrove forests have lower prokaryotic diversity than those in the constant environments, such as the supralittoral and sublittoral regions. Models of biodiversity in dynamic environmental regions, such as intertidal zones, have suggested that oscillating environments may support higher biodiversity than that observed in static environments (Brose and Hillebrand, 2016). Furthermore, studies in anthropogenically impacted mangroves have often found substantial diversity in the intertidal zones as well (Lin et al., 2019; Zhang et al., 2018). An alternative model is that the constant shifts in environmental parameters can create harsh environmental conditions that select for resistant taxa thus reducing the relative biodiversity (Statzner and Moss, 2004), which is consistent with our results. Considering other surrounding environments such as the shallow marine sediments (Supplemental file 2), we may expect the diversity to decrease in comparison to the mangrove sediments, potentially with an increase in Proteobacteria, a trend that has been reported in previous studies which found similar groups prevailing in mangrove sediments (Wang et al., 2012; Behera et al., 2019).

Community functional profiling was undertaken using 16S rRNA gene sequences and identified higher metabolic potentials for the cycles of C, N, P, and S in the sublittoral region. The higher abundances of KOs found in the sublittoral zone is possibly due to the greater taxonomic diversity that was also observed for this region. Notably, this result is in agreement with the understanding that many microbiome functions in ecosystems are redundant (Tolkkinen et al., 2020; Barnes et al., 2020). It is important to note that these functional profiles are not the result of metagenomic sequencing; as such it is not necessarily true that all members of the taxa have the same genome as the reference holotype, nor is it necessarily true that these genes would be expressed equally, even if present, and this should best be understood as a measure of potential metabolic activity. Studies that have performed functional profiling using metagenome sequencing on mangrove sediments resulted in comparatively greater resolution (Zhang et al., 2021), although the use of 16S rRNA gene sequencing has been proved to provide comparable results for the recovery of general patterns of pathway abundances (Jovel et al., 2016).

Generally, diverse communities with organisms possessing redundant metabolic functions may be more stable against environmental perturbations, as the organisms will respond differently to stressors, increasing the likelihood of survival of some taxa (Girvan et al., 2005). Thus, the identification of the intertidal zone as having the lowest diversity suggests it may be in a more precarious position. This concern is amplified as the water itself is frequently the carrier of contamination, both from rivers, as is the case for urban waste (Yunus et al., 2011), and from the oceans through the tides, as the case with oil spill contamination (Cabral et al., 2016). Thus, it is important to consider that different parts of the mangrove tidal zone could be exposed to different levels of contamination and that this could, in turn, affect the organisms in a zone-specific manner.

The identification of enriched taxa with similar functional potentials for the diverse nutrient metabolisms is in agreement with previous observations of the microbial diversity in mangrove sediments (Rocha et al., 2016; Ceccon et al., 2019; Cabral et al., 2016; Mendes and Tsai, 2014; Zhao and Bajic, 2015). Conceptually, one could interpret that this higher-order redundant biodiversity could grant the microbiome a higher level of resilience when challenged by environmental perturbation, as organisms with similar metabolic pathways have the potential to maintain the ecological functions even if some taxa are lost.

In conclusion, the results suggest that communities with distinct diversity and functional profiles occupy the sediments of the different mangrove tidal zones. This variance is significantly associated with the abiotic environmental variables of salinity and organic matter. Our functional potential analysis suggests there could be significant variation in metabolic potentials between tidal zones. We observed that the intertidal region, the most biophysically dynamic zone in the mangrove forest, presents the least biodiversity. One possible explanation for this is that the dynamic, cyclic, tidal environment is itself harsher than that found in the other tidal zones. This study of pristine mangrove tidal zones suggests that some microhabitats may be more sensitive to anthropogenic impact than others. Preservation, conservation, and remediation measures of mangroves should take into account the fragility of the intertidal zone and the increased threat this zone faces from soluble pollutants, such as urban sewage, fertilizer, and pesticide runoff, as well as from buoyant pollution such as plastics and oil spills. However, it remains an open question if other microhabitats, such as the sublittoral zone, have greater resilience to environmental perturbations given their greater abundance of taxonomic diversity and metabolic function. Conceptually, this property of a dynamic environment defining a selective niche, similarly to a physical barrier, is worthy of further study. Further exploration of whether different groups of organisms found to have metabolic potentials are actively participating in the element cycles in these mangrove sediments will also benefit the field.