Nine specimens of A. lixula and nine P. lividus were sampled from each of the sampling locations: La Palma (Canary Islands, Spain) and Cala de Sant Francesc, Blanes (Girona, Spain) in October 2023 and February 2024, respectively (Figure 1 and Supplementary Table S1). In La Palma, nine P. lividus individuals were collected from Fuencaliente, while nine A. lixula specimens were sampled from La Bajita in the eastern part of the island. These individuals were transported in specialized containers to the Observatorio Marino de Cambio Climático (OMACC) in Fuencaliente. Specimens from Blanes were transported and processed at the Universitat de Barcelona. Specimens were separated into individual tanks containing seawater collected from their respective sampling sites during the sample collection process. Once in the laboratory, coelomic fluid was extracted using a 5 mL sterile syringe through the peristomial membrane pre-filled with 2 mL of sterile anticoagulant buffer [80% Calcium/Magnesium Free Artificial Sea Water (CM-ASW), osmolarity 1200 mOsM+20% EDTA stock solution (13.53 g/L)]. The fluid was briefly centrifuged to separate coelomic fluid and coelomocytes, the latter were stored in absolute ethanol and at -20 °C. The remaining coelomic fluid was filtered using sterile Sartorius Minisart™ filters with a pore size of 0.22 µm to retain microbial cells. The filters were fixed with absolute ethanol and stored at -20 °C. Fecal pellets were collected after a 24-hour fasting period to ensure defecation, and samples were preserved in absolute ethanol at -20 °C and used as a non-invasive proxy to investigate the digestive tract-associated microbiota. In La Palma, coelomic fluid (CF), coelomocytes (C), and fecal pellets (F) were collected, whereas in Blanes only coelomic fluid samples were obtained. For La Palma, coelomocytes and feces were collected from multiple individuals, but only five samples per species were selected for sequencing. As controls, two samples of 5 mL of anticoagulant buffer used during coelomocyte extraction were filtered through 0.22 µm membranes and stored in absolute ethanol and at -20 °C. Following processing, the specimens were returned alive to their original sampling locations.

DNA was extracted from all collected samples using DNeasy™Blood & Tissue kit (QIAGEN, Maryland, USA) following a slightly modified manufacturer’s protocol. Coelomocytes and fecal pellets were briefly ethanol-dried, followed by overnight lysis in 20 µL proteinase K and 180 µL ATL buffer at 56 °C. For coelomic fluid samples, half of the filter was cut in small pieces, and an additional 40 µL ATL buffer was added to ensure complete coverage. All DNA was further purified using OneStep PCR Inhibitor removal kit (ZymoBIOMICS™, Zymo Research, Irvine, CA, USA). DNA quantification was evaluated using a Qubit dsDNA HS (High Sensitivity) Assay Kit with the Qubit^® Fluorometer (Massachusetts, USA). A total of 58 samples were selected for sequencing including an artificial seawater and DNA extraction blank as negative controls. We employed Illumina MiSeq sequencing (2×300 bp paired-end reads) to analyse the V3–V4 hypervariable region of the 16S rRNA gene, which was amplified using universal primers 341F (5’-CCTACGGGNGGCWGCAG-3`) and 785R (5’-GACTACHVGGGTATCTAATCC-3`) (Klindworth et al., 2013). Library preparation and sequencing were conducted at the Centre for Genomic Research in Barcelona, Spain. The sequences are available in GenBank Bioproject PRJNA1252121 with accession numbers SRR33214878-SRR33214973.

DNA sequences pre-processing and bioinformatic analyses were performed using the software environment Quantitative Insights into Microbial Ecology version 2019.4 (Bolyen et al., 2019). QIIME2 provides a software environment, data standards, and tool wrappers allowing seamless interoperability between tools used for microbial community analysis. The bioinformatic pipeline used is available at https://github.com/leaschmutsch/microbiota_characterization. Primers were removed without mismatch tolerance using the Cutadapt plugin (Martin, 2011). We used the DADA2 plugin (Callahan et al., 2016) for paired-end merging, trimming parameters were set at 260bp for forward and reverse reads, chimera removal (using consensus method) and clustering into Amplicon Sequence Variants (ASVs). To create a phylogenetic tree, we used the ‘align-to-tree-mafft-fasttree’ pipeline from the ‘q2-phylogeny plugin’. This step consists of doing a multiple-sequence alignment, then filtering the alignment to remove highly variable positions and applying FastTree to generate a phylogenetic tree from the masked alignment. Taxonomic assignment of ASVs was performed against the SILVA v138.1 reference database (Quast et al., 2012). Singletons and ASVs taxonomically assigned to eukaryotes were discarded for statistical analysis in R programming environment (R Core Team, 2024).

