The sampling was carried out seasonally in 2019 (Supplementary Table 1) at five sites of the Venice lagoon, Italy (Figure 1): Chioggia (C), Marghera (M), Palude della Rosa, (P), Sacca Sessola (S), and Tresse (T). Furthermore, C, M, and P were sampled with the same time frequency in 2020. The sites were chosen to cover the four sub-basins, defined according to risk analyses following the international sediment quality guideline of the Lagoon (Apitz et al., 2007; Supplementary Table 2). They are located in the southern (C), in the central (S), in the central-north (M and T), and in the northern (P) sub-basins. The Chioggia site is near the city of Chioggia, very close to the southern lagoon inlet and more subjected to marine influence. The Sacca Sessola site, situated in the center of the Lagoon, south of Venice city, is potentially affected by urban contamination (Picone et al., 2016). The Marghera site is inside the active industrial area of Porto Marghera and is characterized by high anthropogenic impact. The Tresse site is next to an artificial island created as a dumping site for urban and industrial waste. Finally, Palude della Rosa site is less impacted by human activities (Picone et al., 2016).

In Chioggia, Marghera, and Palude della Rosa samples were collected with a Van Veen grab (0.045 m²), while in Sacca Sessola and Tresse with manually operated corers (10 cm diameter). The first centimeter of sediment was collected with a sterile spatula. The sampling was performed in three independent replicates. Water temperature and salinity were recorded with the multiparametric probes HI98128 (HANNA Instruments, Italy) and AP-2000 (Aquasearch, Italy) (Supplementary Table 3).

The water content was estimated by weighing fresh and dried (105°C for 24 h) sediment aliquots (2 cm³). The grain size composition was determined with a laser Particle Size Analyzer (BECKMAN COULTER LS 13 320) after treating 1 g of wet sediment with H₂O₂. The sediment was described using Shepard’s classification (Shepard, 1954) and the data were expressed as percentages of sand, silt, and clay. The limit between silt and clay was fixed to 4 μm, following the Udden-Wentworth classification (Wentworth, 1922).

For Total Organic Carbon (TOC) estimation, freeze-dried sediment was grounded in a ceramic mortar and sieved on a 250 μm iron steel sieve (Endecotts LTD., United Kingdom). Triplicate subsamples of about 8–12 mg were weighed on a micro ultrabalance with an accuracy of 0.1 μg. Before TOC determination, subsamples were treated directly into capsules with increasing concentrations of HCl (0.1 N and 1 N) to remove the carbonate fraction (Nieuwenhuize et al., 1994). Carbon content was determined using a CHNO-S elemental analyzer ECS 4010, Costech, Italy) according to Pella and Colombo (1973). Standard acetanilide (Costech, purity ≥ 99.5%) was used to calibrate the instrument and empty capsules were also analyzed to correct for blank. Measurement quality control was performed using internal standards and it was also verified for carbon against the certified marine sediment reference material PACS-2 (National Research Council Canada). Significant differences in TOC were calculated in the R environment (v. 4.0.3, R Core Team, 2019) by Kruskal-Wallis and Wilcoxon non-paired tests. False discovery rate (FDR) was used for p-value correction. Results with q-value < 0.05 were considered significant.

The prokaryotic abundance was estimated by flow cytometry following Deng et al. (2019) with modifications, using a FACSCanto II (Becton Dickinson) instrument equipped with an air-cooled laser at 488 nm and standard filter setup. A sediment slurry was subjected to mechanical shaking, ultrasonication, and hydrofluoridric acid (1%, final concentration) treatment. Slurry mixed with the Stop solution was filtered through a 10 μm polycarbonate filter, diluted, and stained with SYBRGreen I (2×, final concentration). Negative controls were prepared by autoclaving sediments and running unstained sediments. The acquisition threshold was set to green fluorescence and stained prokaryotic cells were identified in the side scatter versus green fluorescence plot. The flow rate was calibrated daily, by running distilled water and weighing it before and after the run (at least five replicates). Data were acquired and processed with the FACSDiva software (Becton Dickinson). Prokaryotic cells per gram of dry sediment were calculated using the acquired cell counts, the flow rates and the respective sample water content, and statistical analyses were run as described for TOC data.

A Principal component analysis (PCA) was performed to explore and visualize similarities among the samples based on environmental variables (temperature, salinity, grain size, TOC and prokaryotic abundance) with the R package FactoMineR (Le et al., 2008), factoextra (Kassambara, 2015), and ggplot2 (Wickham et al., 2016).

DNA was extracted with DNeasy PowerSoil Pro Kit (Qiagen) following the manufacturer’s instructions. Quality and quantity of the extracted DNA was assessed with Nanodrop spectrometer (Thermo Fisher Scientific) and with Qubit Fluorimeter (Thermo Fisher Scientific).

For the amplicon sequencing analysis, sediments from both 2019 (at each site, one sample per season for a total of 20, in three replicates) and 2020 (at C, M, and P sites, one sample per season for a total of 12, in three replicates) were considered. The V4-V5 region of 16S rRNA gene was amplified using 515-Y (5′-GTGYCAGCMGCCGCGGTAA-3′) and 926R (5′-CGYCAATTYMTTTRAGTTT-3′) primers (Parada et al., 2016). Libraries were prepared following the 16S Metagenomic Sequencing Library Preparation protocol and run on an Illumina MiSeq System for a read length of 2 × 250 bp at the genetic and epigenetic ARGO Open Lab Platform, Area Science Park, Trieste, Italy.

For the metagenomic analysis, sediments from 2019 (at each site, one sample per season for a total of 20, in three replicates) were considered. Libraries were prepared following the Illumina Nextera DNA Flex Library Prep protocol and run on an Illumina NovaSeq 6000 System for a read length of 2 × 300 bp at the genetic and epigenetic ARGO Open Lab Platform, Area Science Park, Trieste, Italy.

