Spatiotemporal dynamics of benthic bacterial communities in the Perdido Fold Belt, Northwestern Gulf of Mexico

This study analyzed the spatiotemporal dynamics of surficial benthic microbial communities in a bathymetric gradient (44 - 3573 m) across four oceanographic campaigns at the Perdido Fold Belt (PFB) in the northwestern Gulf of Mexico (nwGoM). Bioinformatic analysis of 16S rRNA gene amplicons grouped the 27 samples into three clusters according to a longitudinal bathymetric gradient. Differences in community structure among clusters, based on PERMANOVA analysis, were partially explained by cruise, water depth, temperature, salinity, nitrate plus nitrite, silicate, redox potential, Ni, Cd, Pb, and Al, as well by aliphatic and aromatic hydrocarbon concentrations. Into microbial community composition, Gemmatimonadaceae, Planctomycetaceae, and the JTB255 were detected at all depths across the four campaigns. Members of Anaerolinaceae and specific sulfate-reducing bacteria were more abundant in sites located between 43 and 1200 m, and Rhodospirillaceae, wb1-A12, OM1 clade, Desulfurellaceae, Gemmatimonadetes, Nitrospinaceae, and Clostridiaceae 1 were better represented in deeper sites. Alpha diversity was similar between the three groups and remained stable; however, 10 samples presented changes in the community structure across the four campaigns. Finally, a multivariable association analysis revealed 25 bacterial genera positively related with physicochemical parameters that characterized the environment from shallow to deep sea sites. Taken together, these results yield insights into the temporal stability of 17 of 27 sites in the PFB and revealed signature taxa with putatively ecological relevance in sedimentary environments.


Introduction
Marine sediments are highly variable environments in space through time. Their physical and chemical settings are affected by human activities and natural processes such as fluvial inflows and oceanographic phenomena (e.g., vertical eddies, currents, horizontal vortices) so climate and seasons are other important factors for their variation (Hernańdez-Molina et al., 2016;Ward, 2017). The ocean food web depends on the circulation of nutrients and organic matter between surface waters and sediments (Zinger et al., 2011;Orsi, 2018). In marine sediments prokaryotes are the major biological catalysts that largely control organic carbon and nutrients cycling and are involved in the climate change regulation and the degradation of multiple pollutants (Newton et al., 2013;Orsi, 2018;Cavicchioli et al., 2019;Li et al., 2020;Walker et al., 2021). Microbial studies have evidenced temporal and spatial patterns in community structure in response to environmental and biotic factors (Ortmann and Ortell, 2014;Torres-Beltrań et al., 2016;Zorz et al., 2019). Given their essential roles in regulating ecological services, understanding the factors that shape microbial community structure on spatial and temporal scales is crucial to yield insights to anticipate the responses of marine ecosystems to environmental changes afforded (Zinger et al., 2011). Under this idea, deciphering the influence of environmental variables on the community structure of native prokaryotes became an area of great research interest in the Mexican GoM, such as the PFB in the northwestern Gulf of Mexico (nwGoM), which is within the third most important hydrocarbon fields in the basin (Celis-Hernandez et al., 2018).
To date, few studies have assessed microbial communities in the PFB that combined molecular analysis with physicochemical data at a single time point, and none have compared the structure of benthic microbial communities over time. Although the way to analyze, interpret, and present microbial diversity data in these studies are difficult to compare, some interesting trends can be observed. In the first place, the analysis of the beta-diversity of pelagic and benthic communities showed a tendency to cluster the samples by water column depth (Sańchez-Soto Jimeńez et al., 2018;Raggi et al., 2020;Ramıŕez et al., 2020;Sańchez-Soto et al., 2021). Sańchez-Soto Jimeńez et al. (2018; and Ramıŕez et al. (2020) reported the co-occurrence of aerobic and anaerobic metabolisms in the studied sites; however, anaerobes (e.g., sulfate-reducers) showed higher abundances in shallow environments where anoxic conditions predominated. Changes in the community structure in such studies were mostly related to differences in water depth and temperature, that were common factors in determining the betadiversity of benthic communities. Taxonomic and metabolic analysis have revealed indigenous microorganisms to participate in biogeochemical cycling and the presence of putative hydrocarbon degrading bacteria was also a common result in all these studies, suggesting that the PFB could be an area affected by the presence of hydrocarbons in the environment (Ramıŕez et al., 2020;Sańchez-Soto Jimeńez et al., 2018;Raggi et al., 2020;Loza et al., 2022). The study of Sańchez-Soto et al. (2021), based on sequence analysis of the dsrB gene encoding for the beta subunit of the dissimilatory sulfite reductase involved in sulfate-reduction, showed that differences in the community structure of sulfate-reducing microorganisms (SRMs) were related to the sediment concentrations of Ni and Cd.
Diversity measurement is fundamental for understanding community structure and dynamics, but it is particularly challenging for microbes (Lozupone and Knight, 2008). Alpha diversity is often estimated as the number of species in a community (species richness), and beta diversity is usually based on the number of shared species (Lozupone and Knight, 2008). Despite recent progress in assessing community composition and spatial distribution of benthic prokaryotes in the nwGoM, the study of alpha and beta diversity and the factors that shape microbial assemblages continues to be an area of acute research. Analysis of benthic prokaryotes in a bathymetric gradient in the PFB has shown the response of microbial communities to specific parameters variation, including water depth, temperature, OM, salinity, sediment particles, nutrients, and contaminants (Sańchez-Soto Jimeńez et al., 2018;Raggi et al., 2020;Ramıŕez et al., 2020;Sańchez-Soto et al., 2021;Loza et al., 2022). However, how consistently microbial communities respond to such variables over different climatic seasons and in a wider spatial area is still to be investigated. Thus, the objectives of the present study were to 1) analyze beta-diversity, 2) determine the composition and distribution of benthic prokaryotes, 3) characterize the environmental variables of the region identifying their statistical association with changes in the structure of microbial communities across time and space.

