Eukaryotic and prokaryotic diversity was evaluated in three species of bivalves sampled from Meira (Ría de Vigo, NW Spain). The sampled bivalves were mussels (Mytilus galloprovincialis), clams (Ruditapes philippinarum) and cockles (Cerastoderma edule), with 15 individuals of each species collected at different sampling points. The first sampling was performed in the summer of 2016, and the subsequent samplings were performed approximately every three months (i.e., seven sample collections from 2016 to 2018). Mussels were sampled in the intertidal zone, and cockles and clams at 6–7 cm deep in the sediment. Figure 1A shows the experimental layout. After the collection of specimens, they were transported to the laboratory, where a fragment was fixed for further histological procedures. The other part of the samples was preserved in 99–100% ethanol at room temperature (25°C) to enable DNA extraction. For the subsequent analyses, some data from the first sampling point (July 2016) were not available due to the poor quality of the DNA. These samples were: cockle (for 18S rRNA gene) and mussel and cockle (for 16S rRNA genes). R. philippinarum was included despite being non-native to European waters as it constitutes a significant component of the local aquaculture industry, allowing for a comparison of microbiome adaptation in introduced versus native species within the same ecosystem.

The whole-body sampling approach, while enabling robust species comparisons, limits the resolution of tissue-specific microbiome patterns that may exhibit seasonal variations. Nevertheless, this methodology offers a comprehensive overview of host-specific microbial signatures while minimizing potential bias from organ-specific sampling variations.

Oblique transverse sections, approximately 5 mm thick, including the mantle, gonad, digestive gland and gills, were taken from each specimen after fixation in Davidson solution (Shaw and Battle, 1957). Tissue samples were processed in an automatic tissue processor (Leica TP 1020, Leica Microsystems, Germany), embedded in paraffin wax blocks and cut on a microtome. Sections of 5 µm were first deparaffinized, rehydrated and finally stained with hematoxylin-eosin (Howard et al., 2004). Histological sections were examined for the presence of parasites and pathological alterations using an optical microscope (Nikon Eclipse 80i) at 10X and 40X magnification.

For DNA isolation, a longitudinal section of the animal body, including fragments of the gonad, mantle, gill, and digestive gland was collected, pooled, and homogenized with a pestle to obtain a representative DNA sample. DNA extraction was performed in batches as sampling progressed using the Maxwell 16 Blood Purification DNA Kit (Promega, Madison, USA), following the manufacturer´s instructions. DNA quality and quantity were estimated with a NanoDrop™ 1000 spectrophotometer (NanoDrop Technologies, Inc., DE, USA). A pooled DNA sample was obtained for each bivalve species by mixing an equal amounts of DNA (50 ng) extracted from the 15 animals. DNA was stored at -80°C for further analysis.

The V9 region of the 18S rRNA gene (180 bp long) was amplified using the universal eukaryotic-specific primers 1380F/1510R (Amaral-Zettler et al., 2009a). Amplicons were purified using the Qiaquick PCR purification kit (Qiagen). The sequencing library was prepared using the Herculase II Fusion DNA Polymerase Nextera XT Index Kit V2 (Illumina) and the paired-end sequencing (2x300) was performed on an Illumina MiSeq platform (Macrogen, Korea). The raw read sequences obtained were deposited in the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra) under the BioProject accession number PRJNA1083288.

16S rRNA gene PCR amplicon libraries were generated from genomic DNA (Lasa et al., 2019). These primers amplify the V4 hypervariable region. A first target enrichment PCR assay was performed using the 16S rRNA gene conserved primers, followed by a second PCR assay with customized primers that included complementary sequencing adapters (Lasa et al., 2019). Libraries obtained from clam, cockle, and mussel samples collected from July 2016 to June 2018 were sequenced using an Ion Torrent (PGM) Platform (Thermo Fisher Scientific, MA). The raw read sequences obtained were deposited in the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra) under the BioProject accession number PRJNA1082107.

Raw reads were trimmed, and sequencing adapters and low-quality reads were removed. The Microbial Genomics module (Version 25.1) of the CLC Genomics Workbench (Version 25.0.2) was used to perform the analysis.

Following trimming, 18S rRNA gene eukaryotic and 16S rRNA gene prokaryotic reads were clustered into Operational Taxonomic Units (OTUs).

