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 − μM L⁻¹), orthophosphate ( PO 4 3 − μM L⁻¹), and silicate minerals ( SiO 4 4 − μM L⁻¹) which were determined from the bottom water by colorimetric methods following standard protocols using a spectrophotometer UV-VIS (Agilent Technologies Cary 60) (Herguera-García 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 cut-off 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₃ 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).

Organic matter (OM%) was determined by oxidation with potassium dichromate and sand structure (as medium sands, MS%, fine sands, FS%, and very fine sands, VFS%) was analyzed weighing the particles retained in distinct mesh sizes (Hernández-Ávila et al., 2021). Heavy metal concentration (vanadium, V µg/g, nickel, Ni µg/g, cadmium, Cd µg/g, lead, Pb µg/g, and aluminum, Al µg/g) were quantified using ICP-MS (iCAP Q, Thermo Scientific), and hydrocarbons (total aliphatics µg/g and total aromatics, PAHs µg/g) were identified and quantified using a Clarus 500 PerkinElmer gas chromatography coupled to a mass spectrometer detector according to established protocols (Hernández-Ávila et al., 2021; Arcega-Cabrera et al., 2022).

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.

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/μL). 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 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.

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.

All the environmental data obtained from this study are presented in Supplementary Material 1 . Spatial differences of environmental variables (ANOVA, P<0.05) were observed between the bathymetric zones ( Supplementary Table 1 ). Over the four campaigns, the shallow zone (≤200 m of the water column) was characterized by the lower redox values (mean ± standard deviation, -136.8 ± 228.7 mV), the highest temperature (20.7 ± 6.1°C), salinity (39.2 ± 0.4), and the concentrations of Pb (26.1 ± 6 µg/g), Al (15193.6 ± 8009.2 µg/g), and PAHs (0.08 ± 0.04 µg/g), while it presented the lowest concentrations of NO 3 − + NO 2 − (4.3 ± 5.8 μM/L), PO 4 3 − (0.2 ± 0.2 μM/L), SiO 4 4 − (7.2 ± 5.2 μM/L), Ni (22.3 ± 4.9 µg/g), and Cd (0.2 ± 0.1 µg/g) ( Figure 2 ). The deep zone (>1200 m of the water column) presented the higher redox values (188 ± 329.9 mV), and the highest concentrations of NO 3 − + NO 2 − (20.4 ± 6.4 μM/L), SiO 4 4 − (27.5 ± 7.6 μM/L), Ni (32.6 ± 6.1 µg/g), and Cd (0.23 ± 0.1 µg/g), while presented the lowest temperature (6.2 ± 5.2°C), salinity (35 ± 0.01), and the concentrations of Pb (20.4 ± 4.6 µg/g), Al (8343.1 ± 6432.2 µg/g), and PAHs (0.04 ± 0.03 µg/g) ( Figure 2 ). The intermediate zone (201-1200 m of the water column) showed the highest PO 4 3 − concentration (2.3 ± 3.3 μM/L), it was similar to the shallow zone in terms of redox, Pb, and Al, and similar to the deep zone in terms of temperature, salinity, NO 3 − + NO 2 − , and Cd (ANOVA, P>0.05). The content of OM (0.9 ± 0.5%) and aliphatic hydrocarbons (1.2 ± 0.6 µg/g) did not show significant differences between zones (ANOVA, P>0.05). Overall, salinity, water depth, temperature, NO 3 − + NO 2 − , and SiO 4 4 − were among the five variables that contributed most to the spatial env ironmental variability ( Supplementary Table 1 ).

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 ). The highest concentrations of NO 3 − + NO 2 − and SiO 4 4 corresponded to intermediate and deep zones; however, NO 3 − + NO 2 − concentration decreased in P-02 in both zones, while SiO 4 4 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 ).

Community structure analysis grouped the samples into three clusters ( Figure 3A ) differing in terms of the water depth (R² = 0.56, F₂,₁₀₂ = 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).

Proteobacteria was the dominant phyla in all the sites among cruises. At this taxonomic level Acidobacteria, Chloroflexi, Planctomycetes, Nitrospirae (Nitrospirota), Gemmatimonadetes, Actinobacteria, Bacteroidetes, and Bathyarchaeota were detected in high abundances in the studied region ( Supplementary Figure 4 ). The most abundant microbial classes were represented by Deltaproteobacteria, Gammaproteobacteria, Alphaproteobacteria, Nitrospira, and Dehalococcoidia; although SAR202, Phycisphaere, Anaerolineae, and Gemmatimonadetes showed high relative abundances in microbial communities ( Supplementary Figure 5 ). At the family taxonomic level Syntrophobacteraceae, Rhodospirillaceae, Nitrospiraceae, OM1 clade, JTB255 marine benthic group (abbreviated JTB255 MBG), Desulfobacteraceae, Desulfobulbaceae, Anaerolinaceae, and Desulfurellaceae were among the most abundant groups detected across campaigns ( Supplementary Figure 6 ).

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 ).

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 ).

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 (Sánchez-Soto Jiménez et al., 2018; Ramírez 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 (Sánchez-Soto et al., 2018). A strong physical mixing and seasonality were suggested to explain higher diversity in the shelf sediments (Sánchez-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írez 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.

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; Sánchez-Soto Jiménez 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 (Sánchez-Soto Jiménez 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-Fernández 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 (Sánchez-Soto et al., 2018; Ramírez et al., 2020; Sánchez-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 (Sánchez-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.

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 (Salmerón-García et al., 2011). Deep sea is dominated by the Mississippi, Colorado and Río 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; Hernández-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 ammonium-poor 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₂ 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.