Sampling
Four oceanographic campaigns were conducted in the Mexican region of the Perdido Fold Belt (PFB) polygon located in the nwGoM ( Figure 1): Perdido 1 (P-01, May 2016), Perdido 2 (P-02, September-October 2016), Perdido 3 (P-03, June 2017), and Perdido 4 (P-04, September 2017). During each campaign, 27 soft bottom sediment samples were collected using a Hessler-Sandia box corer from four perpendicular transects (B, C, D and F) to the Tamaulipas coastline extended in a bathymetric gradient (~44-3573 m) (Figure 1). Salinity (PSU), temperature (°C), and dissolved oxygen (DO mL L -1 ) were determined in situ from the bottom water using a CTD Seabird 9 plus ® . Bottom water samples were collected using Niskin bottles and kept frozen until their analysis in the laboratory for nutrient determinations such as nitrate plus nitrite (NO − 3 + NO − 2 mM L −1 ), orthophosphate (PO 3− 4 mM L −1 ), and silicate minerals (SiO 4− 4 mM L −1 ) which were determined from the bottom water by colorimetric methods following standard protocols using a spectrophotometer UV-VIS (Agilent Technologies Cary 60) (Herguera-Garcıá et al., 2021). Following the box corer recovery, pH and redox potential were measured directly in the sediment with specific sensors (Extech pH 100 probe and an Extech RE300 probe, respectively) (MA, USA). The first 10 cm of each sediment sample were removed using sterile syringes with cutoff tips and preserved at −20°C for subsequent DNA extraction. A sediment subsample of 400 g was taken and preserved at -4°C in plastic bags for later organic matter (OM%) and grain size determinations (%). Another two subsamples of 100 g of sediment were taken and stored at -4°C for heavy metals and hydrocarbons content determinations. The subsamples for heavy metal analysis were stored in plastic bags previously washed with a solution of HNO 3 1N (pure grade from Sigma-Aldrich) and deionized water, and the subsamples to analyze the hydrocarbon content were stored in glass bottles previously washed with hexane and acetone separately (both of chromatographic grade from Sigma-Aldrich).
Spatial and temporal differences (ANOVA, P<0.05) of environmental variables were tested using Zone and Campaign as grouping factors, respectively (Supplementary Tables 1, 2). Campaign corresponded to cruises P-01, P-02, P-03, and P-04, and Zone corresponded to shallow (≤200 m), intermediate (201-1200 m) and deep (1201-3573 m) categories that were defined along a longitudinal gradient according to the total water depth of each sampling site (Supplementary Material 1). Environmental variability was further explored with a principal component analysis using all the environmental parameters ( Figure 2). Analyses of environmental variables were performed with the InfoStat program (version 2020). Figures were generated using InfoStat and Ocean Data View (version 2022) available in https://odv.awi.de.