Prior to the eukaryotic clustering, a custom reference database was created, using 18S rRNA gene sequences from the Bivalvia Class (available in SILVA database v.132 clustered at 97%) to classify and remove host eukaryotic reads. Reads were mapped to this database, allowing a similarity and length fraction threshold of 0,9. Unmapped reads (i.e. 18S rRNA gene sequences belonging to non-host eukaryotes) were clustered at a 97% level of similarity into OTUs applying two different databases: SILVA v.132 clustered at 97% and an 18S rRNA gene Custom database generated in a previous work (Ríos-Castro et al., 2021).

The 16S rRNA gene reads were also clustered into OTUs at a 97% similarity level using the SILVA v.132 database (clustered at 97%). For downstream analysis, reads from all the samples were subsampled to 12,400 to ensure comparable sequencing depth.

Singletons and chimaeras were detected and removed. Only reads occurring at least twice in the trimmed dataset were retained for analyses.

Eukaryotic and prokaryotic alpha diversity were analyzed using the Vegan R package. The Shannon and Simpson indices, Fisher’s alpha, and OTU richness were calculated for each bivalve sample across different seasons.

Beta diversity analysis was performed by calculating Bray–Curtis distances and applying Principal Coordinate Analysis (PCoA) to the distance matrices. Based on the Bray-Curtis and Jaccard dissimilarity indices, a permutational multivariate analysis of variance (PERMANOVA) was performed using the Microbial Genomics Module of the CLC Workbench 25.1 software.

Both, alpha- and beta-diversity analyses were performed at the OTU taxonomic level.

Environmental data from a previous study (Ríos-Castro et al., 2022) were used to conduct a meta-analysis. For eukaryotes, we explored the potential presence of parasites in bivalves in relation to their surrounding environment (water and sediment) from the same area. For prokaryotes, we compared diversity patterns among bivalves, sediment, and water, and constructed co-occurrence networks between the bivalve microbiome and habitat-associated bacteria. These networks allowed us to identify keystone taxa, which we subsequently examined for functional profiles.

Differences in prokaryotic diversity indices were statistically tested using ANOVA and Tukey’s multiple comparisons test. However, for the Simpson index, the differences were tested with the Kruskal-Wallis test and Dunn’s multiple comparisons test due to the non-compliance with the assumption of homoscedasticity of variances. A square root transformation of the data for Fisher’s alpha diversity was applied to normalize the data distribution.

Network analyses were conducted to detect potential keystone interactions within the prokaryotic communities of each bivalve and its surrounding environment (mussel-water column, clam-sediment and cockle-sediment). To prevent habitat effects, centered log-ratio (CLR) abundances of prokaryotes were corrected using the habitat filtering (HF) algorithm (Brisson et al., 2019). Spearman correlations with an | ρ | > 0.7 were used to generate the co-occurrence network. Networks were visualized through Gephi using the Force Atlas layout. Singleton nodes without edges and loops were removed from the network. Keystone species were identified by calculating the node degree and betweenness centrality (Degree > 12 and Betweenness centrality > 200). These metrics were obtained using Brandes’ algorithm (Brandes, 2001). Bacterial communities were also identified (Blondel et al., 2008).

Functional activities associated with keystone taxa from each bivalve-habitat network were also predicted. For this purpose, a Phylogenetic Investigation of Communities by Reconstruction of Unobserved States analysis (PICRUSt2 v.2.5.0) was performed based on the abundance OTU tables (including taxa with OTUs composed of more than four reads). Subsequently, the most relevant predicted EC pathways were associated with specific keystone taxa. In this way, some of the detected enzymes could be related to the metabolic pathways in which they are involved.

A total of 3,215,210 raw reads were obtained after sequencing the 18S rRNA gene of the three species of bivalves (20 samples, including mussels, clams and cockles). After trimming, reads were clustered at 97% similarity, resulting in 129 OTUs, using SILVA v132 as reference; and 429 OTUs, using the 18S rRNA gene Custom database as reference. Of this, 61 and 62 OTUs were based on the reference databases, while 68 and 367 were de novo OTUs, respectively (Supplementary Table 1).

For 16S rRNA gene sequencing, 1,385,256 reads were obtained after sequencing of bivalve samples (19 samples including mussels, clams and cockles). After trimming, reads were clustered at 97% similarity, resulting in 2,050 OTUs. In this case, 1,455 OTUs were based on the reference SILVA v132 database, while 595 were de novo OTUs. Following read standardization, 1,295 OTUs were obtained and used for the subsequent analyses (Supplementary Table 1).