Bioinformatic analyses were performed with QIIME2 2020.6 (Bolyen et al., 2019). Raw sequences were quality filtered and denoised with DADA2 (Callahan et al., 2016). Alpha-diversity metrics were estimated after samples were rarefied. To assess the reliability and robustness of the replicates, a correlation analysis (Spearman’s ρ) on the amplicon sequence variants (ASVs) table was performed with the R package ggplot2 (Wickham et al., 2016). Taxonomy was assigned to ASVs using the sklearn naïve Bayes taxonomy classifier (Bokulich et al., 2018) against the Silva 138 99% reference database with 7-level taxonomy (Quast et al., 2013). Reads belonging to Eukarya, mitochondria, chloroplast, and with frequency < 2 (singletons) were removed.

Statistical analyses were conducted in the R environment (v. 4.0.3, R Core Team, 2019). Rarefaction curves on ASVs table were plotted with the R package vegan (Oksanen et al., 2019). Significant differences in alpha diversity metrics were calculated by Kruskal-Wallis H test. FDR was used for p-value correction. Results with q-value < 0.05 were considered significant.

To visualize similarity patterns of prokaryotic communities, a Principal Coordinates Analysis (PcoA) using Bray-Curtis dissimilarity matrices was constructed using normalized read abundances in the R package vegan (Oksanen et al., 2019). A permutational multivariate analysis of variance (PERMANOVA) with 4999 permutations was computed on normalized ASVs table to investigate the effect of site and season using the function adonis in the R package vegan (Oksanen et al., 2019).

To further explore the environmental drivers of the community, we used a distance-based redundancy analysis (dbRDA) using the capscale function in the R package vegan (Oksanen et al., 2019). The environmental variables (temperature, salinity, grain size, and TOC) were normalized using a z-score transformation and variance partitioning analysis was used to quantify the contribution made by each variable using the varpart function in the R package vegan (Oksanen et al., 2019). The significance of the variables was tested with ANOVA and 4999 permutations after building a reduced model. As the grain size components (sand, silt, and clay) present multicollinearity, only one component (sand) was included in the analysis.

An indicator species analysis was performed at site level with the R package indicspecies (De Cáceres and Legendre, 2009). The results were visualized as networks by mean of the R package igraph (Csard and Nepusz, 2006). Moreover, an analysis of composition of microbiomes (ANCOM; Mandal et al., 2015) was used to test for differences at the genus level. Microbial diversity visualization at different taxonomic level was done using the average normalized abundance of the three replicates of each sample, with R package phyloseq (McMurdie and Holmes, 2012) and ggplot2 (Wickham et al., 2016). Microbial signatures of fecal- and sewage-associated bacteria (McLellan et al., 2010; Newton et al., 2013; Paliaga et al., 2017) were also evaluated. Statistical Analyses of Metagenomic Profiles 2.1.3 (STAMP; Parks and Beiko, 2010) was used to assess differences in the relative proportion of sequences belonging to the different taxa considering the whole dataset. Significant differences were calculated by ANOVA and Tukey-Kramer post hoc test. FDR was used for p-value correction. Results with q-value < 0.001 were considered significant.

The number and quality of reads were checked with Fast QC (Andrews, 2010). Reads were assembled into contigs using MEGAHIT-1.2.9 (Li et al., 2015) combining the three replicates of each sample after testing the significance of their correlations. Gene prediction was performed with Prodigal 2.6.3 (Hyatt et al., 2010) on contigs longer than 1,000 bp (Van der Walt et al., 2017).

The contigs > 1,000 bp were submitted to Metagenome Analysis Rapid Annotation using Subsystem Technology (MG-RAST; Keegan et al., 2016) for a similarity search using the SEED Subsystem (Overbeek et al., 2005) and the Kyoto Encyclopedia of Genes and Genomes (KEGG; Kanehisa and Goto, 2000) databases, using a cut-off E-value of 1e^–5, minimum identity of 60% and a minimum alignment length of 15 bp. SEED and KEGG annotations were normalized by the number of reads assigned to the prokaryotic single copy gene recA (Acinas et al., 2021) prior to calculating the relative abundances.

Key oceanic marker genes (Sunagawa et al., 2015) and key sediment metabolic genes (Dombrowski et al., 2018) were evaluated to provide an overview of the main biochemical cycling metabolisms present. We have created a novel ad hoc specific microscale marker function list to explore the microbial potential at the small-scale in the sediment environment.

Annotation results were visualized using Statistical Analyses of Metagenomic Profiles 2.1.3 (STAMP; Parks and Beiko, 2010) to assess differences in their relative proportion considering the whole dataset. Significant differences among samples were calculated by ANOVA and Tukey-Kramer post hoc test. FDR was used for p-value correction. Results with q-value < 0.001 were considered significant.

The five stations were characterized by different percentages of the grain size classes of sand, silt, and clay (Supplementary Table 4), in agreement with previous studies (Borin et al., 2009; Picone et al., 2016; Zonta et al., 2018). Based on Shepard’s classification, Chioggia was categorized as sand, Marghera as clayey silt, Palude della Rosa as silty-loam, Sacca Sessola as loam, and Tresse as silty sand.

Total Organic Carbon ranged from 0.87 to 32.8 mg g_d^–1 (Supplementary Table 5). Marghera (16.4 ± 7.5) and Palude della Rosa (14.8 ± 7.5) showed significantly higher concentrations than Sacca Sessola (6.6 ± 0.9), Tresse (5.5 ± 1.0), and Chioggia (3.3 ± 3.2). No significant differences were detected among seasons.

The average sediment prokaryotic abundance was 25.7 ± 21.5 × 10⁹ Cells g_d^–1 (Supplementary Figure 1). Palude della Rosa showed significantly higher abundances (45.4 ± 32.5 × 10⁹ Cells g_d^–1) than Marghera (27.7 ± 24.6 × 10⁹ Cells g_d^–1; p < 0.05), Chioggia (8.3 ± 7.5 × 10⁹ Cells g_d^–1; p < 0.001), Sacca Sessola (7.8 ± 4.3 × 10⁹ Cells g_d^–1; p < 0.001), and Tresse (4.5 ± 4.1 × 10⁹ Cells g_d^–1; p < 0.01). Overall, Autumn data (30.2 ± 26.9 × 10⁹ Cells g_d^–1) were significantly higher than Spring (10.1 ± 8.5 × 10⁹ Cells g_d^–1; p < 0.05) and Summer ones (16.0 ± 12.5 × 10⁹ Cells g_d^–1; p < 0.05). The PCA (Figure 2) showed a separation of sampling sites on PC1 (54.35% variance), mainly due to grain size features and TOC concentration. The effect of seasonality was also observed, even though the projection of temperature on PC2 explained a lower percentage of the total variance (18.16%).