DNA extraction and 16S rRNA amplicon sequencing
Sediment subsamples were thawed and centrifuged for 1 min at 10,000×g to extract the genomic DNA from the settled cells after discharging the remaining water. Genomic DNA was extracted from 1 g of each sediment sample with the DNeasyPowerSoil Kit (QIAGEN, Gilden, Germany) following the manufacturer's instructions and stored at -20°C until needed. The integrity of the extracted DNA was verified by 1% agarose gels.
From the extracted DNA, 16S rRNA amplicons were prepared with the adapters for sequencing by Illumina MiSeq system following the two-step PCR protocol. In the first PCR the primer pair S-D-Bact-0341-b-S-17 and S-D-Bact-0785-a-A-21 (Klindworth et al., 2013) was used to generate approximately 550 bp paired-end reads, covering the V3-V4 16S rRNA region, in a Veriti thermal cycler (Applied Biosystems Veriti ABI Inc., Foster City, CA, USA). The PCR program was as follows: initial denaturation at 95°C for 3 min, 25 cycles including denaturation for 30 s at 95°C, alignment for 30 s at 55°C, and elongation for 30 s at 72°C, with a final extension at 72°C for 5 min. Each PCR reaction mixture (20 µL) included 1 µL of each primer (0.5 mM), 10 µL of 2x phusion High-Fidelity MasterMix (Thermo Scientific, USA) and 2 µL of DNA (~10 ng/mL). PCR products were then purified with AMPure XP beads technology. For the second PCR round, eight cycles were included to attach the Nextera indexes with the XT Index Kit reagents. The size of the amplicons was verified on a QIAxcel Advanced system (QIAGEN, Hilden, Germany) and the Sampling sites in the Gulf of Mexico. Region of study, located from 25.74°to 24°North latitude and from -97.53°to -95.18°West longitude, with the 27 sampling sites distributed in four perpendicular transects (B, C, D, and F) to the Tamaulipas coast that were settled along a bathymetric gradient (43-3573 m). Black crosses indicate the intersection between coordinates. Gray lines indicate isobaths. concentration was quantified by fluorometry on a Qubit 3.0 kit (Life Technology, Shah Alam Selangor, Malaysia). Indexed amplicons were purified with AMPure XP beads and individually diluted in 10 mM Tris (pH 8.5) and combined to an equimolar concentration of 9 pM. The final libraries were loaded on a V3 MiSeq Reagent Kit V3 flow cell (600 cycles) and sequenced on 2x300 Illumina-Miseq platform (Illumina, San Diego, CA, USA) at CINVESTAV Merida.