In both cases, OTUs with a combined abundance across all samples ≤5 reads were excluded from the analysis to reduce non-representative taxa.

Eukaryotic diversity was evaluated in Cerastoderma edule, Mytilus galloprovincialis and Ruditapes philippinarum sampled from Meira, Ría de Vigo (NW Spain).

For 18S rRNA gene, cockles and mussels showed higher abundance than clams, reaching the highest amount in summer 2017 and winter 2018 (29,033 and 10,683 reads, respectively). In contrast, clam reads remained consistently very low in all sampling points, being not sufficient to draw definitive conclusions. The eukaryotic phyla exhibited completely different profiles among the three bivalves (Figure 1B). While some variations in the microbiome pattern were observed across the sampling points, no clear seasonality was detected (Supplementary Figure 1). Clams had a wide variety of taxa, while mussels and cockles showed clear dominant phylotypes. The Arthropoda phylum was predominant in mussels, reaching more than 95% of the total OTU abundance, whereas the Platyhelminthes phylum was predominant in cockles, with more than 90% of the total OTUs (Figure 1B).

For 16S rRNA gene sequencing, the highest number of reads was observed in clam samples (October 2016) (Figure 1C). Reads from clams and cockles varied, increasing in abundance from February to June 2018. The maximum read values for both bivalves were 89,533 and 71,749, respectively. On the contrary, mussel reads were less abundant and remained relatively constant in all the sampling points, reaching a maximum value of 24,923 reads. The bacterial class profiles shared several predominant classes across the three bivalve species, although their relative abundances differed among hosts (Figure 1C). Specifically, the Gammaproteobacteria class was the most predominant in cockles (80,4%) and mussels (52,4%). In contrast, Mollicutes was the most dominant bacterial class in clams (46,4%), but negligible in the other bivalves (~3%). Bacteroides were also abundant in mussels (22,7%) but not relevant in cockles (8,7%) and clams (abundance below 5%). Finally, the abundance of the Alphaproteobacteria class was also remarkable in mussels (11,3%) (Figure 1C). Differences in bacterial class profiles of clams, mussels and cockles throughout the experiment could not be attributed to apparent seasonality (Supplementary Figure 2).

We calculated several alpha diversity indices to compare the diversity among the three bivalve species across the experiment, and no significant differences were found among bivalves or between seasons (Supplementary Figure 3).

Beta diversity analyses revealed distinct eukaryotic and prokaryotic community composition among bivalves. The diversity of both the 16S rRNA gene and 18S rRNA in clams, mussels and cockles was different (Figure 2). Permutational multivariate analysis of variance (PERMANOVA) based on the Bray-Curtis and Jaccard dissimilarity indices considering the three variables of the sampling procedure (organism, year and season) confirmed significant species-specific differences (Figure 2).

The eukaryotic microbiota included several bivalve parasites. The copepod Mytilicola intestinalis (74-99% of relative abundance) was consistently detected in mussels across all sampling points. In contrast, several platyhelminth species, such as Neohexangritrema zebrasomatis and Intromugil mugilicolus, were detected in cockles for most of the sampling period, representing 14% to 48% of the relative abundance. It was also interesting to detect Bucephalus minimus and Prosorynchoides borealis in September 2017 and June 2018. The genus Mycale was mainly detected during the winter season (January 2017 and February 2018). Furthermore, Perkinsus olseni was detected in the three bivalves and Minchinia merecenariae and Marteilia cochillia were present in cockles, although with low coverage (Figure 3A).

Regarding the bacterial profile, Endozoicomonas was one of the most abundant genera detected in all bivalve samples. Its abundance in mussel tissues decreased from January 2017 to June 2018 (from 94,2% to 12% relative abundance). In contrast, Endozoicomonas increased its prevalence in clams (1% to 81,9%) and cockles (7,3% to 95,2%) in the same period. Interestingly, in June 2018 (the last sampling point of our analysis) its abundance was practically non-existent in all three bivalve species. Vibrio genus was present in all bivalve samples, but it prevailed in cockles with relative abundance values ranging from 52,5% to 82,5%. Moreover, Mycoplasma was highly associated with clams since it was dominant for most of the studied period (62,3 to 98,1%), only being detected in mussels and cockles in June 2018 (25,3%) and October 2016 (23,6%), respectively. Finally, the genera Marinifilium and Psychrilyobacter were especially abundant in all mussel samples, but almost undetectable in clams or cockles. The detection of the genus Bacillus in mussels, specifically Bacillus cereus, should also be mentioned, though at relative abundances below 1% (Figure 3B).