Potential contaminants were removed with the decontam R package v1.4 (Davis et al., 2018) using the prevalence method, and the extraction and PCR negative controls. The ASV table was used to calculate rarefaction curves using phyloseq R package v1.46 (McMurdie and Holmes, 2013). We optimized the rarefaction thresholds to balance sequencing depth and sample retention, selecting 7,000 sequences per sample to support robust statistical comparisons. Relative abundance barplots were generated at the family, order and phylum level using the plot_composition function in the phyloseq R package.

We then estimated alpha diversity (taxonomic richness) and beta diversity (community composition) based on pairwise dissimilarity. To test differences, we applied a planned comparisons framework, focusing on a set of predefined, biologically meaningful hypotheses. Seven comparisons were defined a priori, including differences: (i) between A. lixula and P. lividus within the same body compartment; (ii) among compartments within a species; and (iii) between geographic locations for the same species and compartment (see Supplementary Table S2 for details). These comparisons were encoded in a contrast matrix, assigning weights (−1, 0, + 1) to each level depending on the comparison of interest (Crawley, 2012). Each column of the contrast matrix corresponds to a particular comparison. Because our experimental design was partially unbalanced (species × body compartment × location), we set up a dummy variable to uniquely identified each factor combination, and each sample was then assigned to one of these levels. For example, a sample of coelomic fluid of A. lixula from la Palma would have been coded as “Arb.CF.Palma”. This resulted in eight levels for the 56 samples considered, excluding the negative controls. This ensured that all planned comparisons could be tested directly without unnecessary subsetting or multiple testing of the same data. As some contrasts were non-orthogonal (analogous to collinear predictors), we used the inverse of the transposed contrast matrix to calculate the fixed-effects design matrix. By implementing planned contrasts through contrast matrices, and by explicitly modelling fixed effects (species, compartments, locations) and random effects (sample identity), we efficiently used the degrees of freedom in our design and tested only the comparisons of biological interest. Importantly, because all planned contrasts were defined within the models, corrections for multiple testing were not required. This approach provides robust and interpretable estimates while accounting for the hierarchical structure of the sampling design.

Alpha diversity metrics of bacterial communities associated with the sampled sea urchins were estimated using Observed ASVs and Shannon diversity index with microbiome R package v1.24 (Lahti and Shetty, 2018) and Faith’s PD index (phylogenetic diversity) with picante R package v1.8 (Kembel et al., 2010). Statistical comparisons of alpha diversity were conducted using a permutational linear model (PLM), implemented via the lmp() function from the lmPerm R package v2.1.0 (Wheeler et al., 2025), which performs analysis of variance through permutation of residuals. In these models, species, body compartment, and location were treated as fixed effects, while sample identity was included as a random effect whenever multiple compartments were measured from the same sea urchin. This explicitly accounts for repeated measures and avoids pseudoreplication. This method was chosen as it accommodates complex factorial designs while relaxing the assumptions of parametric ANOVA, which are often violated in microbiome datasets.

For the beta diversity analysis of community composition, Unweighted UniFrac, a presence/absence-based metric sensitive to phylogenetic relationships, and Weighted UniFrac, which incorporates abundance data, were calculated using the phyloseq package. For data visualization a Principal Coordinates Analysis (PCoA) was used with the UniFrac distance matrices. Permutational multivariate analysis of variance (PERMANOVA) (Anderson, 2005), was employed to assess differences in microbial community, using both UniFrac distance matrices. Statistical significance was determined with 999 permutations implemented via adonis2 function in the vegan v2.6 package (Oksanen et al., 2007). When repeated measures from the same sea urchin were present, permutations were constrained within individuals to avoid inflating type I error rates. We assessed the assumption of homogeneity of dispersions and analyzed multivariate beta diversity patterns among samples by calculating the average distance to group centroids using the betadisper function from the vegan R package.

