Difference in composition and functional analysis of bacterial communities between Mytilus galloprovincialis gills and surrounding water in a brackish inshore bay, analyzed by 16S rDNA multi-amplicon sequencing

Abstract

Background

Lagoons are important natural systems, with attractive favorable conditions for aquaculture production, such as shellfish cultivation. Despite their socio-economic relevance for human activity, information on the microbial diversity, community composition, and putative functions of gill-associated microbiota and seawater is still limited, particularly as regards the extent of specific taxa enrichment within the gills and the relative effects of the temporal and spatial variables. In this study, we used a 16S rDNA multi-amplicon sequencing approach using an Ion GeneStudio S5 System and a function prediction method (Functional Annotation of Prokaryotic Taxa (FAPROTAX), to inspect the springtime dynamics of microbial communities and their inferred metabolic features in an Adriatic lagoon (Po Delta, Italy).

Results

Mussels and surrounding seawater were sampled in two rearing areas three times between April and June 2021. Sequencing results showed significant (p ≤ 0.05) differences in bacterial community composition and diversity between gills and seawater. Gills were dominated by the Methylobacterium-Methylorubrum and Burkholderia-Caballeronia-Paraburkholderia genera, while in seawater samples Izamaplasma, Planktomarina, and Candidatus Aquiluna were detected as being dominant. The microbiota composition did not differ significantly between the two rearing areas. The sampling time, although limited to a 3-month timeframe, instead revealed a structural variation of the bacterial profile both in gills and seawater for alpha and beta diversities respectively. The functional prediction analysis highlighted an overexpression of human gut-associated bacteria in relation to the season-related increase in seawater temperature.

Conclusions

These findings enhance our understanding of the differences between gill-associated and seawater microbiota composition and provide novel insights into the functions carried out by bacteria inhabiting these niches, as well as on the key host-symbiont relationships of bivalves in lagoon environments.

Background

Mussels (Mytilus galloprovincialis; Lamarck, 1819) are animal feeders able to filter large volumes of water through their gills to satisfy their food requirements (Wathsala et al. 2021). During this powerful filter-feeding activity, organic particles such as algae, bacteria, viruses, inorganic contaminants, and detritus can be retained in animal tissues (Orban et al. 2002; Vernocchi et al. 2007). Indeed, mussels have long been utilized in the biomonitoring of environmental quality (Gagné et al. 2019). Preliminary studies suggest that despite the environment influencing mollusk bacterial composition, the process of selective filtering determines a tissue-specific microbial profile (Franzellitti et al. 2016; Freitas et al. 2019; Musella et al. 2020; Pathirana et al. 2019; Schoinas et al. 2023; Trabal Fernandez et al. 2014). These bacterial communities are assembled progressively from cells ingested through gills and can potentially act in metabolite synthesis and/or excretion, osmoregulation, and invasion success, playing an important role in the processing of nutrients and defense against pathogens (Jahnke et al. 1995). However, despite its importance, knowledge of the mussels’ microbiota is still limited.

In particular, the microbiota of bivalve gills is considered autochthonous, due to the higher presence of symbiotic bacteria than in other organs (Balbi et al. 2020; Cappello et al. 2015).

Examples of symbioses include Campylobacteria and Epsilonbacteria in Bathymodiolus platifrons bivalves (Lin et al. 2023), sulfur-oxidizing bacteria of the genus Beggiatoa (Imhof et al. 2003), and the Methylomonaceae group (Lin et al. 2021). These symbionts can provide mussels with probiotic functions, with final consequences on animals’ health and productivity (Ikeda-Ohtsubo et al. 2018; Simon et al. 2019). Consequently, the gill microbiota should be considered an important target for understanding mussel physiology and management (Lokmer et al. 2016; Roterman et al. 2015; Zurel et al. 2011).

In Italy, the Po River Delta system is one of the most productive areas of the Mediterranean with a famous aquaculture activity (Zoppin et al. 2020). Among farming areas, the Scardovari lagoon presents environmentally favorable conditions for the development of thriving mussel farming (Caroppo et al. 2018; Cibic et al. 2019). The area (32 km²) is characterized by low water levels with a constant change linked to the cycles of tides and a significant supply of fresh water and nutrients. These delta areas are subjected to different stress factors of natural and anthropogenic origin that can influence and modify the microbial community composition of the surrounding environment and mussels, with possibly severe consequences for the lagoon ecosystem and seafood production (Franzo et al. 2019; Mascolo et al. 2019).

For this reason, the study of bacteria profiles of mussels and seawater in relevant geographic locations could cast some light on their status and dynamics.

The aim of this project was to explore the microbiota of M. galloprovincialis at the gills’ level and in comparison to that of their surrounding water in samples collected at Scardovari lagoon, and their temporal succession, to provide basic knowledge on the respective microbial community structures, and the possible presence of specific bacteria hosted in the animal tissues. Moreover, by characterizing the bacterial content in two different areas of the lagoon (north-west and south-east) between April and June 2021, a further aim was to better understand the successional prokaryotic abundance in gill tissues with respect to the seawater reservoir source also as a function of the seasonal environmental changes of physical and chemical parameters. To explore as deeply possible the complexity of microbial communities, we used the Ion GeneStudio S5 technology to simultaneously amplify seven hypervariable regions of the bacterial 16S rRNA gene. Indeed, all the regions of this target gene are considered significantly necessary to obtain an accurate estimation of taxonomic diversity (Claesson et al. 2010; Liu et al. 2008). The working hypotheses under investigation, yet unexplored by prior studies, were the following: (1) the selective fractionation that leads to a community structure in gills which differs from the surrounding water source varies during the season and (2) the temporal effects are hierarchically more effective than spatial ones in shaping community composition of either water or gills. As innovative elements of the present approach, in the first instance, no studies have dealt with the time-related variation of mussel gills microbiota through a consecutive monthly timeframe during the transition from spring to summer, and second, our report addressed the difference between the animal tissue and its surrounding habitat to assess overlaps, enrichment, and divergences in the comparison of the two connected microbiomes of seawater and gills.

Results

Physical-chemical seawater parameters

Spatial and temporal dynamics of the physical-chemical parameters are shown in Fig. 1. Statistical differences between the northwest and southeast of the lagoon were detected for temperature and salinity at all sampling times (Supplementary Table S1). The chlorophyll was significantly higher (p ≤ 0.05) in the southeast compared to the northwest site in April and May. As regards the pH, statistical differences between the two areas were found during April and June (Supplementary Table S1). No differences were found for the dissolved oxygen at any of the sampling times.

Fig. 1

Dynamics of the monitored physical-chemical parameters. Bars indicate standard errors

Looking at the April–June temporal dynamics, in the northeast area the temperature, salinity, dissolved oxygen, and chlorophyll showed statistically significant differences across the three sampling times (p ≤ 0.05) (Supplementary Table S2).

The highest temperature and chlorophyll values were detected in June, with a mean value of 28.2 °C and 7.4 μg/l, respectively. Salinity reached its maximum value in May, with a mean of 28.7‰. Different dissolved oxygen concentrations were also found during the sampling period, with the highest value in April. No differences were found for the pH parameter in relation to time.

In the southwest area, the temperature and salinity were higher in June with respect to April and May (Table S2). pH reached its minimum in June, with a mean of 8.3. Statistical differences were also detected for dissolved oxygen and chlorophyll, with the highest values during April.

To further explore the relative role of the physical-chemical parameters in shaping the bacterial communities redundancy analysis (RDA) was used, the results of which are shown in Supplementary Figure S2. Redundancy analysis (Legendre & Legendre 1998) takes as input a site/data matrix where each site has given values for one or more environmental/explanatory variables as well as a number of response (dependent) variables. The strongest role was displayed by temperature, followed by salinity, while pH had a minor incidence. The positions of dissolved oxygen and chlorophyll resulted in between.

Microbiota of Mytilus galloprovincialis and the surrounding water

A total of 5,219,296 reads were obtained from the sequencing. The mean number of reads per sample was 175,872 + 116,536. As previously mentioned, any samples which failed to satisfy customarily sufficient reads abundance were removed. Specifically, two-gill samples (one from April in the northwest area and one from April in the southeast area) and one seawater sample (from May in the northwest area) were excluded from the dataset. High-quality reads yielded 19,078 amplified sequence variants (ASVs).

The overall microbiota profile of M. galloprovincialis gills is represented in Fig. 2. The phylum Proteobacteria (88%) was the most abundant. The subdominant phyla were Cyanobacteria (8%) and Actinobacteria (2%). Further phyla were represented by Firmicutes and Bacteroidetes (both with an abundance of ~1%). In particular, Proteobacteria were represented especially by the classes of Alphaproteobacteria (60%) and Gammaproteobacteria (28%). The families of Sphingomonadaceae (31%) and Beijerinckiaceae (28%) were mostly represented in Alphaproteobacteria, and Sphingomonas and Methylobacterium-Methylorubrum were the dominant genera. The class of Gammaproteobacteria was represented by the family Burkholderiaceae (26%), and specifically, the Burkholderia-Caballeronia-Paraburkholderia genera complex was highly featured. Other minor genera from different families (<1%), including Dyella, Nocardioides, Brevundimonas, Pseudomonas, and Staphylococcus, were detected.

Fig. 2

Bar plots summarizing the microbiota composition of the 10 most abundant bacteria taxa for gill and seawater samples. For gills and seawater, the samples with different sampling times and sites are averaged together within the same histogram bar. Other taxa are grouped into “Other”

As for seawater, Proteobacteria (48%) was the dominant phylum (Fig. 2). The subdominant phyla were Bacteroidetes (17%), Actinobacteria (10%), Firmicutes (9%), Cyanobacteria (7%), Desulfobacteriota (2%), and Verrucomicrobiota (2%). Proteobacteria were represented especially by the classes of Alphaproteobacteria (30%) and Gammaproteobacteria (30%). The most represented families within these classes were Rhodobacteraceae (18%), Flavobacteriaceae (9%), and Microbacteriaceae (7%). Halieaceae (2%), Sphingomonadaceae (2%), and Vibrionaceae (1%) were minor components. At the genus level, Izamaplasma (9%), Planktomarina (8%), and Candidatus Aquiluna (6%) were the dominant taxa. The most represented subdominant genera were Rhodococcus and Litoricola (2%), followed by other genera with an abundance of 1%, such as Vibrio.

Beta diversity between gills and seawater microbiota was evaluated using a PCoA analysis. The samples’ separation into distinct clusters was clearly evident (PERMANOVA p ≤ 0.05) (Fig. 3A). Species richness and abundance were higher in seawater than in mussels (p ≤ 0.001) (Fig. 3B). The differential abundance analysis revealed the presence of 509 genera differentially expressed between the two groups (p ≤ 0.05), in particular, Rhodoplanes, Streptomyces, Oerskovia, Nocardia, Blastococcus, Comamonas, Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium, Pedobacter, Ralstonia, Planococcus, Methylobacterium-Methylorubrum, Arthrobacter, and Bulkholderia-Caballeronia-Paraburkholderia were up-represented in gill samples compared to water.

Fig. 3

Multivariate analysis and species diversity plots for the bacterial communities in gill tissues and seawater samples. Principal Component Analysis (PCoA) (A). Box plot of alpha diversity using the Shannon index (B). The significance level (Mann-Whitney test) is reported above the B graph

Comparison of the bacterial communities between the northwest and southeast lagoon sites

To investigate the differences in microbiota structure in the two different rearing sites, a multivariate analysis of gills and seawater microbial compositions was conducted. The PCoA analysis showed that the bacterial profile of the two sites does not cluster separately. Indeed, the PERMANOVA test provided no significant p value for either gill or seawater samples. As for the alpha diversity, a comparison between the two sites (northwest and southeast) did not show significant differences in species richness and abundance (Figure S1).

Impact of sampling time on microbiota profile

Aiming to evaluate the sampling time effect on the microbiota profile of seawater and gills, we compared the bacterial profiles of samples collected at three distinct times, i.e., April, May, and June 2021. In gill samples, the PCoA analysis, based on Bray-Curtis distances, did not evidence a clear pattern of cluster separation (Fig. 4A). Indeed, the PERMANOVA pairwise test provided no significant p value between groups. Alpha diversity comparison of groups revealed a significantly lower value in community diversity in April compared to May and June (p ≤ 0.05) (Fig. 4B). The differential abundance analysis showed 9 genera differentially expressed between June and May, 3 genera between May and April, and only one genus with a statistically different abundance between June and April (Fig. 5A).

Fig. 4

Principal coordinate analysis of gill samples between different sampling times (A). Box plot of alpha diversity of gill samples with Shannon index (B). Principal coordinate analysis of seawater between different sampling times (C). Box plot of alpha diversity of seawater with Shannon index (D). The significance level (Kruskall-Wallis test) is reported above the B and D graphs

Fig. 5

Clustered bar chart with the Log₂ fold change values of different (Wald test, p ≤ 0.05) genera between the three sampling times in gill samples (A) and seawater (B). Significant differences are highlighted with the (*) symbol

The seawater samples were clustered into three different groups (PERMANOVA p ≤ 0.002) (Fig. 4C). No differences in bacterial abundance and richness were detected during sampling times by the Shannon index (Fig. 4D). At the genus level, the differential abundance analysis highlighted 15 bacterial genera with significantly different abundances between May and June, 16 between June and April, and 6 between May and April (p ≤ 0.05) (Fig. 5B).

To further explore the effect of this as well as those from the other variables in relation to all aspects of beta diversity in terms of taxa distribution and their respective prevalence, we conducted a similarity percent (SIMPER) analysis, the results of which are shown in Supplementary Table S3. Inspecting the scores to assess which taxa are primarily responsible for the observed difference between samples, it can be seen that both for the gill vs. water comparison and for that between the northeast vs. southwest shore sites, the top of the score was occupied by Burkholderia-Caballeronia Paraburkholderia, while in the contrasts between sampling months, from May onwards, the leading taxon resulted as Methylobacterium-Methylorubrum. The beta-diversity was further inspected by performing a Morisita-Horn similarity measurement in non-metric multidimensional scaling and focusing on four contrasts: the difference between South and North shore site for gills, between South and North shore site for seawater; among months for gills, and among months for seawater. Results are shown in Supplementary information Fig. S3. The months-by-seawater contrast resulted the least overlapping for the involved sets.

Functional prediction of the microbial communities of M. galloprovincialis and seawater

Functional phenotype prediction based on the taxonomical information with FAPROTAX allowed us to infer different predicted pathways between gills and seawater. Overall, 40 functional groups were identified in gills; the dominant ones pointed to the category of human pathogens and associated, oxidation of sulfur compounds, animal parasites, and symbionts. In seawater, the dominant functional groups included methanotrophy, chlorate reducers, fish parasites, and the respiration of sulfur compounds. Twenty-one functional groups showed statistically significant differences between gill samples and seawater bacterial communities (Wald test, p ≤ 0.05) (Fig. 6). The functional groups that were more abundant in gill samples consisted of the oxidation of sulfur compounds, animal parasites, symbionts, and intracellular parasites.

Fig. 6

Clustered bar chart with the Log₂ fold change values of significantly different (Wald test, p ≤ 0.05) functional groups between gill samples and seawater. Genera with positive values are overrepresented in seawater. Genera with negative values are overrepresented in gill samples

Different relative abundances of functional groups in gill samples were evidenced between June and May. The gills microbiota in June showed enrichment in the oxidation of sulfur compounds, animal parasites, symbionts, intracellular parasites, and non-photosynthetic cyanobacteria. The Wald test also revealed the presence of differences (p ≤ 0.05) in functional bacteria groups between sampling times in seawater samples. Eighteen differentially functional groups were evidenced between April and May, 12 between June and April, and 14 between May and June (Fig. 7).

Fig. 7

Clustered bar chart with the Log₂ fold change values of different (Wald test, p ≤ 0.05) functional groups between the three sampling times in seawater. Significant differences are highlighted with the (*) symbol

Discussion

In this study, we aimed to characterize and compare the composition of gills’ bacterial communities of M. galloprovincialis and those of the surrounding water in a brackish inshore bay (Scardivari lagoon), with a multi-amplicon sequences technology. In agreement with other studies on mussels’ gills microbiota (Brito et al. 2018; Lin et al. 2021; Musella et al. 2020), sequencing results revealed the dominant presence of Proteobacteria and Actinobacteria at the phylum level. Within the group of Proteobacteria, our results highlighted a high relative abundance of the classes of Alphaproteobacteria and Gammaproteobacteria, which are well-known highly abundant marine microorganisms, as reported in previous studies on the 16S rRNA gene (Fernandez-Piquer et al. 2012; Zurel et al. 2011). As for seawater, the dominant phyla were those of Proteobacteria and Bacteroidetes, followed by the groups of Actinobacteria and Firmicutes. These findings are consistent with other studies on microbial communities of coastal lagoons (Musella et al. 2020; Trombetta et al. 2022), which highlighted a relative abundance of Proteobacteria from 50 to 70%, and Bacteriodetes 15–30%; with a minor presence of Actinobacteria. Interestingly, the Firmicutes, which usually have an abundance of less than 1% (Musella et al. 2020; Trombetta et al. 2022), reached 9% in the environment that we analyzed. In this respect, the phylum being typically endowed with spore-forming capabilities, the differential yield across reports can be linked to different efficiencies in the DNA extraction particularly in the mechanical and chemical lysing steps. In our case, the 10-min tissue lyser treatment may have warranted a thorough extractive power. According to our alpha and beta diversity results, the gills and seawater microbial profiles presented a markedly different structure. Consistently, Musella et al. (2020) showed the Mytilus galloprovincialis microbiota well differentiated from that of seawater. As an advancement of prior knowledge, in our study, we also show that the microbial gill composition is affected by seasonal dynamics and it changes upon sampling times distributed in a continuous series within the same season. We found gills microbiota was dominated by Methylobacterium-Methylorubrum, Burkholderia-Caballeronia-Paraburkholderia, and Sphingomonas. The differential abundance analysis revealed the first two genera were over-represented in gill tissues. The SIMPER analysis (Supplementary Table S3) allowed us to see the respective occurrence of these dominant taxa and their shifts of prevalence in dissimilarity when comparing temporal or spatial aspects, as well as the gills versus the surrounding water habitat. Methylobacterium-Methylorubrum is a facultative methylotrophic symbiont that can utilize one-carbon atom compounds during its metabolism as an energy source. Species of this genus have often been found associated with plants, mussels, and soil (Jia et al. 2020; Jiang et al. 2020). The abundance of these bacteria in mussels has been reported in different environments and their presence has been recognized to contribute to the optimization of energy yield (Bergquist et al. 2004; Duperron et al. 2007; Lin et al. 2021). Burkholderia is a ubiquitous bacterium that can produce secondary metabolites with antagonistic effects on pathogenic bacteria (Trabal Fernandez et al. 2014). This bacterial genus can also act as a symbiont due to its capacity to use methylotrophic, nitrogen, and hydrocarbon compounds during energy metabolism (Mills et al. 2012; Pérez-Pantoja et al. 2011).

The functional prediction of microbial community evidenced the dominant presence of sulfur-oxidizing bacteria, animal parasites and symbionts, and intracellular parasites in gills tissues.

Aerobic sulfur-oxidizing prokaryotes belonging to different genera such as Methylobacterium and Bathymodiolus are known for symbiotic associations that sustain host metabolism (Friedrich et al. 2001). Interestingly, FAPROTAX also evidenced functional groups related to human pathogens. Among the dominant taxa, the Burkholderia and Sphingomonas genera can contain possible opportunistic pathogens that play a role in human and animal diseases (Espinosa-Victoria et al. 2020). Nevertheless, it needs to be remembered that this and any other tools basing functional inference on taxonomical data yield solely predictive indications. Timewise, as shown in Supplementary Table S3, Sphingomonas showed an over 30-fold increase from April to May followed by a 3-fold reduction from May to June. Moreover, the Wald test also detected an increase in Ralstonia and Comamonas in gills tissue with respect to seawater, which can contain some possible pathogenic bacteria. However, some members of these genera are also capable of nitrogen reduction and methylotrophy symbiotic metabolism (Ryan and Adley 2010). Therefore, further studies are required to better understand the relations of these bacterial genera with gill tissues.

In our study, we also investigated the effect of sampling sites on microbiota communities. The alpha and beta diversity analysis did not show differences in species richness and abundance in either seawater or mussels as regards the two sides of the bay. Recently, Lin et al. (2021) found similar gill bacteria communities in Bathymodious mussels from the adjacent cold seep and hydrothermal vent environments. In contrast, previous studies that characterized the influence of location on the mollusks found significant differences in microbial communities (Valenzuela et al. 2021).

Data collected in this study identified several differences between the two sites, suggesting that in the farming area of the Scardovari lagoon, although successional changes occur, the overall ecology indexes for the corresponding different communities of gills and seawater maintain similar levels.

As mentioned, we specifically explored the effect of sampling time on microbial communities. Previous studies on marine organisms evidenced a significant temporal heterogeneity in microbial composition (Li et al. 2019; Zurel et al. 2011). These findings are consistent with our results, where a sampling time effect was detected in both gill samples and seawater.

With regard to gill samples, the bacterial diversity resulted significantly lower in April than in May and June, when it reached similar and apparently stable values. This pattern of alpha diversity under our conditions differed from a previous study (Li et al. 2019), where warmer temperatures led to a temporal decrease of the microbial diversity in the gut of Mytilus coruscus. However, such an effect has to be considered valid for the limited quarter of the year that has been examined.

Interestingly, the seawater Shannon index did not differ among sampling times, but the PCoA showed a specific and significantly different clustering (PERMANOVA, p < 0.05) of the three groups. During our study, sampling time shifts in environmental conditions were dominated by changing temperature, salinity, dissolved oxygen, and chlorophyll. Some reports have demonstrated that bacterial community structure is influenced by several environmental factors, such as pH, salinity, dissolved oxygen, nutrients, and especially water temperature (Gerpe et al. 2021; Lokmer et al. 2016).

The functional predictor FAPROTAX detected sampling time-related changes also in functional groups, with enrichment of bacteria endowed with the oxidative potential of sulfur compounds in gill samples in the June sampling. Previous studies demonstrated the presence of sulfur-oxidizing symbionts that could coexist in the deep-sea Bathymodiolus mussels (Ling et al. 2020). As for seawater, the functional groups were affected by the sampling time, with enrichment of thiosulfate respiration, and aromatic and hydrocarbon degradation pathways, suggesting a possible higher anthropogenic pollution by aromatic and hydrocarbon compounds. Moreover, human and mammal gut functional groups were detected in May and June, suggesting a higher anthropic impact in these 2 months. Moreover, gut-derived bacteria are well-known as potential pathogenic microorganisms. These findings are in agreement with previous reports that showed that warm temperatures can alter microbial communities’ composition, favoring the shift from beneficial bacteria to opportunistic pathogens (Erwin et al. 2012; Zhang et al. 2021).

Conclusions

Our study provides innovative information about the composition and putative function of microbiota present in the gills of Mytilus galloprovincialis and the surrounding water in the lagoon ecosystem, to advance our understanding of their dynamics and composition. Our results highlighted a characteristic gills’ bacteria profile, with a sharply different structure from the seawater micro-ecosystem. Moreover, the identification of three main bacterial groups in gills’ tissues, extends the knowledge of the symbiotic association between bacteria and mussels.

The study also provides novel information about the environmental influence on the microbiota of mussels and seawater in mussel-farming natural lagoons. The farming areas are not the primary factor that affects the microbiota profile of gills and seawater. On the contrary, sampling time differences were reflected by bacterial community structure with an increase across the sampling times of potential opportunistic bacteria.

In relation to our two tested hypotheses, we found that in relation to either alpha or beta diversities: (1) not only the animal-mediated selective enrichment leads to a bacterial community structure in gills different from that of the surrounding water source but also that such mussel-specific assemblage varies during the season both in terms of taxonomy and as predictable phenotypic traits and (2) the hierarchical ranking of the variables driving the change within the chosen basin, either in water or within gills, sees the temporal dimension dominate over the spatial one.

These findings suggest that the knowledge of autochthonous bacteria communities associated with mussels and seawater in the Po River Delta areas is essential to predict their changes and dynamics, to better manage mollusk production, and to preserve the lagoon environment. Future studies should focus on the seasonal dynamics of the mussel microbiota in different organs, to better understand bacterial variability at the tissue scale and also provide more information about the environmental effects on the Mytilus galloprovincialis microbiota.

Materials and methods

Study area and physical-chemical seawater quality parameters

Mussels and seawater were collected in two areas of the Scardovari lagoon (position: north-west 44° 51′ N, 12° 43′ E and southeast 44° 51′ N, 12° 40′ E) (Fig. 8) for three times between April and June 2021 on the following dates: April 21, 2021, May 18, 2021, and June 16, 2021; inside a mussel farm located in the Po Delta, Porto Tolle (RO) (Italy). The farm is situated between two terminal branches (Tolle and Gnocca) of the Po River. The collection area extends over 535 ha and has an average depth of 3 m. Samples were taken along the mussel ropes within a 1-m-long stretch spanning between 2 and 3 m of depth from the water surface.

Fig. 8

Geographical localization of the two rearing areas (black stars; N: northwest 44° 51′ N, 12° 43′ E; S: southeast 44° 51′ N, 12° 40′ E) in the Scardovari lagoon

The data related to temperature, salinity, dissolved oxygen, pH, and chlorophyll of water were obtained from multiparameter ARPAV probes measured every 30 min/day (ARPAV 2021); for every parameter, the mean daily values recorded were analyzed on both the north and south side at every sampling time.

Sampling and tissue preparation

At every harvest, 8 kg of mussels were collected, and 96 individuals were randomly picked to be immediately transported to the laboratory under refrigeration in coolers (4°C). After transport, the mussels were brushed, washed under water, and cleaned with ethanol (96%); a sterile scalpel was used to pry them open, which allowed draining off by gravity and discarding the intra-valvular fluid, after which the flesh was dissected under microbiologically controlled conditions. For each animal, 200 mg was taken from the gill tissue, placed in a 2 mL Eppendorf tube, and stored at −80 °C.

One liter of seawater was collected at a depth of 3 m using a sterile tube with a diameter of 10 cm and a length of 0.5 m that was made to rise along the water column. From that amount, 50 mL of seawater was taken and placed separately in two Falcon tubes (total 100 mL), transported in coolers (4 °C) to the laboratory, and then stored at −80 °C until analysis.

DNA extraction

For mussels, the bacterial DNA extraction from gills tissues was performed with the following protocol: 400 μL of RTL buffer (guanidine thiocyanate 0.12 M) (Qiagen, Hilden, Germany) was added to 200 mg of tissue in a 2 mL Eppendorf tube, together with one metallic bead. The homogenization step was performed with the use of a Tissue Lyser (Qiagen, Hilden, Germany) for 10 min at 30 Hertz. After the supplement of 400 μL of RTL, the samples were subjected to enzymatic treatment with lysozyme (50 mg/mL, Thermo Fisher Scientific, Waltham, MA, USA) and proteinase K (20 mg/mL, Thermo Fisher Scientific, Waltham, MA, USA). At first, 20 μL of lysozyme was added; samples were then incubated for 1h at 37 °C. After this step, tissues were treated with 20 μL of proteinase K and incubated for 1h at 50 °C. At the end of the enzymatic treatment, samples were centrifuged for 6 min at 2100 × g, after which the supernatant was collected.

Seawater samples were centrifuged in the Falcon tubes at 5500 rpm for 10 min. The supernatant was removed, and 300 μL RLT buffer was added to resuspend the pellet. The samples were then transferred to a 2 mL Eppendorf tube and subjected to enzymatic treatment with Proteinase K at 50 °C for 1h. In the end, samples were again centrifuged for 5 min at 6000 × g, after which the supernatant was collected.

DNA purification involved, for all samples, the use of a Biosprint 96 automated workstation extractor (Qiagen, Hilden, Germany): 300 μL of lysate was transferred into an S-block with 200 μL isopropanol and 20 μL MagAttract magnetic beads suspension (Qiagen, Hilden, Germany) and loaded into the instrument. The protocol involved the use of five other plates: one S-block with 500 μL of RPW buffer (guanidine hydrochloride, 1.31M; Qiagen, Hilden, Germany), two plates with 500 μL of ethanol (96%), another S-block with 500 μl of tween solution at 0.02% and a last flat 96-well plate with 100 μL of nuclease-free water where the DNA was eluted. The final concentration of DNA was measured with a Qubit 3.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) using the Qubit™ DNA HS Assay Kit Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA).