A total of 31,189,504 raw sequences were generated. After the denoising procedure, 20,191,333 were retained with an average of 219,471 ± 92,086 sequences per sample. The average number of ASVs was 5,230 ± 1,536. The rarefaction curves indicated that the sequencing effort was enough to assess the sample biodiversity (Supplementary Figure 2). The replicates significantly correlated (p < 0.05) with an average ρ of 0.86 ± 0.09.

Shannon’s diversity index showed an average of 10.7 ± 0.4 and Pielou’s evenness showed an average of 0.89 ± 0.02 (Supplementary Figure 3). Overall, Autumn had a significantly higher index than the other seasons. No significant differences in biodiversity indices were detected among sites and site-season combinations.

The PERMANOVA highlighted that site, season, and their interaction had a significant effect (p < 0.001) on the prokaryotic community structure (Supplementary Table 6). The highest amount of variance was explained by the site (37%), followed by season-site interaction (11%) and season (6%).

The PCoA (Supplementary Figure 4), in agreement with the PERMANOVA, grouped the samples by site.

The variance partitioning analysis showed that environmental parameters explained ∼20% of the variance (Supplementary Table 7). The dbRDA (Figure 3) revealed that grain size, TOC, and salinity had a significant role (p < 0.001) in shaping the communities, while temperature was not a significant variable (p > 0.05).

Regarding the taxonomic composition, at the domain level most reads were assigned to Bacteria (average on the overall dataset 96.1 ± 2.6%), and a smaller percentage (3.8 ± 2.6%) to Archaea.

For Bacteria (Figure 4), the most represented phyla were Proteobacteria (29.3 ± 4.4%), followed by Bacteroidota (15.5 ± 4.4%) (Figure 4). Desulfobacterota were significantly less abundant in Chioggia than the other sites (9.2 ± 2.4% vs. 12.8 ± 0.8%), while Actinobacteriota presented the opposite distribution pattern (4.8 ± 1.2% vs. 2.2 ± 0.5%) (Figure 4). The most represented classes were Gammaproteobacteria (23.9 ± 3.7%), followed by Bacteroida (14.7 ± 4.3%) (Supplementary Figure 5). The most represented orders were the bacterial Flavobacteriales (6.8 ± 3.8%), Desulfobacterales (6.1 ± 2.2%), and Pirellulales (5.2 ± 2.3%) (Supplementary Figure 5). The most represented families were Flavobacteriaceae (6.2 ± 3.6%), Desulfosarcinaceae (5.2 ± 1.8%), and Pirellulaceae (4.5 ± 0.9%). The most abundant genus was Woeseia (4.4 ± 0.9%) (Figure 4).

For Archaea (Figure 4), the most represented phyla were Crenarchaeota (1.3 ± 0.6%), followed by Thermoplasmatota that were lower in Chioggia in respect to the other sites (0.14 ± 0.08% vs. 0.9 ± 0.6%). The ammonia-oxidizing archaeon Candidatus Nitrosopumilus was significantly higher in Chioggia (1.2 ± 0.7% vs. 0.1 ± 0.1%) (Figure 4).

Among the genera related to microscale processes, we detected the cable bacteria Candidatus Electrothrix at all sites (0.1 ± 0.2%); these bacteria have long, multicellular filaments that can conduct electric currents over centimeter-scale distances (Trojan et al., 2016). Likewise, lignin-degrader genera (Bacillus, Novosphingobium, Pseudomonas, Rhodococcus, and Streptomyces; Lee et al., 2019) were found at all sites (0.006 ± 0.008%). The Indicator Species analysis showed that each site hosted a set of indicator bacterial genera (Figure 5, Supplementary Table 8), with Chioggia showing the highest number (n = 21). On the contrary, excluding Chioggia, the highest number of genera (n = 27) were shared among Marghera, Palude della Rosa, Sacca Sessola and Tresse. Chioggia was characterized by the presence of aerobic (or strictly aerobic), halophilic bacteria such as Aquimarina, Aureisphaera, and Pseudoalteromonas and by the absence of genera of sulfate-reducing Desulfobacteraceae that were shared among the other sites, like Desulfoconvexum, Desulfonema, and Desulfobacter. ANCOM analysis (Supplementary Figure 6) confirmed the significant higher abundance ammonia-oxidizing archaeon Candidatus Nitrosopomilus in Chioggia and a lower presence of the Desulfobacteraceae member Desulfobacter. Moreover, it highlighted a higher presence of the nitrite-oxidizing Nitrospira and of the aerobic Ahrensia and Granulosicoccus. Palude della Rosa showed significantly higher presence of the halophile Halochromatium and the glycogen-accumulating Candidatus Competibacter.

Microbial signatures of selected anthropogenic-related taxa (Figure 6) presented overall 0.1–0.3% relative abundance for fecal-associated bacteria and 0.01–0.04% for sewage-associated bacteria. In Sacca Sessola and Chioggia, fecal-associated microbes showed the highest presence, mainly due to Clostridiaceae, while Marghera showed the highest proportion of Bacteroidaceae and Ruminococcaceae. For the sewage-associated bacteria, Marghera showed on average the highest proportion of reads due to Acinetobacter and Trichococcus while Palude della Rosa showed the lowest proportion.

A total ∼120 million of raw reads (30 Gbp) were generated, specifically 14.8 ± 3.7 millions of paired-end reads with quality score of 35 ± 2 for each metagenome. On average, contigs derived from the three co-assembled replicates for each sample were 5,208,341 ± 376,210 contigs > 1,000 bp, and 652,755 ± 174,565 coding sequences.

