Corynactis viridis individuals were collected at the “Il Paguro” site by the diving association “Dive Planet” (Rimini, Italy). “Il Paguro” is the wreck of a methane extraction platform located around 12 nautical miles offshore Ravenna, Italy. The wreck extends from a depth of 10 m in its most superficial part to 35 m, and it hosts a great variety of marine wildlife. In order to regulate the scuba activity on the site and preserve its coral reef-like ecosystem, an Italian association responsible for the site (“Associazione Paguro”) has been established and, starting from July the 21st, 1995, the area was declared “Area of biological protection” by the Italian Ministry of Agricultural Resources. In February 2010, “Il Paguro” was declared “Site of community interest” by the Emilia-Romagna region (Rinaldi and Rambelli, 2004). “Associazione Paguro” kindly provided permission to sample C. viridis within the “Il Paguro” site for scientific purposes only and Italian Cost Guard was made aware of the sampling activity. We used satellite maps from NASA (available at https://giovanni.gsfc.nasa.gov/giovanni/) to provide some data on the main environmental variables of the sampling site, including temperature and chlorophyll a to evidence eutrophication occurrence (Ignatiades, 2005). Ten anemone individuals were collected at three different time points in December 2018 (winter sampling), May 2019 (spring sampling) and July 2019 (summer sampling), for a total of 30 individuals (Supplementary Figure 1). For each site, 2 L of seawater were also sampled at the same depth as the anemones, by using previously sterilized polypropylene bottles. Anemones were collected using knives previously sterilized with alcohol, by placing 1–3 rocks hosting different individuals in sterile disposable plastic containers to avoid contamination. The picking of anemones from the respective rocks was performed once returned to the laboratory (within a few hours of sampling at sea), under a vertical laminar airflow cabinet using sterile tweezers and scalpels. Each individual was placed in a single 1.5-mL microtube and stored at −80°C until further processing. Seawater samples were filtered using vacuum filtration with MF-Millipore (Darmstadt, Germany) Membrane filters with 0.45-μm pore size, 47-mm diameter and mixed cellulose esters membrane. Each 2 L sample of seawater was filtered on a single filter, which was rolled on itself using sterilized forceps and placed in a 2-mL microtube, then stored at −80°C until further processing.

Total microbial DNA extraction was performed for each individual anemone (weighing approximately 0.15–0.20 g) using the DNeasy PowerBiofilm Kit (Qiagen, Hilden, Germany) with a modified protocol. In particular, the homogenization step was performed using a FastPrep instrument (MP Biomedicals, Irvine, CA, United States) at 5.5 movements per sec for 1 min and, before the elution step, an incubation of 5 min at 4°C was performed. Seawater microbial DNA extraction was carried out using the DNeasy PowerWater Kit (Qiagen) following the manufacturer’s instructions. DNA samples were quantified using NanoDrop ND-1000 (NanoDrop Technologies, Wilmington, DE, United States) and stored at -20°C until further processing.

Amplification of the V3-V4 hypervariable region of the 16S rRNA gene was performed in 50-μL volumes containing 25 ng of microbial DNA, 2X KAPA HiFi HotStart ReadyMix (Roche, Basel, Switzerland), and 200 nmol/L of 341F and 785R primers carrying Illumina overhang adapter sequences. The thermal cycle consisted of 3 min at 95°C, 30 cycles of 30 s at 95°C, 30 s at 55°C, and 30 s at 72°C, and a final 5-min step at 72°C (Turroni et al., 2016). PCR reactions were purified using Agencourt AMPure XP magnetic beads (Beckman Coulter, Brea, CA, United States). Nextera technology was used to prepare indexed libraries by limited-cycle PCR reaction. After a further clean-up step as described above, libraries were normalized to 4 nM and pooled. The sample pool was denatured with 0.2 N NaOH and diluted to 6 pM with a 20% PhiX control. Sequencing was performed on Illumina MiSeq platform using a 2 × 250 bp paired end protocol, according to the manufacturer’s instructions (Illumina, San Diego, CA, United States).

Raw sequences were processed using a pipeline combining PANDAseq (Masella et al., 2012) and QIIME 2 (Bolyen et al., 2019). High-quality reads (min/max length = 350/550 bp) were binned into amplicon sequence variants (ASVs) using DADA2 (Callahan et al., 2016), as previously described (D’Amico et al., 2019). Taxonomy was assigned using the VSEARCH algorithm (Rognes et al., 2016) and the Silva database (December 2017 release) (Quast et al., 2013). The sequences assigned to the eukaryotic host (chloroplasts and mitochondria) as well as the unassigned ones were discarded. The ASV table was rarefied up to 1733 reads per sample, corresponding to the lowest read number observed among the samples. Alpha diversity was evaluated using three different metrics: Faith’s Phylogenetic Diversity (PD_whole_tree) (Faith, 1992), Chao1 index for microbial richness and number of observed ASVs. The unweighted UniFrac distance was used to build Principal Coordinates Analysis (PCoA) plots.

Statistical analysis was performed using R software version 3.6.1¹. PCoA, bacterial co-abundance groups (CAGs), box-and-whisker plots and Venn diagrams were generated using “Stats”², “Made4” (Culhane et al., 2005), “vegan”³, “boxplot” (Murrell, 2018) and “VennDiagram”⁴ packages. Unweighted UniFrac distance was used for PCoA, and the significance of separation was tested by permutation test with pseudo F-ratio (function “Adonis” in “vegan”). Bacterial groups with the largest contribution to the ordination space were found by using the function envfit of the R package vegan on the family relative abundances. Significant data separation was assessed by Kruskal–Wallis test or Wilcoxon rank-sum test, based on the data. When necessary, p-values were corrected for multiple testing with Benjamini-Hochberg method, with a false discovery rate (FDR) ≤ 0.05 considered as statistically significant.

For graphical representation, anemone microbial phyla whose relative abundance was less than 0.5% in at least 10% of the samples were filtered out under the name “other phyla.” The same was done for microbial families and genera with a relative abundance lower than 2% in at least 10% of the samples (“other families” and “other genera,” respectively). For the generation of CAGs, only bacterial genera present in at least 10% of the samples with relative abundance ≥ 0.5% were considered. Co-abundance was evaluated by the Kendall correlation test and displayed with hierarchical Ward-linkage clustering based on the Spearman correlation coefficient (Claesson et al., 2012). Wiggum plot networks were created using Cytoscape software⁵ (Smoot et al., 2011), where the circle size represents the bacterial abundance and the connections between nodes represent significant Kendall correlations between genera (FDR ≤ 0.05). For Venn diagrams, the ASVs present in at least one sample of each season group were taken into account. We processed all the 4,623 ASVs obtained from DADA2 binning of high-quality reads through the R package VennDiagram, in order to produce a visual output of the ASVs sharing between seasons. The FASTA files of the ASV sequences shared between all seasons or couples of seasons were blasted against the 16S rRNA database through BLASTN algorithm (Altschul et al., 1990).

A total of 30 individuals of C. viridis were sampled in three different seasons in December 2018 (winter), May 2019 (spring) and July 2019 (summer) at the “Il Paguro” site in the North Adriatic Sea. The monthly satellite maps of temperature and surface chlorophyll a from NASA for each sampling period are provided in Supplementary Figure 2. Briefly, seawater temperature during the three sampling times variated between 10–12°C in winter, 15–17°C in spring, and 26–28°C in summer. Chlorophyll a concentration was always > 2 mg/m³ (ranging between 2.2 and 3.6), classifying the sampling region as eutrophic (Ignatiades, 2005; Ferreira et al., 2011). For each time point, 10 individuals and 2 L of seawater were collected at a depth of 20 m. From each sample, the microbiome compositional structure was obtained by NGS sequencing of the V3–V4 hypervariable region of the 16S rRNA gene, resulting in 252,160 high-quality reads and an average of 7,641.212 ± 5,177.348 (mean ± SD) reads per sample. High-quality reads were binned into 4,623 ASVs.

At first, we compared the alpha diversity between C. viridis and seawater communities. Although not significant, our data show a trend toward an overall higher microbial richness in seawater than in anemone, as measured by Chao1 metrics and number of observed ASV (Wilcoxon rank-sum test controlled for multiple testing using FDR, p-value > 0.05) (Supplementary Figure 3).

We next explored the community composition in C. viridis and seawater at different phylogenetic levels (Figure 1). At the phylum level, the C. viridis microbiota is characterized by five dominant phyla, Proteobacteria (mean relative abundance ± SD, 35.3 ± 14.4%), Firmicutes (15.3 ± 20.8%), Tenericutes (15.3 ± 21.3%), Planctomycetes (10.3 ± 8.4%) and Bacteroidetes (10.0 ± 9.1%). Less abundant phyla are represented by Actinobacteria (4.0 ± 4.3%), Cyanobacteria (2.7 ± 4.4%), Verrucomicrobia (1.9 ± 1.6%) and Chloroflexi (1.7 ± 2.2%). The most represented families are Spiroplasmataceae (15.2 ± 21.3%), Planctomycetaceae (9.8 ± 8.1%) and Rhodobacteraceae (7.1 ± 5.5%). Subdominant families are Flavobacteriaceae (4.5 ± 3.0%), Clostridiales family XII (2.9 ± 5.0%), Vibrionaceae (2.9 ± 3.8%) and Clostridiaceae (2.5 ± 4.0%). At the genus level, the dominant taxon is Spiroplasma (15.2 ± 21.3%), while less represented genera are Fusibacter (2.9 ± 5.0%), Blastopirellula (2.5 ± 2.6%) and Pir4 lineage (2.5 ± 3.5%). On the other hand, seawater shows a distinct microbial structure already at the phylum level, being dominated by Proteobacteria (57.1 ± 0.75), Bacteroidetes (19.6 ± 2.3%), Cyanobacteria (7.6 ± 2.7%) and Verrucomicrobia (6.2 ± 2.8%), with Planctomycetes (2.6 ± 2.2%), Actinobacteria (2.5 ± 2.5%) and Firmicutes (1.6 ± 2.8%) as subdominant taxa. The most represented families are Flavobacteriaceae (10.9 ± 3.3%), family I of Cyanobacteria subsection I (7.6 ± 2.8%) and Rhodobacteraceae (6.9 ± 2.9%), while Halieaceae (3.5 ± 3.8%), Verrucomicrobiaceae (3.0 ± 2.0%) and Rhodospirillaceae (3.0 ± 1.2%) are subdominant. Synechococcus (5.3 ± 5.3%) is the most abundant classified genus, while Tenacibaculum (2.5 ± 2.2%) the subdominant one.

In order to highlight the overall compositional differences between C. viridis and seawater microbial communities, a PCoA of the unweighted UniFrac distances was carried out. As expected, the C. viridis microbiome significantly segregates from that of seawater and such a segregation was found to be robust to seasonality (permutation test with pseudo-F ratio, p-value ≤ 0.01) (Supplementary Figure 4). To identify the bacterial families most contributing to the separation, their relative abundance was superimposed on the PCoA plots. According to our findings, the families Planctomycetaceae, Spiroplasmataceae, and Clostridiales family XII are the most characteristic of the C. viridis microbiome, but with a seasonality-specific pattern. Conversely, Flavobacteriaceae and family I of subsection I in the phylum Cyanobacteria are the most distinguishing of the seawater samples.

We tested whether the C. viridis-associated microbial communities varied by season. To do this, we first assessed the existence of significant differences among seasons for alpha diversity. All metrics consistently indicated greater microbial diversity in winter and summer, compared to spring (Kruskal–Wallis test controlled for multiple testing using FDR, p-value ≤ 0.05) (Figure 2A). We next explored seasonal patterns in the community structure by analyzing beta diversity. Interestingly, the unweighted UniFrac-based PCoA revealed a sharp segregation of the microbiome structures through seasonality (permutation test with pseudo F-ratio, p-value ≤ 0.001) (Figure 2B). Conversely, seawater showed only a limited variation across the different seasons. When comparing the relative abundance of C. viridis microbial families among seasons (Figure 3), we found that during the winter, Caldilineaceae, Nitrospiraceae, and Planctomycetaceae were significantly more abundant, whereas the Vibrionaceae family was largely underrepresented. On the other hand, the spring samples were characterized by a higher proportion of the families Rhodobacteraceae, Rhodospirillaceae, Campylobacteraceae, Oceanospirillaceae, and Spiroplasmataceae. Finally, the relative abundance of Marinilabiaceae, family I of Cyanobacteria subsection III, Clostridiaceae 1, family XII of Clostridiales, Lachnospiraceae, Peptostreptococcaceae, Desulfobulbaceae, and Desulfovibrionaceae was significantly higher during the summer.

To gain further insights into the C. viridis-associated microbiome, we explored its topological variation by clustering the bacterial genera into CAGs. Three CAGs were identified, and the prevalence and connections between the represented genera were obtained by a Wiggum plot network analysis (Supplementary Figure 5). Interestingly, we observed a peculiar declination of the C. viridis CAGs based on seasonality. CAGs were named according to the genus showing the higher overabundance in a seasonal-dependent pattern: “Mycobacterium CAG,” “Colwellia CAG,” and “Desulfovibrio CAG” (Figure 4). Winter communities were characterized by the “Mycobacterium CAG”, which also included Nitrospira, Planctomyces, Blastopirellula, Pir4 lineage, OM60 (NOR5) clade and Lysinibacillus. Conversely, the “Colwellia CAG,” including Pseudoalteromonas, Aquibacter, Tenacibaculum, Spiroplasma, Vibrio, Sulfitobacter, Pseudofulvibacter, and Arcobacter, dominated the spring microbiomes. Finally, the summer communities were characterized by the “Desulfovibrio CAG,” including Peptoclostridium, Marinifilum, Cetobacterium, Maribacter, Rhodopirellula, Fusibacter, Streptomyces, and Truepera.

Finally, we investigated whether C. viridis possesses a core microbiome, which is stably associated with the host and thus does not change with the seasons. As shown in Figure 5, only four ASVs were shared among all three seasons, whereas nine ASVs were shared between winter and spring samples, 18 ASVs between winter and summer samples, and seven ASVs between spring and summer samples, for a total of 38 shared ASVs in the three sampling seasons. According to a BLAST analysis, the four ASVs shared among the three seasons were all assigned to the bacterial family Spiroplasmataceae (Table 1). More detailed information on BLAST alignment for each ASV is provided in Supplementary Table 1.