Gills and seawater bacterial 16S metabarcoding

For mussels, 5 μL of the extracted DNA solution was taken from every sample and used to create three pools of thirty-two samples for each site and time for a total of 18 pools. Twelve seawater samples were analyzed.

The library preparation started with the amplification of hypervariable regions (V2, V4, V8 and V3, V6-7, V9). To perform this first step, we utilized the 16S Ion Metagenomics Kit (Thermo Fisher Scientific, Waltham, MA, USA) and the following amplification program: 95 °C for 10 min, 25 cycles of 95 °C for 30 s, 58 °C for 30 s, 72 °C for 20 s, and a hold stage at 72 °C for 7 min. After PCR a cleaning step was performed, and libraries were normalized to 30 ng/μl. For barcode ligation, the Ion Xpress Plus 9 Fragment Library Kit (Thermo Fisher Scientific, Waltham, MA, USA) and Ion Express Barcode Kit (Thermo Fisher Scientific, Waltham, MA, USA) were used, and each sample was ligated with a distinctive barcode. Another amplification step was performed with the following program: 95 °C for 5 min, 7 cycles of 95 °C for 15 s, 58 °C for 15 s, and 70 °C for 1 min, then stored at 4 °C. Amplified libraries were quantified using a Qubit 3.0 Fluorometer with Qubit™ DNA HS Assay Kit to pool libraries at a final concentration of 100 pM. The pool was then processed by the Ion 520^TM & Ion 530 ^TM Kit – OT2 400bp (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s instructions. In the end, the sample was loaded on the Ion 520 chip (Thermo Fisher Scientific, Waltham, MA, USA) and the sequencing run was performed with Ion™ GeneStudio S5 System (Thermo Fisher Scientific, Waltham, MA, USA). Sequencing data have been deposited in the NCBI Sequence Reads Archive SRA under project Code PRJNA954590.

Bioinformatics and statistics

Raw reads were processed to trim 20 base pairs on both ends to remove primers using cutadapt v3.4 (Martin 2013) and analyzed using QIIME2 v2021.4 (Bolyen et al. 2019). The “qiime dada2 plugin” was used to denoise and dereplicate high-quality reads into amplicon variants (ASVs). An alpha rarefaction plot was produced using the “qiime alpha diversity” plugin to determine if the sequencing depth was sufficient. The SILVA SSU v138.1 database was additionally utilized as a reference for the taxonomic assignation of ASVs (Quast et al. 2013). Bacterial diversity and composition, including the Morisita-Horn beta diversity distances, were estimated with the MicrobiomeAnalyst online utility (https://www.microbiomeanalyst.ca/) (Chong et al. 2020); before the analysis, any samples not meeting customarily sufficient reads abundance were discarded and subsequently data were filtered to exclude reads occurring less than 4 times as standard safe practice to minimize biases related to sequencer error rate. Finally, a TSS (total sum scaling) data transformation was performed. Alpha diversity was evaluated using the Shannon diversity ecological index and assessing the significance of differences by the Kruskal-Wallis/Mann-Whitney statistical test (FDR-corrected p ≤ 0.05). Beta-diversity among samples was analyzed with principal coordinate analysis (PCoA) based on the Bray-Curtis distances to which the permutational MANOVA (PERMANOVA) statistical analysis (FDR-corrected p ≤ 0.05) was subsequently applied in order to investigate the statistical significance of sample division. Differential abundance analysis, using the Wald test and TMM (trimmed mean of M values) normalization, was performed with the CLC Genomics Workbench (23.0.1) (Qiagen, Hilden, Germany) dedicated tool, to highlight significant differences (Bonferroni-corrected p ≤ 0.05) within bacteria taxa. The functional prediction was performed with FAPROTAX (Functional Annotation of Prokaryotic Taxa; Louca et al. 2016) using the default parameters with TSS normalization, to predict the main ecological processes of the microbial communities associated with gills and the surrounding seawater. The input was the genus level counts file. The Wald statistical test with the CLC Genomics Workbench (23.0.1) dedicated tool was performed to highlight significant differences (Bonferroni-corrected p ≤ 0.05) within bacteria functional groups.

The physical-chemical seawater parameters were processed using the STATISTICA 7 software package (StatSoft). Comparisons among means were performed by the Wilcoxon test (Bonferroni-corrected p ≤ 0.05).

The similarity percent SIMPER and multivariete redundancy analysis RDA were calculated using the PAST Software (Hammer et al. 2001).