In the KEGG functional characterization (Supplementary Figure 7), the most represented level 1 categories were Metabolism (63.2 ± 0.4% on the overall dataset), Genetic Information Processing (17.0 ± 0.2%), and Environmental Information Processing (14.2 ± 0.5%); at level 2, Amino acid metabolism (21.5 ± 0.6%), Carbohydrate metabolism (12.9 ± 0.3%), and Translation (7.9 ± 0.1%); at level 3, ATP-binding cassette (ABC) transporters (6.3 ± 0.4%), Two-component system (4.7 ± 0.2%), and Aminoacyl-tRNA biosynthesis (4.2 ± 0.1%).

Key oceanic marker genes (Sunagawa et al., 2015) and key sediment metabolic genes (Dombrowski et al., 2018) have been selected and clustered as follows: Carbon fixation, Energy and metabolism, and Carbon-nutrient coupling.

The category of Carbon fixation included: oxygenic photosynthesis, anoxygenic photosynthesis, and other autotrophic pathways (Wood-Ljungdahl, Hydroxypropionate bicycle, hydroxypropionate/hydroxybutyrate, C1-metabolism). The category of Energy and metabolism included: aerobic respiration, fermentation, manganese-related metabolism, nitrogen metabolism, iron-related metabolism, sulfur metabolism, methanogenesis, butane metabolism, and hydrogenase. The category of Carbon-nutrient coupling included: phosphorus metabolism, peptide degradation, motility and chemotaxis, transport, and fatty acids degradation. We also assessed Anthropogenic-related metabolisms regarding hydrocarbons and xenobiotics (Figure 7, Supplementary Table 9).

Overall, in the sediment microbial environment, Peptide degradation was the most represented category (4.5 ± 0.1%), followed by Transport (2.7 ± 0.1%), Aerobic respiration (1.8 ± 0.1), Other autotrophic pathways (1.7 ± 0.01%), and Motility and chemotaxis (1.3 ± 0.1%). Here we present the significant differences detected in the above-mentioned categories across the sites (Figure 7). Within Carbon fixation, Oxygenic photosynthesis was higher in Chioggia (0.93 ± 0.02% vs. 0.86 ± 0.01%), while Anoxygenic photosynthesis was higher in Palude della Rosa (0.008 ± 0.001% vs. 0.003 ± 0.002%).

Within Energy and metabolism, Chioggia presented a lower number of genes coding for functions related to Fermentation (0.50 ± 0.03% vs. 0.72 ± 0.03%), Sulfur metabolism (0.62 ± 0.01% vs. 0.72 ± 0.03%), and Butane metabolism (0.66 ± 0.01% vs. 0.58 ± 0.03%) in respect to the other sites. Marghera was higher in Iron-related metabolism (1.13 ± 0.03% vs. 0.97 ± 0.03%) whereas Palude della Rosa, Tresse and Marghera harbored the highest number of Hydrogenase genes (0.190 ± 0.005% vs. 0.147 ± 0.01%). Within the Carbon-nutrient coupling category, Phosphorus metabolism functions were less present in Marghera and Tresse (0.582 ± 0.007% and 0.692 ± 0.008%, respectively). Regarding Nitrogen metabolism, we detected significant amounts of genes coding for ammonia monooxygenase subunits (amoA, amoB, and amoC) in Chioggia in respect to the other sites (0.0007 ± 0.0003% and 0.00008 ± 0.00002%, respectively), in accordance with the presence of Candidatus Nitrosopumilus (Walker et al., 2010; Bayer et al., 2019). Chioggia presented the highest number of genes related to Fatty acid degradation (0.713 ± 0.001% vs. 0.662 ± 0.002%). Furthermore, in the Anthropogenic-related metabolisms category: Hydrocarbon degradation was significantly higher in Chioggia (0.183 ± 0.009% vs. 0.118 ± 0.003%), mostly due to hydrolases involved in the degradation of chlorocyclohexane, chlorobenzene fluorobenzoate, including toluene. The Xenobiotic metabolism and degradation in the same site was also significantly higher (0.524 ± 0.051% vs. 0.410 ± 0.002%), mostly due to haloalkane dehalogenase (dhaA).

This sandy site was significantly different from the others for Motility and chemotaxis (1.5 ± 0.2% vs. 1.3 ± 0.3%, flagellum-related genes), and Transport (2.5 ± 0.1% vs. 2.7 ± 0.1%, lipid and vitamin).

In the SEED functional characterization (Supplementary Figure 8), the most represented level 1 categories were Clustering-based subsystems (14.7 ± 0.1%; functional coupling evidence with unknown function), Carbohydrates (11.0 ± 0.3%), Amino Acids and Derivatives (8.8 ± 0.1%), and Protein metabolism (8.6 ± 0.1%), whereas at the level 2, Plant-Prokaryote DOE project (6.8 ± 0.1%), Protein biosynthesis (5.3 ± 0.1%), and Central carbohydrate metabolism (3.7 ± 0.1%).

Regarding anthropogenic-related functions, the “Virulence, Disease and Defense” category was significantly higher in Marghera and Tresse in respect to the other sites (3.36 ± 0.02% vs. 3.04 ± 0.03%), mainly due to a significantly higher presence of genes in the “Resistance to antibiotics and toxic compounds” group (3.08 ± 0.04% vs. 2.72 ± 0.03%).

Antibiotic resistance annotation (Figure 8) showed that Multidrug and Beta-lactamase resistance were the most present at all sites (0.93 ± 0.01% and 0.58 ± 0.03%, respectively), with Vancomycin resistance being significantly higher in Chioggia (0.016 ± 0.002% vs. 0.007 ± 0.002%).

Among Toxic compound resistance annotation, Cobalt-zinc-cadmium was the most present at all sites (0.676 ± 0.002%), with Mercury resistance being significantly higher in Marghera and Tresse (0.009 ± 0.001% vs. 0.004 ± 0.002), and Arsenic resistance being significantly higher in Chioggia and Palude della Rosa (0.221 ± 0.002% vs. 0.196 ± 0.004).