Despite sharing a common environment, the microbiomes of the three species exhibited minimal seasonal variability and strong species-specific signatures, indicative of robust host selection.

Eukaryotic and prokaryotic diversity in the water column and sediment from the same geographical area in Meira, Ría de Vigo (NW Spain) was previously studied (Ríos-Castro et al., 2021, Ríos-Castro et al., 2022). A meta-analysis was conducted to compare the bivalve microbiome with those of the surrounding environment.

The raw abundance of eukaryotic reads obtained from clams, mussels and cockles was considerably lower than that in water column and sediment samples primarily due to the detection hosts 18S rRNA gene (Figure 4). However, the sequencing depth was sufficient to detect several bivalve parasites and potential pathogens. These OTUs corresponded to protozoans, platyhelminths and copepod parasites. Most parasites found in bivalves were also present in the environment, including the water column, sediment, or both. Parallel histological examinations of the same specimens corroborated the presence of some parasites detected by metagenomics.

Protozoans of the genus Perkinsus, specifically Perkinsus olseni and the trematode Bucephalus minimus were commonly found in the two environmental samples and the bivalves. The presence of P. olseni was further confirmed histologically in clams and cockles. Other protozoans, such as Minchinia mercenariae, Marteilia cochilia, or the platyhelminth Prosorhynchoides borealis present in the environment, were only detected in cockles with M. cochilia additionally corroborated by histology. Similarly, the platyhelminth Urastoma cyprinae was found only in clams.

Only a few parasites were present exclusively in bivalves. These included the platyhelminth Macrovestibulum obtusicaudum, the genus Rhinopericolidae and the copepod Herrmannella longicaudata in cockles, and the copepod Mytilicola orientalis in mussels. The trematode Labicola cf. elongata was found only in clams. Finally, some parasites such as Grafilla buccinicola and Mytilicola intestinalis were present in various bivalves, with M. intestinalis confirmed histologically in mussels (Figure 4).

Despite the detection of several parasites by molecular and histological approaches, all examined individuals of the three bivalve species exhibited good health, with no histopathological signs of disease (Supplementary Figure 4).

Significant differences were observed among bivalves, sediment, and water column in terms of bacterial communities (Figure 5). Four different metrics were used to compare bacterial diversity between samples: Shannon index, Simpson index, Fisher alpha and number of species. Sediment and seawater diversities were significantly higher than those found in bivalves (Figure 5A). Principal Coordinate Analysis (PCoA) of bacterial diversity showed clear separation between bivalve and environmental samples. However, a certain degree of similarity existed between bivalve microbiomes on one side and environmental microbiome on the other. Permutational multivariate analysis of variance (PERMANOVA) based on the Bray-Curtis and Jaccard dissimilarity indices revealed significant differences among all groups except for clam and cockle samples (Figure 5B).

There were also marked differences between the bivalve and environmental bacterial profiles (Figure 6A). The five most abundant genera previously described in the three bivalves (Endozoicomonas, Mycoplasma, Vibrio, Marinifilum and Psychrilyobacter) were either not detected or present only rarely in the water column and sediment (in both cases, with relative abundance below 2%).

In environmental samples, OTUs exhibited a more balanced abundance with no single genus standing out. The most abundant genera in the water column were Pseudoalteromonas, Colwellia and Glaciecola. In the case of sediment, the most abundant genus was Arcobacter, followed by Woeseia, Lacinutrix and Aquibacter.

Water and sediment shared 91 OTUs, while all the bivalves and environmental samples shared 54 OTUs. However, only 10 bacterial taxa were exclusively shared among bivalves, indicating a substantial difference between the bivalve microbiome and the environmental microbial ecosystem (Figure 6B; Supplementary Table 2).

A co-occurrence network was performed to examine the interaction between bivalve and environmental microbiome (Figure 7). Significant correlations were detected between some bacterial taxa within each bivalve and their respective environmental habitats.