Unique and shared taxa between species and between sample types of each species were counted at multiple taxonomic levels and visualized using Venn diagrams with ggvenn R package. Differential abundance analysis was conducted using Analysis of Composition of Microbiomes with Bias Correction (ANCOM-BC) from the Ancombc R package (Lin and Peddada, 2020; Nearing et al., 2022) on unrarefied ASV tables and taxonomic assignments were reported at genus, family and phylum level. Log fold changes (LFCs) in taxa abundances between samples were computed, with p-values adjusted using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995) to control for false discovery rates (FDR). Results were visualized using ggplot2 R package (Wickham, 2016) as heatmaps showing taxa-specific LFCs when three factors were compared (i.e sample types within each species) and as waterfall plots when two factors were compared (i.e species and locations), highlighting significant differences between sample groups.

To infer the potential functionality of microbial taxa, the Functional Annotation of Prokaryotic Taxa (FAPROTAX) software (Louca et al., 2016) was used, which maps prokaryotic taxa to ecological functions based on curated literature. For each functional group in a given sample, the count reflects the cumulative occurrences of the number of ASVs associated with that group. The functional table was normalized to relative proportions, and filtered to remove functional groups with zero abundance. Mean and standard error of functional abundances were calculated for each factor of the desired comparison, and visualized with stacked bar plots comparing functional groups with abundances ≥0.01%. To test for significant differences in predicted functional profiles, we applied PERMANOVA using Jaccard dissimilarity and 999 permutations to test for the effects of species, sample types or location as fixed factors depending of the comparison. When repeated measures from the same sea urchin were present, permutations were constrained within individuals. The functional profile table, derived from FAPROTAX annotations, was normalized to relative abundances prior to analysis. To identify the most influential functions driving observed differences, we performed a SIMPER (Similarity Percentage) analysis with vegan package in R, which estimates the contribution of individual functional groups to dissimilarities between groups.

To study the differences in microbial composition between the two species, all sample types of each species from La Palma were pooled, to capture the overall microbial profile of each species (Supplementary Figure S2). Alpha diversity comparisons between A. lixula and P. lividus revealed significant species-level differences for Observed ASVs and Shannon diversity. P. lividus harbored a significantly higher number of observed ASVs and Shannon diversity compared to A. lixula. In contrast, Faith’s phylogenetic diversity did not differ significantly between species (Figure 2A; Supplementary Table S5).

Beta diversity analysis based on PCoA of Unweighted and Weighted UniFrac distance matrices revealed that the first principal coordinates explained 49.4% and 45.5% of total variation, respectively, and clearly separated samples according to species (Figure 2B; Supplementary Figure S3). PERMANOVA models confirmed significant differences in microbial community composition between species (Unweighted: p = 0.001; Weighted: p = 0.001). Group dispersion analyses (PERMDISP) showed significant variation in dispersion using the Unweighted UniFrac matrix (p = 0.001), but not for the Weighted UniFrac matrix (p = 0.356) (Supplementary Tables S6). These dispersion differences were mainly driven by variation among tissue types, particularly in A. lixula, where fecal samples were more distinct compared to coelomic fluid and coelomocyte samples (Figure 2B).

The bacterial identified comprised 29 phyla across all sample types, locations, and species. The most abundant bacteria phyla were Proteobacteria, Bacteroidetes and Fusobacteriota comprising the 32.11%, 27.02% and 22.20% of relative abundance across all samples, respectively (Figure 2C; Supplementary Figure S4). The comparison of bacterial taxa between the two species revealed significant differences across several taxonomic levels and their relative abundances as identified through ANCOM-BC analysis (Supplementary Table S7, S8). Interestingly, more taxa were exclusively found in A. lixula than P. lividus (Supplementary Figure S5). The most abundant microbial phylum in P. lividus was Fusobacteria and particularly the family Fusobacteriacea (31.43%). In A. lixula the most abundant phylum was Bacteroidota (36.7%), within this phylum, the families Flavobacteriaceae and Cryomorphaceae were significantly more abundant in A. lixula (15.6% and 18.9%, respectively) than in P. lividus (0.2% and 5.9%, respectively). The phylum Proteobacteria was more predominant in P. lividus microbiota than in A. lixula (26% and 16.8%, respectively). Within this phylum, families such as Psychromonadaceae and Colwelliaceae were more abundant in P. lividus (3.4% and 2.8%) than in A. lixula (<0.5%), whereas other families such as Shewanellaceae, Pseudoalteromonadaceae, and Rhodobacteraceae were significantly more abundant in A. lixula. Five phyla, Desulfobacteriota, Verrucomicrobiota, Cyanobacteria, Planctomycetota and Campilobacterota showed higher significant abundances in P. lividus compared to A. lixula, exceeding log fold change of 2.3. Particularly, taxa within the phyla Desulfobacteriota and Verrucomicrobiota were present in P. lividus at higher abundances than 5% and in A. lixula in relative abundances lower than 1.4%. The phylum Firmicutes was more abundant in A. lixula (9.74% relative abundance) compared to P. lividus (3.80%). Notably, the genus Candidatus Hepatoplasma (LFC = 4.46), within Firmicutes, was significantly more abundant in A. lixula, reaching up to 70% relative abundance in some samples, whereas its maximum relative abundance in P. lividus samples was only 7.5%. Moreover, the phylum Spirochaeta represented 11% of average relative abundance in A. lixula microbial profiles while only 1.6% average relative abundance in P. lividus (Figure 2D; Supplementary Table S8).