The key specific microscale marker functions, identified with the SEED database, have been clustered as follows: DNA metabolism (CRISPs and DNA uptake and competence), Mobilome (Gene Transfer Agents, integrons, pathogenicity islands, phages and prophages, plasmids, and transposable elements), Regulation and Cell Signaling (biofilm formation, phenazine biosynthesis, osmotic stress, sporulation, siderophores, toxic-antitoxic systems, and two-component systems) and Nanomachines (protein secretion system Type IV and Type VI) (Figure 9).

Overall, in the sediment microbial environment, Phages and prophages was the most represented category (0.58 ± 0.07%), followed by Osmotic stress (0.27 ± 0.01%), Pathogenicity islands (0.26 ± 0.01), Toxic-antitoxic systems (0.13 ± 0.02%), Transposable elements (0.12 ± 0.09%), and Two-component systems (0.01 ± 0.09%) (Figure 9).

Here we present the significant differences detected in the above-mentioned categories across the sites. Within Mobilome, Phages and prophages and Transposable elements were lower in Chioggia (0.47 ± 0.02% vs. 0.62 ± 0.04% and 0.1 ± 0.01% vs. 0.12 ± 0.01%, respectively). Within Regulation and Cell Signaling, Biofilm formation was higher in Marghera, Sacca Sessola, and Tresse (0.04 ± 0.01% vs. 0.02 ± 0.01%), Osmotic stress was higher in Chioggia, Palude della Rosa, and Sacca Sessola (0.28 ± 0.01% vs. 0.25 ± 0.01%), and Toxic-antitoxic systems was lower in Chioggia (0.08 ± 0.01% vs. 0.14 ± 0.01%). Within Nanomachines, Protein secretion system type VI was lower in Chioggia (0.03 ± 0.01% vs. 0.06 ± 0.01%).

Microbial communities in marine sediment are extremely diverse due to the combination of different factors such as habitat temporal stability, high niche diversity, and resource partitioning due to the complex physical-chemical gradients (Zinger et al., 2011; Acosta-Gonzalez and Marques, 2016). As such, they represent a major reservoir of genetic variability in the marine environment (Polymenakou et al., 2005).

Our analysis detected a pronounced diversity in the communities, in agreement with 16S studies in the Venice Lagoon (Borin et al., 2009; Quero et al., 2017) and in other transitional and coastal benthic systems (Bianchelli et al., 2020; Miksch et al., 2021). The sampling site specific features (grain size, salinity and TOC content, Figure 2) were important drivers in shaping the prokaryotic assemblages, while seasonality had a significant but minor role (Figure 3). Depending on the area, seasonality in the sediment microbial communities is not as regularly and clearly found as in the water column (Tšertova et al., 2011; Gobet et al., 2012; Miksch et al., 2021), even in the Venice Lagoon (Quero et al., 2017). Quero et al. (2017), in a study including pelagic and benthic communities in the central part of the Venice Lagoon, found seasonality in both compartments but less pronounced in surface sediments than in the overlaying water, where this feature is well documented (Celussi et al., 2009).

In marine sediments, it is now known that physical-chemical characteristics including salinity, chlorophyll a, TOC content, grain size, and benthic fauna can influence the microbial communities (Sapp et al., 2010; Pala et al., 2018; Hoshino et al., 2020). Generally, the composition of the sediment bacterial communities can present high dissimilarities between sites, reflecting the heterogeneity of these environments, mostly due to limited physical mixing and the geographical distance (Acosta-Gonzalez and Marques, 2016).

We found that salinity, grain size, and TOC significantly contributed to the structure of the microbial communities in the different sites of the Venice Lagoon (Figure 3).

Salinity is recognized as a major environmental determinant of microbial composition in benthic communities (Lozupone and Knight, 2007; Bolhuis et al., 2013) in coastal (Bolhuis and Stal, 2011; Severin et al., 2012; Bolhuis et al., 2013) and lagoon systems (Tsuboi et al., 2013; Pavloudi et al., 2016; Behera et al., 2017). In our study, the more marine-influenced site (Chioggia, located at the southern inlet of the Venice Lagoon) was confirmed as the most diverse one, forming a separated group in the multivariate analyses (see Figure 3).

It is known that grain size can affect microbial abundance and activity by providing different surface properties (e.g., fine sediments provide greater areas that can be colonized) and influencing permeability, water chemistry, hydraulic conductivity, and pore space between sediment particles (Dale, 1974; Zhang et al., 1998). A positive correlation between bacterial abundance, sediment grain size, and organic carbon content is well established (Dale, 1974) and reported in different studies (Santmire and Leff, 2007; Legg et al., 2012; Tsuboi et al., 2013; Fazi et al., 2020). In another Italian Adriatic lagoon, it was found that lowest salinity, higher percentage of fine grain sediment and organic matter corresponded to the highest microbial abundances (Cibic et al., 2019). This was confirmed also in our study, in which the sites enriched in silt, Marghera and Palude della Rosa, presented the highest prokaryotic abundance and TOC concentration.

Sediment grain type influences indirectly the microbial assemblage at the microscale, due to the variable water fluxes (Ahmerkamp et al., 2020) that changes the redox conditions thus driving more diverse and able-to-adapt communities in sandy sites compared to the clay and silty ones.

Being and important source of energy for heterotrophs the organic carbon (TOC) pool is an important structuring agent for benthic microbial communities (Oni et al., 2015). As such, in our study the content of TOC was important in structuring the differences among the communities (e.g., Chioggia vs. other sites), even though we did not analyze its quality (biopolymers’ concentration, labile vs. recalcitrant matrices, etc.), which is also known to affect the composition of bacterial communities (e.g., Fazi et al., 2020). Likewise, ecological factors that were not investigated in this study can contribute to the differentiation of the sites: interspecific bacteria competition, grazing pressure, viral lysis, vegetation, contaminants, and also stochastic events (Böer et al., 2009; Luna et al., 2013; Acosta-Gonzalez and Marques, 2016; Bianchelli et al., 2020).