Bioinformatic analyses
Sequences were processed with the QIIME2 pipeline (Bolyen et al., 2019). The amplicon sequence variants were resolved with the DADA2 plugin (Callahan et al., 2016). Due to low quality detected in the last 150bp of some R1 sequences, merging with R2 was not possible and only R1 reads were used for this study. Trimming and truncating of R1 reads were made on positions 30 (end 3') and 150 (end 5'), removing chimeric sequences with the "consensus" parameter. The representative sequences of ASVs were classified with the "classify-consensus-vsearch" plugin (Rognes et al., 2016), using the SILVA database (v.128) as a taxonomic reference. A phylogeny of the representative sequences of ASVs was calculated with the "align-to-tree-mafft-fasttree" plugin (Katoh and Toh, 2008;Price et al., 2010). The feature table and phylogeny were exported to the R environment.
The data were processed in R with the phyloseq package (McMurdie and Holmes, 2013). The alpha diversity was evaluated using the number of ASVs and the Shannon diversity index. The bar plots of taxonomic abundance were performed and prettified with the ggplot2 packages (Wickham, 2010). A PCoA was calculated based on the weighted UniFrac distance. The gap statistic to determine the optimal number of sample clusters in the PCoA was performed with the cluster package (Maechler et al., 2014) using the partition around medoids (pam1) method and 1000 bootstraps. To determine variables that influence changes on microbial community structure among sites a PERMANOVA (Permutational Multivariate Analysis of Variance) with the adonis function analysis was used. LEfSe analysis was used to identify the features most likely to explain differences between clusters of samples based on statistical different abundances (Segata et al., 2011). The Maaslin2 package (Mallick et al., 2021) was used to associate the genera with the environmental variables using the parameters as default (method:LM, transformation: log2, normalization:TSS). The threshold P-value to consider a result significant was set up to<0.05.
Temporal changes (ANOVA, P>0.05) were also observed in the environmental parameters (Supplementary Table 2; Figure 2). Water temperature significantly increased during P-01 and P-02 in the intermediate and deep zones, respectively, while it decreased in the shallow zone in P-02 ( Figure 2). In P-02, salinity decreased only for the shallow zone, while the redox state decreased in all the region (mean ± standard deviation, shallow zone = -511.2 ± 77.2 mV, intermediate zone = -443.4 ± 104.8 mV and deep zone = -304 ± 327.1 mV). Differences in the OM% were detected only in the zones intermediate and deep, with the highest content in cruises P-01 and P02 compared to cruises P-03 and P-04 (Supplementary Table 1 decreased in P-01 and P-02 in the deep sites ( Figure 2). An increased content of Pb, Al, Ni, and Cd was recorded in P-02, while they decreased in P-01 ( Figure 2). V content showed temporal differences with the highest content in P-01, and the lowest in P-03. The highest concentration of PAHs was recorded in P-02, followed by P-01. The lowest concentration of PAHs was detected in P-03 and P-04. Aliphatic hydrocarbons showed statistical differences among cruises in the shallow and deep zones, the lowest concentrations were quantified in P-02 (Supplementary Table 2; Supplementary Figure 1). Medium (MS%), fine (FS%) and very fine sands (VFS%) presented a similar spatial distribution during the four cruises (Supplementary Table 1). However, they changed in time (Supplementary Table 2). The highest percentage of MS was obtained in P-01 and P-04 (mean ± standard deviation, 53.4 ± 14.4%), and the lowest proportion was obtained in P-02 (9 ± 14%). Meanwhile, P-02 presented the highest percentage of FS (68.4 ± 10.7%) and P-01 the lowest proportion (34 ± 13.1%) (Figure 2). The highest VFS percentage (20.8 ± 8.6%) was observed in cruise P-02, while the lowest proportion (4.9 ± 2.3%) corresponded to P-04 ( Figure 2). pH values ranged from 5.01 to 9.25 showing the highest values in P-01 and the lower ones in P-02 (Supplementary Figure 1). Overall, MS, FS, Cd, redox, and Al contributed most to the environmental variability between cruises according to the ANOVA results (Supplementary Table 2).

Microbial diversity and PERMANOVA
Community structure analysis grouped the samples into three clusters ( Figure 3A) differing in terms of the water depth (R 2 = 0.56, F 2 , 102 = 64.48, P<0.0001) at which the analyzed sediments were collected. In P-01 the samples were clustered following a clear bathymetric gradient; cluster I included the samples collected from 44 to 107 m of the water column, cluster II the sediments collected from 372.5 to 826 m of water, and cluster III the samples collected from 1000 to 3548 m of water (Supplementary Figure 2). Based on the community structure analysis, 17 of 27 samples showed temporal stability, since they were in the same cluster campaign after campaign (Figure 4; Supplementary Figure 2). The remaining 10 samples showed community structure shifts, and thus, they were grouped in different clusters in some campaigns. Most of these 10 samples were in transect D, which showed larger differences than the others, especially in P-02 and P-03 compared with P-01 (Supplementary Figure 2). Cruises, water depth, temperature, salinity, NO − 3 + NO − 2 , SiO 4− 4 , redox potential, the concentration of Ni, Cd, Pb, and Al, and the content of aliphatic and aromatic hydrocarbons were related to community structure changes. Nonetheless these parameters explained only a small proportion of these changes ( Figure 3B). Shannon diversity index (H') varied in the range from 5.4 to 7.5 (mean ± standard deviation, 6.7 ± 0.4) and was similar (P>0.05) between clusters across the four campaigns according to the ANOVA test (Supplementary Table 3). The species richness (S) based on the number of ASVs varied in the range from 402 to 2581 (1247 ± 391). P-01, P-03, and P-04 showed similar species richness between the three clusters. However, cluster II showed a significant increase (P<0.05) in the number of species during P-02 (Supplementary Table 3). Both alpha diversity metrics showed differences only in P-02 between transects with the highest values observed in transect D (H' = 1344 ± 128, S = 6.9 ± 0.06) (data not shown).

Differential abundances in microbial communities
LEfSe analysis identified ASVs with consistent differences between clusters. The ASVs with differential abundances between clusters were related to 57 microbial families (Figure 4). The most differentially abundant bacterial members in clusters I and II belonged to Nitrospiraceae, Anaerolinaceae, and potential sulfate-reducing bacteria affiliated to Syntrophobacteraceae, Desulfobacteraceae, Desulfarculaceae, and Desulfobulbaceae. However, Nitrospiraceae and Desulfobulbaceae were specially enriched in cluster I and Syntrophobacteraceae was the most abundant phylogenetic unit in cluster II with percentages higher than 10% in some samples ( Figure 4). Although clusters I and II had some families in common, LEfSe analysis indicated differences between their population structure. wb1-A12 and OM1 clade were detected as differentially abundant groups in clusters II and III. However, they were consistently more abundant in cluster III. Also of note was the relevant representation of sulfate-reducing bacteria from the Desulfurellaceae family in cluster III, as well as the abundance of Rhodospirillaceae with percentages higher than 10% in some samples. A similar distribution was found for Gemmatimonadaceae and Nitrospinaceae in cluster III samples. The group JTB255 MBG was relatively abundant at all depths, and therefore it was determined by LEfSe to be characteristic of all sites (Figure 4).

Bacterial genera and environmental variables
Analysis of multivariable association of microbial features revealed 25 genera correlated with specific environmental parameters ( Figure 5). The correlation of environmental variables with specific taxa was done at genus level because this kind of analysis must be performed at the lower taxonomic level to secure correspondence between metabolisms related with the environmental parameters. ASVs affiliated to Desulfatiglans and Sva008, as well as ASVs affiliated to uncultured bacteria from Syntrophobacteraceae, Latescibacteria, Anaerolineaceae, Nitrospiraceae, and an unassigned Dehalococcoidia genus were positively related with temperature, salinity, and the content of Pb, Al, and PAHs ( Figure 5). ASVs affiliated to wb1-A12, Rhodospirillaceae, Chloroflexi classes SAR202 and S085, NB1-j and SAR324 clade marine group B from Deltaproteobacteria, Acidobacteria class Subgroup22, Gemmatimonadaceae, Actinobacteria clade OM1, and the Desulfurellaceae genus H16 were positively related with water depth, NO − 3 + NO − 2 , SiO 4− 4 , redox, Cd, and Ni ( Figure 5).

Discussion
Two years of sampling at two different seasons allowed to evaluate the spatiotemporal dynamics of surficial benthic microbial communities along a bathymetric gradient (44 - Analysis of multivariable associations of microbial genera with environmental variables related with the community structure changes in PFB sediments. These genera were present in 25% of the total samples. Differential abundances analysis. Labels in y axis indicate the sample ID followed by the sampling campaign. Clusters I, II, and III, based on beta diversity analysis, are indicated on top. 3573 m). Alpha diversity analysis based on Shannon's diversity index showed spatiotemporal stability although the number of ASVs suggested that an increase in the species richness may reflect a larger environment variation. Beta diversity analysis revealed three clusters of samples and indicated that community structure changes were associated with sampling time, water depth, temperature, salinity, redox, and the concentration of nutrients (silicates and nitrate+nitrite), metals (Pb, Al, Ni, and Cd), and PAHs. Based on beta diversity analysis, 17 out of 27 samples showed temporal stability in the community structure and composition during the four sampling campaigns. However, 10 samples showed temporal variation in the community structure and composition suggesting underlying processes that influence microbial assemblages in the studied region. Moreover, microbial community analyses allowed the detection of genera with putative ecological relevance that correlated with specific environmental parameters. These aspects are discussed below.

Alpha diversity stability
Previous studies in the PFB estimating alpha diversity, have reported differences in microbial communities between shallow and deep sites based on relatively small data sets at a single time point (Sańchez-Soto Jimeńez et al., 2018;Ramıŕez et al., 2020). In the present study despite the environmental changes observed Shannon diversity index was similar between the bathymetric zones and remained relatively stable in the four sampling moments (Supplementary Table 3). The lack of spatial variability in the Shannon diversity between the bathymetric zones contrasts with a previous study in which higher microbial diversity was observed in shelf sediments (Sańchez-Soto et al., 2018). A strong physical mixing and seasonality were suggested to explain higher diversity in the shelf sediments (Sańchez-Soto et al., 2018). In the present study environmental differences were found between the bathymetric zones and the four campaigns, even between P-01 and P-03 conducted at the end of the dry season and between P-02 and P-04 conducted at the end of the rainy season (Figure 2;  Supplementary Tables 1, 2). This result suggests that microbial communities might develop independently of environmental changes as reported in other studies (King et al., 2013). It is possible that a stable diversity indicates a stable state of the ecosystem, since diversity is considered crucial in maintaining ecosystem functions (e.g., organic matter decomposition, nutrient cycling, and resistance and resilience to perturbations) (Loreau et al., 2001;Lozupone and Knight, 2008;Deng, 2012).
Temporal and spatial differences observed in the number of ASVs (species richness) indicated that species richness increased in cluster II during P-02 (Supplementary Table 3), and in transect D in this campaign (data not shown). The correspondence of a greater number of species with significant changes in environmental conditions observed in P-02 (Figure 2), suggests that changes in the species richness may reflect larger environmental variability (Walsh et al., 2016). Perhaps the increase in the number of species allows maintaining the stability of the system since, there seems to be a positive relationship between species richness and ecosystem temporal stability (Xu et al., 2021). The temporal environmental variation that influences microbial richness could be related to the latitudinal location of the transects, which has been found to explain sedimentary differences resulting from the prevailing hydrographic conditions throughout the year (Salcedo et al., 2017). The set of samples analyzed at the end of the dry season (P-01 and P-03) and at the end of the rainy season (P-02 and P-04) for two years, along with the relatively larger ocean area into the nwGoM, provided higher spatial and temporal resolution than previous works (e.g., Ramıŕez et al., 2020). However, further long-term studies are prompt to explore microbial benthic communities using water depth as an explanatory variable. This may allow to recognize whether the alpha-diversity metrics such as species richness can be used as sensitive indicators of environmental changes from shelf to the deep-ocean over time in the PFB.

Beta diversity in the bathymetric gradient
As documented previously, beta-diversity analysis at local and regional scales in the GoM has shown that benthic microbial communities living under comparable environmental conditions (e.g., water depth, temperature, nutrients, electron acceptors, etc.) may have greater similarity among them (Devereux et al., 2015;Sańchez-Soto Jimeńez et al., 2018;Overholt et al., 2019). In the present study, PCoA clustered the samples in three groups ( Figure 3A). The environment variation between the shelf, slope, and deep sea (Figure 2) may result in the community structure differences between these zones as indicated by the beta diversity analysis ( Figure 3A). As water depth, temperature, salinity, redox, and the concentration of nutrients (silicates and nitrate+nitrite), metals (Pb, Al, Ni, and Cd), and PAHs changed along the bathymetric zones, the abundance of prokaryotic species better adapted to different conditions may increase, while the abundance of less adapted species may decrease ( Figure 5). Weighted UniFrac distances suggested that benthic prokaryotes that changed in abundance between shelf, slope, and deep zones were distantly related ( Figure 3A). The low relatedness of benthic prokaryotes that changed in abundance between these zones may reflect alternative community assemblages conformed by species with different metabolic and functional requirements, and thus, with different functions in the ecosystem (Ortmann and Ortell, 2014).
Shelf and slope sites (43-1200 m) presented lower redox conditions and higher Al, Pb, and PAH concentrations, especially in P-02 (Figure 2). Higher metal and hydrocarbon concentrations in continental margin sediments are generally attributed to riverine inputs, that also affect the deposition of fine grain sediments (Trefry and Presley, 1976;Celis-Hernandez et al., 2018), which limit pore water exchange, and thus, favoring a lower redox state in these sites (Huettel et al., 2003;Burdige, 2007). Desulfarculaceae, Desulfobacteraceae, Syntrophobacteraceae, Nitrospiraceae, Anaerolinaceae, and Dehalococcoidia were significantly more abundant in these sites (Figure 5), and are recognized as anaerobes (Matturro et al., 2017;Hoshino et al., 2020;Walker et al., 2021). Desulfarculaceae, Desulfobacteraceae, and Syntrophobacteraceae are involved in the sulfur cycle through sulfate reduction, the dominant process in the anaerobic decomposition of OM in margin sediments (Parkes and Sass, 2007;Jørgensen et al., 2019). They may play a relevant role in PAHs degradation and metal precipitation (Geets et al., 2005;Shin et al., 2019) in the shelf and slope sediments. Nitrospiraceae members are nitrite-oxidizing bacteria and recently have been found to perform complete nitrification, oxidizing ammonia to nitrate (Daims et al., 2015). Anaerolinaceae has been related to higher OM loadings (Walker et al., 2021), although OM% was similar along the bathymetric zones in the present study (Supplementary Table 1), the lower redox state in shelf and slope (Figure 2) would indicate a higher activity of the oxidation processes (Alewell et al., 2008). Metabolic predictions from Dehalococcoidia cells have indicated their capacity to oxidize aromatic compounds, to be resistant to metals, and to use multiple electron acceptors that enable these bacteria to develop in low redox conditions (Wasmund et al., 2016). These taxa comprised the assemblage in most sites in the shelf and slope during the period of observation (Supplementary Figure 6; Figure 4). Their putative metabolisms seemed to correspond to the environmental conditions.
Deep sediments (>1200 m) presented higher redox conditions and nitrate+nitrite, silicate, Ni, and Cd concentrations (Figure 2). Higher redox values suggested oxic environmental conditions as reported previously for these sites (Sańchez-Soto Jimeńez et al., 2018). Under oxic conditions, nitrite usually accumulates because microorganisms may perform nitrification coupled to ammonia oxidation (Orsi, 2018). Silicate is associated with biological production (Sarmiento et al., 2004), and potential toxic Ni and Cd are a growing concern in the GoM owing to anthropogenic discharges (Ruiz-Fernańdez et al., 2019;Arcega-Cabrera et al., 2022). The abundant (resident) taxa in this zone were represented by Rhodospirillaceae, that has been related to an increase of inorganic nutrients (Shi et al., 2017), Gemmatimonadaceae, with aerobic bacteria capable to remove phosphorus (Zhang et al., 2003), Clostridiaceae that includes bacteria related to respiration of sulfur species (Muyzer and Stams, 2008), and OM1 clade, related to nitrogen fixation to produce ammonia (Dworking et al., 2006). In addition, a high abundance of wb1-A12 was detected which is capable of nitrite oxidation under oxic and anoxic environmental conditions (Orsi, 2018), and Desulfurella, related to sulfate-reducers that has been detected in high abundances in the deep biosphere (Orsi et al., 2016). Except for Desulfurellaceae, these taxa have been detected in high abundances in the study region (Sańchez-Soto et al., 2018;Ramıŕez et al., 2020;Sańchez-Soto et al., 2021). For Desulfurella this is the first report in the GoM, but sulfate-reducers seem to have an important role in oxic and metal rich deep sediments in the PFB (Sańchez-Soto et al., 2021). These abundant microbes found in most of the deep sites throughout the observation period, may perform relevant microbial functions under the prevailing environmental conditions, as posed for microbial communities in other marine ecosystems (Gobet et al., 2012). Thus, the detection of these abundant taxa with putative ecological relevance at most sites in the deep zone in each campaign suggested a relative temporal stability in the deep microbial assemblages in the PFB.

Temporal environment variation and beta diversity
No seasonal patterns were observed when comparing sampling cruises conducted at the end of the dry season (P-01 and P-03) and at the end of the rainy season (P-02 and P-04) (Figure 2). Nonetheless ANOVA test showed environmental differences between campaigns (Supplementary Table 2). These results were unexpected as it is well known that during the rainy season local rivers dominate the shelf and upper slope (Zavala-Hidalgo et al., 2003), while in the dry period Mississippi, Atchafalaya, and Texas rivers affect the shelf (Salmeroń-Garcıá et al., 2011). Deep sea is dominated by the Mississippi, Colorado and Rıó Grande rivers (Balsam and Beeson, 2003). The fluctuation in oceanographic processes that affect the circulation and distribution of the physicochemical properties may explain the apparent lack of seasonality in the sedimentary environments analyzed during this study as it has been reported elsewhere (Jochens and DiMarco, 2008).
Most sites (17 of 27) showed temporal stability in microbial assemblages since they were in the same cluster across campaigns. However, 10 sites showed changes in the community structure in different campaigns, mainly the transect D sites during P-02 and P-03 (Supplementary Figures 2, 3). The multivariable association analysis revealed positive and negative relations between the abundances of bacterial genera and specific environmental parameters that did exhibit temporal variation (Figures 5, 2; Supplementary Table 2). The complexity lies in that the parameters involved in structuring microbial communities (i.e., temperature, redox, salinity, nutrients, metals, hydrocarbons) changed differently between the bathymetric zones in each campaign (Figure 2; Supplementary Figure 1). This may be related to the fact that the sedimentary physicochemical properties that drove the changes in benthic prokaryote communities are closely linked with the hydrological and oceanographic dynamics of the study site according to its latitudinal and longitudinal location (Jochens and DiMarco, 2008;Hernańdez-Molina et al., 2016;Salcedo et al., 2017). Likely, the biochemical properties in surface sediments can show spatial and temporal variability ranging from millimeters to decimeters and seconds to seasons (Huettel et al., 2003). The combination of changes in community-structuring variables may partially explain a deep microbial assemblage being present on the shelf and a shelf microbial assemblage being present in the deep zone, as we observed in this study at 10 specific locations ( Figure 2B; Supplementary  Figure 2). In this regard, community assembly can be interpreted under contrasting perspectives such as deterministic niche paradigm and stochastic or neutral community models which are not necessarily exclusive (Pholchan et al., 2013). An understanding of the temporal and spatial scales over which these community changes occur might indicate underlying processes of microbial assemblages, and whether changes in microbial communities affect the functioning of the system or whether it is a mechanism to maintain ecosystem functions and services.
Environmental changes over time and space may provide an array of niches occupied by metabolically diverse microorganisms (Anderson et al., 2015;Yu et al., 2016). For instance, niche separation of ammonia-oxidizing ecotypes has been proposed by finding a dominant ecotype under oxygen-and ammonia-rich conditions, and a dominant ecotype in oxygen-and ammoniumpoor conditions (Muck et al., 2019). In the present study the presence of Thaumarchaeota, Nitrospinae, Planctomycetaceae, Nitrospira, and Nitrosomonadaceae (Figure 4; Supplementary  Figures 4, 6), provided evidence of nitrification and denitrification processes in the studied sites according to literature (Yu et al., 2016;Muck et al., 2019;Overholt et al., 2019). However, putative ammonia oxidizers related to Nitrospirota genera were associated to different physicochemical variables ( Figure 5) revealing a preference distribution, likely due to a niche differentiation. A genus related to the wb-A12 family (Nitrospira) was positively associated with depth, nitrite+nitrate, silicate, Cd, and Ni, suggesting a higher abundance toward slope and deep sites. wb-A12 has been related to nitrite oxidation (Orsi, 2018). Under oxic environmental conditions, such as the deep sites, nitrogen oxides commonly accumulate because microorganisms may perform ammonia oxidation to nitrite. Ammonia in oxic sediments is usually low because it can be used as an energy source when OM is scarce (Orsi, 2018), however, ammonium was not measured in this study. On the contrary, a genus from Nitrospiraceae was positively associated with temperature, salinity, lead, aluminum, and PAHs, suggesting its higher abundance in the shelf microbial assemblages ( Figure 5). Lower nitrite+nitrate concentration in these sites may be attributed to denitrification processes (Braker et al., 2012). These results propose the existence of different microorganisms likely being keystones in ecosystem functions under different conditions. The positive correlation between bacterial relative abundance and the concentration of silicate, nitrate plus nitrite, Ni, and Cd represented a feature in the deep zone ( Figure 5). It has been found that Ni and Cd may affect microbial growth (Cabrera et al., 2006). However, it has been proposed that in environments poor in organic matter (<2%), as the analyzed sediments ( Supplementary Figure 1), microorganisms can take advantage of the precipitation of minerals and metals that are necessary for microbial growth and metabolism. For instance, H 2 production by nickel-containing hydrogenases may be influenced by metal ions in the environment (Martins and Pereira, 2013). There seems to be a mutual influence between the precipitation of silica minerals/metals and microbial community composition which is influenced by environmental geochemistry (Sauro et al., 2018). The increase in metal concentration and inorganic cations has been found parallel to the silica minerals increase as in the present study ( Figure 2). Precipitation of silica minerals and metals is influenced by alkalinization and the production of extracellular polymeric substances (EPS) that may result from microbial processes such as sulfate reduction, inorganic nitrogen transformation, and MO degradation (Braissant et al., 2007;Sauro et al., 2018). In this respect, sediments in the deep zone showed an abundance of microorganisms able to perform nitrogen fixation to produce ammonia such as OM1 clade, nitrite oxidation such as wb1-A12, sulfate reduction such as Desulfurellaceae, and organic carbon degradation such as heterotrophic bacteria related to SAR324 (Dworking et al., 2006;Sheik, Jain, and Dick, 2014;Orsi et al., 2016;Orsi, 2018). Thus, it is possible that temporal changes in microbial community composition may affect the dynamic of silica minerals and metals in these environments.
Temperature, salinity, and PAHs can induce changes in microbial communities (Roussel et al., 2015;Chen et al., 2017;Chen et al., 2018;Shin et al., 2019) as was observed in the current study ( Figure 3B). Temperature changes may affect microbial activity, growth rates, and survival (Degerman et al., 2013;Orsi, 2018), higher salinity has been found to stimulate carbon metabolism and microbial processes such as sulfate reduction (Sivakumar et al., 2019), and PAHs can be degraded by aerobic and anaerobic bacteria or also inhibit microbial growth (Haritash and Kaushik, 2009). In marine sediments anaerobic degradation of PAHs by sulfate reducers has been demonstrated (Hamdan et al., 2017;Shin et al., 2019). This allows proposing that bacteria positively related with PAHs such as Desulfatiglans, Sva-0081, and Syntrophobacteraceae genera ( Figure 5), detected in high abundance in the shelf and slope may represent a biological mechanism to control the concentration of PAHs in the environment, that showed higher concentrations in this study (0.002 to 0.2 µg g -1 ) than reported (0.01-0.07 µg g -1 ) (Botello et al., 2015). Otherwise, PAHs could be toxic compounds (Haritash and Kaushik, 2009). In this regard, PAH concentrations seemed to affect the abundance of other microorganisms such as SAR324 genus which seemed to be better represented in deep sites ( Figure 5). SAR324 members have been linked to sulfur oxidation, nitrite reduction and alkane degradation, thus, these bacteria could control the concentration of simple hydrocarbons in the environment over time. In future works it would be interesting to address the effect of hydrocarbons in structuring benthic microbial communities under different environmental conditions and to explore the ecological relevance of benthic prokaryotes to control the concentration of different hydrocarbons in the PFB.

Conclusion
A two-year study at the end of the dry (P-01 and P-03) and the end of the rainy (P-02 and P-04) seasons was undertaken to compare the diversity, structure, and composition of benthic microbial communities to determine whether changes in benthic microbial assemblages are associated with spatial and temporal environment variation. The alpha diversity based on the Shannon index was similar in the bathymetric gradient over time. However, differences in the number of ASVs, suggested that higher species richness may reflect larger environmental changes. Species richness increase could be crucial to maintain ecosystem stability and ecosystem functions. Beta diversity analysis revealed that samples were grouped in three different clusters according to the depth gradient. The differences in the community structure between these clusters were related to the campaign, and changes in the water depth, temperature, salinity, redox, nitrate plus nitrite, silicate, Ni, Cd, Pb, Al, and aliphatic and aromatic hydrocarbons. Shelf and slope microbial assemblages showed high abundances of anaerobes related to sulfate-reduction, nitrate respiration, anaerobic degradation of OM, and the oxidation of PAHs. Deep sea microbial assemblages were integrated by members related to sulfur respiration, nitrogen fixation, and nitrite oxidation among other processes. Microbial composition in the shelf, slope, and deep sea showed relative temporal stability in most sites. However, environmental heterogeneity over time was complex and could provide different niches for ecotypes better adapted to specific environmental conditions, which could explain community structure changes in 10 out of 27 sites. Changes in benthic microbial communities may reflect the temporal fluctuation in the physicochemical properties of sedimentary environments, which may result from the prevailing oceanographic conditions according to the latitudinal and longitudinal location of the study site. Taken together, these results yield insights into the spatiotemporal dynamic of benthic bacterial communities in the PFB, nwGoM. However, long-term studies are needed to explore benthic microbial communities using water depth as an explanatory variable. This could allow the identification of sensitive bioindicators of environmental changes. Also, it is necessary to address in more detail the diversity of functional groups (including abundant and rare taxa) and their ecological relevance concerning biogeochemical cycles, nutrient supply, and metals and hydrocarbon transformations. An important issue is to understand whether changes in microbial communities affect system functioning or whether this is a mechanism for maintaining ecosystem functions and services under environmental variability.

Data availability statement
The original contributions presented in the study are publicly available. This data can be found here: NCBI, PRJNA885157.

Author contributions
MJ analyzed environmental data, interpreted bioinformatic results, wrote the manuscript, and accepted the final draft of the paper; DC-G performed bioinformatic analysis of the sequence reads, wrote sections of Materials and Methods, reviewed the manuscript, and accepted the final draft of the paper; MA-M and JG-M conceived the work, reviewed and contribute to the overall manuscript, and accepted the final draft of the paper. All authors contributed to the article and approved the submitted version,

Funding
This research was supported by the Mexican National Council for Science and Technology -CONACyT -Mexican Ministry of Energy-Hydrocarbon Fund, project 201441. CONACYT awarded MJ with a Ph.D. scholarship. This is a contribution of the Gulf of Mexico Research Consortium (CIGoM).