Functional inference using FAPROTAX and Jaccard-based dissimilarity revealed that both host species significantly differed in the presence/absence of microbial functional traits. PERMANOVA results showed effects of species (p = 0.004) (Supplementary Table S9). SIMPER analysis identified both widespread and rare functional groups contributing to differences between A. lixula and P. lividus microbiomes. General functions such as chemoheterotrophy, fermentation, and aerobic chemoheterotrophy were shared across species. In A. lixula, microbial communities were enriched in nitrogen-cycling functions, including nitrate/nitrogen respiration (p < 0.002), as well as metabolism pathways such as methanol oxidation and methylotrophy (p < 0.001). These functions correspond to the higher relative abundance of families such as Shewanellaceae, Rhodobacteraceae, and Flavobacteriaceae, which are classified by FAPROTAX as contributors to nitrogen cycling, and in the case of Rhodobacteraceae, also to methylotrophy and methanol oxidation. Additionally, A. lixula harbored a greater abundance of bacteria associated with intracellular parasitism and potential pathogenicity (e.g., intracellular parasites, human pathogens; p < 0.01). In contrast, P. lividus microbiomes were functionally distinct due to a higher prevalence of sulfur-related metabolism, including sulfate respiration and respiration of sulfur compounds (p < 0.001). These functions correspond to the higher relative abundance of taxa such as Desulfobacterota and Campylobacteriota, which are classified by FAPROTAX as sulfate reducers and sulfur-respiring bacteria (Figure 2E; Supplementary Table S10).

Among the microbial communities of the three A. lixula sample types (coelomic fluid, coelomocytes, and feces) from La Palma, fecal samples exhibited significantly higher alpha diversity. Planned contrast confirmed that feces harbored significant greater diversity than coelomic fluid (Shannon diversity: p < 0.001 and Observed ASVs: p < 0.001) and coelomocytes (Shannon: p=0.001 and Observed: p<0.001), while phylogenetic diversity did not differ significantly across tissues (p = 0.094) except between coelomocytes and feces (p=0.043). Differences between coelomic fluid and coelomocytes were not significant for any diversity metric (p > 0.15). (Figure 3A; Supplementary Table S11). Ordination of samples using PCoA based on Unweighted and Weighted UniFrac distances revealed two distinct clusters, one corresponding to fecal samples and another encompassing coelomocyte and coelomic fluid samples, indicating a separation in microbial community composition between these two groups (Figure 3B; Supplementary Figure S6). Venn diagrams also supported distinctiveness of fecal sample compared to the samples from the coelomic compartment (Supplementary Figure S7). PERMANOVA revealed significant differences in assemblage composition between feces and coelomic fluid for Unweighted UniFrac (p = 0.040), but not for Weighted UniFrac (p = 0.078). Similar trend between feces and coelomocytes for Unweighted UniFrac (p = 0.017), but not for Weighted UniFrac (p = 0.543). Differences between coelomic fluid and coelomocytes were not significant for either distance metric. Tests of homogeneity of multivariate dispersions (PERMDISP) showed no significant differences among groups (all p > 0.12; Supplementary Table S12), indicating that the observed compositional separation reflects true differences in community structure rather than unequal variability. Taxa within the phylum Proteobacteria, such as the families Vibrionaceae (20.5%), Alteromonadaceae (3.5%) and Shewanellaceae (4.7%), the phylum Fusobacteriota such as the family Fusobacteriaceae (27.9%) and the family Desulfocapsaceae (2%) within the phylum Desulfobacteriota were significantly more abundant in feces than coelomic fluid and coelomocytes, with log-fold changes exceeding 2.86 (Supplementary Table S13, S14). Conversely, the phylum Spirochaetota (family Spirochaetaceae) was primarily enriched in coelomic fluid and coelomocytes, reaching 15% of relative abundance in coelomic samples and only 0.02% in fecal samples. Within the phylum Bacteroidota, the family Flavobacteriaceae was more abundant in coelomic samples (33.9%) than in fecal samples (1.3%). Additionally, the phylum Firmicutes (order Izemoplasmatales), even at low abundances, were found in an order of magnitude higher in coelomic compartments than in feces. The heatmap visualization of the log-fold changes highlighted distinct clustering patterns, with fecal microbiota displaying stronger differential abundance patterns than the other tissue types (Figure 3C; Supplementary Tables S14).