The microbial communities were characterized by few dominant taxonomic groups, the core microbiome, and a high number of low abundant taxa, the so-called rare biosphere (Sogin et al., 2006), that accounted for most of the phylogenetic diversity.

For instance, among the 84 phyla detected in the whole dataset, 49 were considered rare, as represented by an average relative abundance < 0.01% (Galand et al., 2009; Gobet et al., 2012), as well as for 1181 out of the 1587 detected genera. The presence of such rare biosphere is of ecological importance as these taxa, often overlooked, can have an over-proportional role in biogeochemical cycles and can represent a hidden driver of ecosystem functions (Jousset et al., 2017).

Even if we detected significant differences among sites in high abundant taxa, from the phylum to the genus level, the numerous low-abundant taxonomic units were those that more contributed to their site-specific fingerprint (see Figures 3, 5). These taxa were mostly uncharacterized (e.g., belong to uncultured bacteria typically detected only with molecular techniques and are still not defined in terms of physiology and metabolism), and more studies are needed to define their ecological roles.

In accordance with Jousset et al. (2017) we suggest that the rare biosphere in the Venice Lagoon can be considered as an ensemble of keystone species that provide stability and seeds for ecosystem functions, specifically in terms of water purification, nutrients, and global climate regulations (da Mosto et al., 2020).

Regarding the most abundant taxa, the globally distributed (Gamma)Proteobacteria were confirmed to be prevalent at all sites (see Figure 4), as generally detected in sediment in both the costal and the deep seafloor (Fry et al., 2008; Zinger et al., 2011; Acosta-Gonzalez and Marques, 2016), including the Venice Lagoon one (Borin et al., 2009; Quero et al., 2017). Among these, the ubiquitous genus Woeseia was the most abundant (see Figure 4). Recent molecular characterization suggested that it is likely to grow on proteinaceous matter, potentially derived from detrital cell membranes, cell walls, and other organic remnants (Mußmann et al., 2017; Hoffmann et al., 2020).

Several significant differences in taxonomic composition were highlighted by ANOVA/Tukey-Kramer and ANCOM analyses in sulfur-related bacteria (Desulfobacteria, Desulfobacterales) involved in sulfate reduction. In our study, these taxa were less represented in Chioggia where, differently form the inner area of the Lagoon where the sediment is mostly anoxic, hypoxic conditions are commonly found (Borin et al., 2009). Chioggia sandy sediment is in fact more exposed to the effect of winds, tidal currents, and ship movements, that favor partial oxygenation (Borin et al., 2009). Permeable sandy sediment pore space is constantly flushed by the overlying water and could trap detritus and living cells from the water column (Boudreau et al., 2001; Gobet et al., 2012). In the other sites, the higher amount of organic matter contributes to the oxygen depletion and to the activity of sulfate-reducing bacteria in sulfide production (Zaggia et al., 2007). This is further supported by the metagenomics analysis, which also highlighted a significantly lower number of genes related to sulfur metabolism in Chioggia (see section “Wide Range of Carbon Fixation and Carbon Flow Strategies”).

Sediment bacteria play a key ecological and biogeochemical role in marine ecosystems. They present high abundance and genetic variability, critical role in the transformation and speciation of major bioactive elements (e.g., carbon, nitrogen, phosphorus, oxygen, and sulfur), and in the degradation of organic pollutants (Polymenakou et al., 2005; Wu et al., 2008). The metagenomics approach adopted in this study allowed the identification of the genes involved in biogeochemical cycles providing insights into the potential metabolic functioning of the microbial communities (Doney et al., 2004; Dombrowski et al., 2018; Acinas et al., 2021). Following the “carbon currency” in the microbial sediment compartment we were able to map the carbon fixation processes, the energy and metabolism processes, the carbon-nutrient coupling, the carbon reworking processes, which includes the degradative processes, and the behavior and anthropogenic-related metabolisms (see Figure 7).

In our sites, we found that oxygenic photosynthesis (by Prokaryotes) was the main carbon fixation pathway. Anoxygenic photosynthesis was also detected together with five additional carbon fixation pathways which differ in reducing compounds, energy source, and oxygen sensitivity of enzymes alternative pathways (Wood-Ljungdahl, Hydroxypropionate bicycle, hydroxypropionate/hydroxybutyrate, C1-metabolism). These findings highlighted the metabolic flexibility of the community in fixing CO₂ from different sources thus suggesting high diversity in C fixation modes in the sediment. Once the carbon is turned into biomass and organic carbon, it can be metabolized via aerobic and anaerobic respiration or fermentation. The surface sediment is affected by the light diurnal cycle and by the presence of oxygen (as a by-product of photosynthesis and from the oxygenated seawater column) (Petersen et al., 1994; Nealson, 1997). Aerobic respiration was the main pathway to respire the organic matter and secondly, anaerobic respiration (N, Mg, Fe, S; i.e., changing the terminal electron acceptor according to their availably, Nealson, 1997) was detected. Furthermore, often understudied metabolisms, such as methanogenesis, butane, proton degradation (hydrogenase, with a great diversity of enzymes), and fermentation were also important in enriching the metabolic repertoire of the microbial communities across sites (Falkowski et al., 2008). Within the concept of shared metabolisms, cable bacteria (family Desulfobulbaceae, Kjeldsen et al., 2019) were present in our samples. These microbes are formed by centimeter long filaments that span from the oxic water-sediment interface to the sulphidic deeper sediment layer. Within their filament, they couple sulfide oxidation with oxygen or nitrate reduction via long-distance electron transport.

Sulfur-related genes were found in high abundance in our sites with the exception of the sandy site of Chioggia. This is possibly due to the sediment texture that is more permeable to water and oxygen thus causing the community metabolism to switch more often from anoxic to oxic functionalities. Sulfur is essential for building amino acids (e.g., methionine and cysteine), cellular components (sulfolipids) and metabolites (DMSP) (Wasmund et al., 2017). It has a paramount role as electron acceptor (e.g., sulfate and thiosulphate) and electron donor (sulfite and elemental sulfur) (Jørgensen et al., 2019). We found genes related to sulfate-reduction and sulfur-oxidation processes that are coupled in the surface sediment. Overall, sulfur cycling is interconnected with the other elements, the availability of organic matter and redox species present in the lagoon system. In a high-temperature and low-oxygen lagoon scenario, sulfur-related metabolisms might become dominant in the ecosystem functioning.