The bacterial network for the mussel and the seawater column was the most complex, with 135 nodes and 766 edges. Positive correlations (437 co-occurrences) were slightly more abundant than negative correlations (329 co-exclusions) (Supplementary Table 3). The most noteworthy keystone taxa were Vibrio sp. (Degree 60, Betweenness Centrality 3291), which had the highest number of interactions with other taxa. Moreover, Loktanella sp. (Degree 42, Betweenness Centrality 2071), Shewanella sp. (Degree 40, Betweenness Centrality 2654), Algibacter sp. (Degree 39, Betweenness Centrality 1401) and Arcobacter sp. (Degree 35, Betweenness Centrality 1521) were also identified as keystone taxa in the mussel-seawater system. Bacterial communities within the network were clustered by prioritizing stronger interactions within clusters over those with the rest of the network. These groups likely represent functionally related taxa or those sharing similar niches. We detected 10 different bacterial communities, each characterized by its relevant species, interacting within the mussel-seawater environment (Figure 7A).

On the other hand, the network of clam-sediment bacteria consisted of 69 nodes and 197 edges in which co-occurrences (127) were also slightly over co-exclusions (68) (Supplementary Table 3). In this case, Woeseia (Degree 20, Betweenness Centrality 1728), Tenacibaculum (Degree 16, Betweenness Centrality 1893), Halioglobus (Degree 14, Betweenness Centrality 565), Shewanella (Degree 14, Betweenness Centrality 728) and Subsaxibacter (Degree 12, Betweenness Centrality 1238) were identified as keystone genera. Among these keystone genera, the most notable co-occurrences were observed between the Gammaproteobacteria Woeseia sp. and the Bacteroidota Psychroserpens sp. Eight different bacterial communities, with their influenced by the most important species were identified (Figure 7B).

Finally, the cockle and sediment bacterial network comprised 36 nodes and 103 edges. Co-occurrences (61) were slightly higher than co-exclusions (42) (Supplementary Table 3). Lutimonas (Degree 14, Betweenness Centrality 203) and Vibrio (Degree 13, Betweenness Centrality 376) were identified as the keystone species. The most relevant co-occurrences were observed between Lutimonas sp. and Thioclava sp. At the same time, relevant co-exclusions occurred between bacteroides such as Lutimonas sp. - Polaribacter sp. or Vibrio sp. - Aquibacter sp. In this case, eight different bacterial communities and the most important species for every community were also detected (Figure 7C).

PICRUSt2 was used to predict metabolic pathways associated with the most important keystone bacteria identified across all the samples (bivalves, seawater, and sediment). As a result, we obtained 415 metabolic pathways and 2,185 predicted EC numbers (enzyme functions) among all the OTUs. After selecting the most representative EC numbers for each keystone species, eight enzymes were identified as being involved in 23 different metabolic pathways (Figure 8).

Examining the predicted enzymatic functions associated with bacterial profiles, DNA polymerase (EC:2.7.7.7) was one of the most representative enzymes shared among all keystone taxa (except in Algibacter sp.). However, the specific metabolic pathways associated with this enzyme could not be clearly resolved. The same occurred for histidine kinase, DNA helicase and peptidoprolyl isomerase enzymes (EC:2.7.13.3, EC:3.6.4.12 and EC:5.2.1.8), which were shared among fewer bacteria, and their specific metabolic functions also remained unclear.

Nevertheless, NADH: ubiquinone reductase (H+-translocating) enzyme (EC:1.6.5.3) was shared among Loktanella sp., Arcobacter. sp., Woeseia sp., and Tateyamaria sp. and it would be involved in metabolic pathways related to energy metabolism as “aerobic respiration I and II”, “Fe (II) oxidation”, “NAD(P)/NADPH interconversion”, and “NADH to cytochrome bd oxidase electron-electron transfer I”. Otherwise, Acetyl-CoA C-acetyltransferase, enoyl-CoA hydratase and medium-chain Acyl-CoA dehydrogenase (EC:2.3.1.9, EC:4.2.1.17 and EC:1.3.8.7) were three exclusive enzymes detected in Halioglobus sp. These enzymes were predicted to be involved in metabolic pathways related to biosynthesis of secondary metabolites and several routes related to biological processes, carbon biogeochemical cycle, energy metabolism and amino and nucleic acid metabolism.