Among the microbial communities of the three P. lividus sample types, fecal samples of P. lividus supported higher microbial diversity than coelomic fluid and coelomocytes, whereas coelomic fluid and coelomocytes did not differ significantly. Specifically, feces exhibited significantly higher Shannon diversity and Observed ASVs, while Faith’s phylogenetic diversity did not differ between compartments (p > 0.453) (Figure 3D; Supplementary Table S15). The PCoA based on Unweighted and Weighted UniFrac distances revealed two main clusters in P. lividus, one corresponding to fecal samples and another encompassing coelomocyte and coelomic fluid samples (Figure 3E; Supplementary Figure S8). Venn diagrams also supported distinctiveness of fecal sample compared to the samples from the coelomic compartment (Supplementary Figure S9). PERMANOVA indicated no significant differences between coelomic fluid and coelomocytes across any distance metric (all p > 0.526). However, feces differed significantly from coelomic fluid in Unweighted UniFrac (p = 0.049), but not in Weighted UniFrac (p > 0.073). Tests of homogeneity of multivariate dispersions (PERMDISP) revealed no significant differences in variability between compartments across metrics (all p > 0.06), indicating that the observed differences reflect true compositional changes rather than dispersion effects (Supplementary Table S16). Significant differential abundant taxa were observed across sample types using ANCOM-BC analysis, particularly when comparing between coelomic fluid and fecal samples. However, coelomocyte samples notably presented a similar presence and abundance of microbial taxa as coelomic fluid. Microbial taxa within the phylum Proteobacteria were consistently more abundant in fecal samples (41.6% of relative abundance) than in coelomic fluid (22.05%) and coelomocytes (20.03%). Within this phylum, the families Alteromonadaceae and Psychromonadaceae were significantly more abundant in the fecal samples, with a log fold change exceeding 1.7. Additionally, the family Vibrionaceae was found at twice the abundance in fecal samples (21.5% relative abundance) compared to the other sample types (average 11.5%). In contrast, taxa within the phylum Firmicutes, particularly the order Izemoplasmatales, were significantly more abundant in the coelomic compartment. Besides that, the phylum Bacteroidota is one of the most abundant across all sample types of P. lividus, representing an average of 25% of relative abundance. The family Flavobacteriaceae presented significantly higher abundance in the coelomic and coelomic fluid samples. The relative abundance of the phylum Fusobacteriota in coelomic fluid (32.7%) and coelomocytes (35.5%) samples was twice as high compared to fecal samples (13%) (Figure 3F; Supplementary material S17, S18).

To gain deeper insight into the predicted microbial functional profiles across body compartments, ASV-assigned taxa from both A. lixula and P. lividus were annotated using the FAPROTAX database. In both species, distinct functional signatures were observed across tissue types, with fecal samples consistently enriched in diverse metabolic pathways compared to the more host-associated coelomic fluid and coelomocytes. In A. lixula, PERMANOVA did not detect statistically significant differences among sample types (p = 0.058), and the relevant pairwise comparisons were also non-significant (Supplementary Table S19). In P. lividus, however, functional profiles significantly differed among sample types (p = 0.047), with strong pairwise differences between feces and coelomocytes (p = 0.004) and coelomocytes and coelomic fluid (p = 0.027) (Supplementary Table S20). SIMPER analyses revealed consistent functional enrichments in fecal samples of both species, including aerobic chemoheterotrophy and multiple steps of the nitrogen cycles, such as nitrate reduction, nitrate/nitrite respiration, functions largely assigned by FAPROTAX to Rhodobacteraceae and Shewanellaceae. In A. lixula, fecal samples also exhibited enrichment in sulfur-related metabolic functions, including respiration of sulfur compounds, sulfate respiration, dark oxidation of sulfur compounds, and dark sulfite oxidation which correspond to the presence of Desulfobacteraceae and other sulfur-reducing lineages. In P. lividus, feces additionally featured photoheterotrophy and anoxygenic photoautotrophy, associated with members of the Rhodobacteraceae, suggesting potential light-driven metabolic processes. In contrast, the coelomic fluid and coelomocyte microbiota in both species exhibited reduced functional diversity and a higher prevalence of traits associated with animal symbionts, intracellular parasites, and putative pathogens (e.g Vibrionaceae) (Figure 4; Supplementary Table S21, S22).