Phosphorus-cycle-related genes were also highly abundant in our samples, being especially linked to the involvement of nucleic acids in biofilm structures (Karygianni et al., 2020) and organic P degradation. Furthermore, approximately 4.5% of genes (KEGG annotation) in the metagenomes were addressed to peptide degradation. Proteolysis is in fact the fastest hydrolytic activity within surface brackish and coastal sediments, aiming to provide prokaryotes with both organic C and N via amino acids, as highlighted by direct measurements (e.g., Cibic et al., 2012; Franzo et al., 2019). Together with carbohydrates, proteins and lipids are indeed the major constituents of the labile organic matter in the sediments (Pusceddu et al., 1999) and microbes are adapted to utilize these compounds heterotrophically (Oni et al., 2015), thus supporting the concept of sediments as hot spots for organic matter degradation (Lipka et al., 2018).

Within this framework of reference, we have mined our data set to test hypotheses on the putative functional microscale behavior of the microbes in the lagoon surface sediments (see Figure 9). Starting from DNA, we have asked questions about whether the community was able to display adaptive strategies such as competence and gene exchange to improve their persistence in the environment. Among the mobile genetic elements, what were the most abundant, and was it connected with the site features?

DNA is an integral part of the sediment, carrying genetic information and providing 3D structure by binding to other molecules (Dell’Anno and Corinaldesi, 2004; Corinaldesi et al., 2008; Torti et al., 2015). It is known that viruses are abundant in the sediment and are major structuring agents, keeping the microbial community under check (Danovaro et al., 2001; Breitbart et al., 2004; Breitbart and Rohwer, 2005). Among the microbial defense systems (Makarova et al., 2013), we detected the CRISPR-Cas (Doudna and Charpentier, 2014) associated genes at all the stations, despite the sediment grain size. This suggests that there is a persistent strategy of Bacteria and Archaea to fight off the viral invasions using this adaptive immunity against foreign genetic elements while creating a viral data bank for further viral encounters. Recent studies showed that in Global Ocean Sampling Expedition (GOS) and Tara Ocean data set, CRISPR-Cas were highly diverse and specific to the viral architecture (Cai et al., 2013; Nasko et al., 2019). In addition to CRISPR-Cas, Competence and Sporulation related functions were found throughout the Lagoon. DNA can be taken up during the competence window (Lennon, 2007; Dubnau and Blokesch, 2019), while a population of microbes is undergoing major environmental stress such as oxygen depletion/enrichment or nutrient scarcity. Among that population, a few cells can promote the lysis of the other individuals exploiting the released DNA for ensuring success in the future. Temporally subsequent to competence is the sporulation, an important unique microbial mechanism in marine sediment for surviving burial (Wörmer et al., 2019; Morono et al., 2020).

An ability that microbes use to move information relies on the highly diverse architecture of the mobile genetic elements, the Mobilome (Gillings, 2013; Carr et al., 2021). Prophages and Transposons were lower in Chioggia, possibly due to the site features such as higher water flow and lower TOC that might select against these mobilome strategies. Overall, mapping the sediment mobilome could give insights on rapid adaptation processes to the presence of pollutants in marine habitat due to anthropogenic activities.

Within the Regulation and cell signaling and Nanomachines categories, we have mined and detected genes related to physiology, adaptations and interactions of the microbial communities and consortia in the sediment. Biofilm in the marine environment is a pervasive structuring feature (de Carvalho, 2018). Biofilm formation is essential in keeping the microenvironment stable for the establishing of the microbial communities thus providing protection from grazers and phages, toxic substances (e.g., antibiotics), and accumulating nutrients (Flemming et al., 2016). Toxin-antitoxin (TA) systems regulate dormancy and persistence in microbes (Spoering and Lewis, 2001; Keren et al., 2004; Shah et al., 2006). TA genes are abundant on plasmids, phages, and chromosomes. They are formed by a toxin, that arrests growth by interfering with a vital cellular process, and a cognate antitoxin, that neutralizes the toxin activity to resume normal growth (Page and Peti, 2016). TA plays a major role in supporting microbial resistance to antibiotics by developing a subpopulation of persistent cells. It is interesting to note that in Chioggia, biofilm and TA were lower than the other sites, suggesting that these strategies were not favored in a sandy sediment habitat and or were not so successful as in the other sites. Despite it is known that in sandy sediments, microbes colonize grains via biofilms (Probandt et al., 2017), the biofilm might be more widespread in muddy sediment to create an organic bridge among clay and silt particles.

In terms of competition for space, nutrient, energy sources and ability to resist osmotic stress, phenazine (Price-Whelan et al., 2007; Schiessl et al., 2019), siderophores (Sandy and Butler, 2009; D’Onofrio et al., 2010), two-component systems (Held et al., 2019), osmotic stress genes, and nanomachines for cell-cell warfare were identified in the data set. Microbes constantly under siege could express phenazine, a redox-active pigment produced by Pseudomonas aeruginosa that affect gene expression, metabolic flux, and redox balancing to promote their success. Pseudomonas despite being a multi-drug resistant, biofilm former, quorum sensing human pathogenic microbe is highly versatile and has been found in the marine environment (Kimata et al., 2004). Siderophores are structures produced by microbes that chelate ions, such as Fe to facilitate the acquisition of the insoluble species thus harshening the competition for this element. We identified, within the two-component signal transduction systems (Yamamoto et al., 2005) single transduction for anaerobic respiration, nitrate and nitrite utilization, nitrogen assimilation metals efflux system, fumarate respiration, outer membrane proteins, phosphate assimilation, magnesium transport, and osmotic and envelop stress response, chemotaxis modulation, and flagellar motor switch. This great diversity found in the sediment microbial communities highlighted the adaptive strategies that microbes need to exploit while sensing the external environment at the microscale in order to activate genes that promote their success in the microniche. In the water column, a high abundance of these regulative network systems was found to be linked to copiotrophy and diazotrophy (Held et al., 2019).

The last microscale features we have identified in our dataset were the Type IV (Craig et al., 2019) and Type VI (Cianfanelli et al., 2016) secretion nanomachines, that are involved in cell-cell interactions by allowing nucleic information being shuttle from one individual to the other and attaching with effectors the other organisms. Type IV was found in microbes from acetate-stimulated aquifer sediments (Kantor et al., 2013), thus suggesting conjugation being potentially important to exchange genes and associated metabolic potentials. In the ocean, Type VI secretion systems were found to be globally distributed and associated with particle-attached microbes (Kempnich and Sison-Mangus, 2020).

Finally, motility and chemotaxis patterns, highly abundant in Chioggia, suggested a specific strategy in sandy sediments, that are more permeable to water in comparison to clayey and silty one and where a more diverse microenvironment can develop within the grains.

In sum, microscale-related genes highlight a very elaborated fine structure of interactions and behavior in the surface sediment microbial communities that then results in the large-scale biogeochemical processes on the lagoon system.

DNA metabarcoding and metagenomic analyses allowed the exploration of microbial features related to different anthropogenic pressures, of particular importance in the monitoring of fecal and sewage bacteria, antibiotic and metal resistant, hydrocarbon degraders, and xenobiotic degraders microbes.

The use of DNA-based approach is increasingly applied to microbial profiling in wastewater plants, coastal areas, and in tracking and identifying microbial signatures in the environment; such studies contributed to the identification of bacterial taxa associated with sewage- and fecal- contamination (Newton et al., 2013; Tan et al., 2015). Marine sediments are a reservoir of fecal bacteria (Luna et al., 2010) and have the potential to favor the persistence of fecal microbes and to contaminate the overlying water through resuspension, which likely poses important public health and environmental threats (Luna et al., 2016). In our study, microbial signatures showed the highest proportion of sewage-associated taxa in Marghera and the lowest proportion in Palude della Rosa, considered less impacted by human activities. Sacca Sessola and Chioggia had a higher percentage of fecal-associated taxa, being mostly affected by urban contamination (Picone et al., 2016). Overall, we detected several indicators of fecal and sewage contamination (see Figure 6), and their relative abundance and taxonomic assignment are in line with a previous study in central part of the Venice Lagoon (Luna et al., 2016), confirming a diffuse contamination and accumulation of these bacteria in the Lagoon sediment.

Antibiotic resistance spread is one of the main threats to global human, fish, shellfish, and in overall ecosystem health (Yang et al., 2012; Ventola, 2015; Mackenzie and Jeggo, 2019; Helsens et al., 2020). For this reason, the detection, tracking, and assessment of the increasing antibiotic resistant microbes and genes, is necessary in terms of risk management (Heß et al., 2018). The marine environment is considered a global reservoir of antibiotic resistance genes, but we still have little understanding of their presence and diversity in both pelagic and benthic environments (Chen et al., 2013; Hatosy and Martiny, 2015). Several antibiotic resistance genes were detected in our dataset (see Figure 8), with a significantly higher presence in the two more impacted sites of Marghera and Tresse. The Multi-drugs and Beta-lactamase resistance genes were found as the most abundant at all sites, as previously detected in seawater (Hatosy and Martiny, 2015) and sediment (Huang et al., 2019).

Surface sediment metal contamination is another consequence of anthropogenic pressure in transitional environments, representing a potential threat to both public health and marine ecosystems (Han et al., 2011). Metals not only accumulate in the sediment but can be subjected to complex cycling that can make toxic compounds available to benthic organisms; sediment mobilization and resuspension under particular pH and redox conditions may release them in the water column becoming a secondary source of pollution (Turner and Millward, 2002; Eggleton and Thomas, 2004; Zonta et al., 2018). Furthermore, microbes that are antibiotic resistant can carry information for metal resistance and vice versa, thus exacerbating the problem of resistance spreading (Chen et al., 2019). Several toxic compound resistance genes were detected in our dataset, with a significantly higher presence in the two more impacted sites of Marghera and Tresse. These sites presented the highest number of genes related to Cobalt-zinc-cadmium resistance, and to Mercury resistance, due to the release of such metals from industrial plants. In particular, the mercury contamination mainly originates from the chloralkali discharge in Porto Marghera, which occurred until the 1980s (Bloom et al., 2004). Mercury is one of the most severe pollutants, and its flux to Venice Lagoon water is mostly due to the resuspension of contaminated sediment (Bloom et al., 2004), underlying the importance of its regulated management. Chioggia and Palude della Rosa presented the highest number of genes coding for Arsenic resistance. Even if the presence of this contaminant in the lagoon is mainly due to the geochemical characteristics of the sediments, human activities can contribute to its accumulation (Zonta et al., 2018).

Chioggia presented the highest number of genes related to hydrocarbon degradation. Interestingly, in this site, both surface and sub-surface sediment were recently found to harbor the highest concentrations of PAH and PBC in the whole Venice Lagoon (Cassin et al., 2018), due to the important urban waste discharge, road/boat traffic, and being located close to the harbor of the second Italian largest fishing fleet.

It was worth noting that among the identified hydrocarbon degraders some taxa were also able to utilize lignin. Lignin is a very complex and hard-to-degrade polymer (Hedges and Mann, 1979) and Bacteria, Archaea and Eukaryotes have evolved adaptive strategies to gain energy from it (Cragg et al., 2015). In the Venice Lagoon, wood has been historically and extensively used for maritime needs, like mooring piles and navigation channel delimitation (named “briccole,” approximately 7,000, formed by two or three poles) (Ghirardini et al., 2010). Briccole are subjected to physical and biotic degradation, the latter mainly due to marine borers (Ghirardini et al., 2010). Since our data highlighted that the sediment represents a reservoir of potential lignin-degrader bacteria, we hypothesize a role of these taxa (alone or in consortia with other microorganisms) in the turnover/recycling of human-introduced wood structures in the Lagoon.