To assess geographic variation in microbiota, we focused exclusively on coelomic fluid samples, as this was the only compartment sampled for both species at both locations. Geographic comparisons showed that alpha diversity varied by location only in A. lixula (Figure 5A; Supplementary Table S23). Samples from Blanes exhibited significantly higher Shannon diversity than those from La Palma (p = 0.03), while Observed ASVs and Faith’s phylogenetic diversity did not differ significantly between locations (p > 0.428). In contrast, no significant geographic differences were detected for P. lividus across any of the alpha diversity metrics (all p > 0.48; Supplementary Table S24). However, Venn diagrams at different taxonomic levels showed that there were more unique taxa in Blanes than in La Palma in both species (Supplementary Figure S10).

Geographic location also influenced microbial community composition, although patterns differed between species (Figure 5B; Supplementary Figure S11). In A. lixula, PERMANOVA indicated significant location effects across both distance metrics, Weighted UniFrac (p = 0.001), and Unweighted UniFrac (p = 0.008) (Supplementary Table S25A). By contrast, in P. lividus, location effects were significant for Unweighted UniFrac (p = 0.001), but not for Weighted UniFrac (p = 0.182; Supplementary Table S26A). In both species, PERMDISP tests showed no significant differences in within-group variability (p > 0.18), confirming that observed effects represent genuine compositional shifts rather than differences in dispersion (Supplementary Tables S25B, S26B).

The differential abundance analysis in A. lixula revealed that the phyla Proteobacteria (including the family Vibrionaceae) and Firmicutes (order Izemoplasmatales and family Mycoplasmataceae) were more abundant in Blanes (15.17% and 21.9% average relative abundance, respectively) than in La Palma (8.15% and 11.21%). In contrast, the phyla Bacteroidota and Spirochaeta were more abundant in La Palma, exceeding log fold change of 1 (Figure 5C; Supplementary Figure S12; Supplementary Table S27). However, no significant differences were found with ANCOM-BC analysis (Supplementary Table S28). In P. lividus the phylum Proteobacteria was more abundant in Blanes (47.9% average relative abundance) than in La Palma (20.23%), where the order Rickettsiales, and the families Colwelliaceae, Rhodobacteracea, Shewanellaceae and Vibrionaceae were significantly more abundant, exceeding log fold changes of 2.5. Microbial taxa within the phyla Firmicutes (families Hungateiclostridiaceae, Christensenellaceae and the order Izemoplasmatales) and Verrucomicota (Kiritimatiellaceae and the class Lentisphaeria “P.palmC41”) were found significantly more abundant in individuals from La Palma compared to Blanes. The phylum Bacteroidota was more abundant in La Palma (19.86% relative abundance) than in Blanes (11.46%), particularly the family Cryomorphaceae (Supplementary Tables S29, S30).

To explore the influence of geographic location on microbial functional potential, we used FAPROTAX annotations of ASV-assigned taxa from coelomic fluid samples of A. lixula and P. lividus collected from La Palma and Blanes (Figure 6). In A. lixula, PERMANOVA revealed no significant differences in predicted functions between locations (p = 0.303). However, SIMPER analysis identified functional groups contributing to location-based variation, with higher representation of nitrate and nitrogen respiration, ureolysis, and dark sulfur-related metabolisms (e.g., dark sulfite/sulfur oxidation) in Blanes samples, functions largely assigned to Colwelliaceae, Rhodobacteraceae and Shewanellaceae (Supplementary Table S31). In contrast, P. lividus exhibited strong functional differentiation by location, with PERMANOVA showing a significant effect of geographic origin (p = 0.002) (Supplementary Table S32). SIMPER results highlighted functional enrichments in Blanes, particularly in aerobic chemoheterotrophy, nitrate reduction, and sulfur-related metabolisms such as sulfate respiration and dark hydrogen oxidation which were associated with Desulfobacteraceae. Moreover, functions related to intracellular parasitism appeared more abundant in Blanes, largely associated with Rickettsiales.