Bacterioplankton did not exhibit any time or treatment-specific pattern.
We first assessed the patterns observed in the bacterioplankton community, a factor known to covariate with fish skin mucus communities [28, 31, 32]. ANOSIM tests performed on GUniFrac distances (Table 1), and Kruskal-Wallis tests performed on Simpson index (Supplementary data, Figure S1, Tables S2, S3) showed that phylogenetic structure and alpha diversity of bacterioplankton did not exhibit any time or treatment-specific pattern. This result suggests that bacterioplankton was not significantly associated to the microbial community restructuring observed in clownfish and anemones from physical (PI) and remote (RI) interaction groups from T1 to T3, as detailed below.
Gradual convergence of the epithelial microbiome of clownfish and anemones.
GUniFrac dissimilarity analyses. Analysis of the anemones’ epithelial microbiota (Figure 2a) shows that prior to contact with clownfish, GUniFrac dissimilarity between test and control groups was minimum (0.49 ± 0.02) and not significantly different between RI and PI (MW test p = 0.85) (T0: after three weeks of acclimation). At T1, 15 minutes after the clownfish test individuals were transferred from the fish control tank system into their respective two-tanks systems for remote interaction (RI) (i.e. six biological replicates of one anemone tank connected with one fish tank), and after the first 15 min of physical contact between physical interaction (PI) clownfish individuals with their respective anemone (i.e. six biological replicates of physical interaction), dissimilarity between test (remote and physical interaction) and control anemones was significantly higher (mean 0.3955 ± 0.0502 for RI, MW test, Bonferroni corrected p = 0.01 and 0.55 ± 0.01 for PI, MW test, Bonferroni corrected p = 0.003) relatively to that of T0. Then, the dissimilarity between test and control anemones remained high and stable during the interaction period (mean from T1 to T3: 0.58 for RI and 0.4357 for PI± 0.03), as there was no significant difference between time nor group (T1 to T3). From T4 (one week after PI / RI clownfish individuals were retrieved), to T5 (two weeks after PI / RI clownfish individuals were retrieved), dissimilarity between test and control anemones reached a plateau for RI anemones, whereas increased again for PI anemones to become significantly more differentiated from controls than RI anemones (MW test, Bonferroni corrected p = 0.01). Patterns of dissimilarity based on Thetayc distances were similar to that of GUniFrac (Mantel test, r = 0.8885 P = 0.001).
Regarding clownfish skin microbiota (Fig. 2b), a similar pattern to that of the anemones occurred from T0 to T3: GUniFrac dissimilarity at T0 between PI / RI clownfish test and control groups was minimum (0.040 ± 0.01 for RI and 0.43 ± 0.01 for PI) and significantly different between RI and PI (MW test, Bonferroni corrected p = 0.035) prior to fish contact with their respective anemone. At T1, after the first 15 min of physical contact with their anemone for PI clownfish test individuals and 15 min after RI clownfish test individuals were placed into their respective tank systems, dissimilarity between PI / RI test and control clownfish was significantly higher (0.64 ± 0.01 for RI, MW test Bonferroni corrected p = 2E-16 and 0.56 ± 0.01 for PI, MW test Bonferroni corrected p = 4.1E-15) relatively to that of T0. Then, the dissimilarity between PI / RI test and control clownfish increased further during the interaction period (T1 to T3) to reach 0.68 ± 0.01 for RI and 0.69 ± 0.01 for PI at T3. Interestingly, a significantly higher dissimilarity with the control group was detected at both T1 and T2 for RI (MW tests, Bonferroni corrected p = 3.8E-4 at T1 and p = 2.6E-11). From T4 (one week after PI / RI clownfish individuals were retrieved and moved back to the control clownfish water system), to T5 (two weeks after), dissimilarity between PI / RI test and control clownfish groups remained stable (mean from T3 to T5: 0.69 for RI and 0.67 for PI) and not significantly lower compared to that of the contact period (T1-T2-T3) (MW test, p = 1). Patterns of dissimilarity based on GUniFrac distances were similar to that of Thetayc index (Mantel test, r = 0.8406 P = 0.001).
Finally, regarding GUniFrac dissimilarity between fish and anemone microbiota (Fig. 2c; Table S1), it was similar at T0 in all groups (mean 0,65 ± 0.01 for RI, 0.66 ± 0.01 for PI and 0.66 ± 0.01 for control, MW test p = 1). Then, at T1, T2 and T3, the dissimilarity in both PI and RI test groups was significantly lower than in control groups (Mann-Whitney test, Bonferroni corrected p < 0.015). At T2, dissimilarity dropped down to 0.58 ± 0.01 in PI, whereas it increased up to 0.73 ± 0.01 in RI, although being significantly lower than dissimilarity between control fish and anemones (MW test, Bonferroni corrected p = 0.015). From T2 to T3 test groups dissimilarity in PI and RI converged to a mean value of 0.65, significantly below the stable dissimilarity (0.79 ± 0.01) observed between fish and anemone control groups (MW test, Bonferroni corrected p < 6.3E-13). At T4, one week after fish-anemone pairs' separation, the dissimilarity values of PI and RI were still significantly lower than in control fish anemone pairs (MW tests with Bonferroni correction, p = 1.2E-05 for PI and p = 2.7E-11 for RI). At T5, two weeks after fish-anemone pairs' separation, the dissimilarity values converged to that of control for both RI and PI test groups, although being still significantly lower to that of the control group (MW tests with Bonferroni correction, p = 1.3E-04 for PI and p = 5.9E-05 for RI).
This partial recovery is likely explained by the stability of the dissimilarity between PI / RI test and control clownfish groups from T3 to T5 (Fig. 2b) and in RI anemones (Fig. 2a). In addition, the dissimilarity between the host microbiota and the bacterioplankton (Fig. S2) was never significantly different between the test and control groups. Thetayc dissimilarity between fish and anemone microbiota exhibited a similar trend to that of the GUniFrac for both RI and PI test groups during the whole experiment (i.e. from T0 to T5). (Mantel test, r = 0.835 P = 0.001).
Bacterial taxa mostly associated to fish-anemone epithelial microbiota convergence peaked after two weeks of interaction.
ASV tables detailing the most differentially abundant taxa at each time point in each group have been provided in the Supplementary Materials (Tables S4 and S5). At T0, there were only 5 differentially abundant taxa between interaction and control fish-anemone pairs. At T1, after the first 15 minutes of clownfish-anemone interaction, differentially abundant taxa increased to 10 ASVs. At T2, after one week of interaction, differentially abundant taxa doubled to reach 21 ASVs. At T3, after two weeks of interaction, differentially abundant taxa peaked at 30 ASVs. At T4, one week after separation of interaction fish-anemone pairs, the number of differentially abundant ASVs remained at 30. At T5, two weeks after separation of interaction fish-anemone pairs, the number of differentially abundant taxa dropped to 17 ASVs. From T2 to the end of the experiment (T5), three ASVs (2, 49, 177) matching to Cellulophaga tyrosinoxydans strain EM41 (95% identity, 100% coverage, 1.31E-119 to 2.81E-121 e-values) exhibited an interesting dynamic: they peaked at T3, with the three highest Bonferroni corrected p-values, with a fold change ranging from 9 to 12, then decreased gradually after separation of fish-anemone pairs: fold change ranging from 6 to 11 at T4, and from 3 to 8 at T5, relatively to fish-anemone control group. Therefore, these three ASVs related to Cellulophaga tyrosinoxydans were further analyzed in terms of abundance dynamics in water, sea anemone and clownfish mucus.
Differential abundance analysis (DESeq2) onCellulophaga sp. in water, sea anemone and clownfish. The monitoring of the three ASVs (2, 49, 177) related to C. tyrosinoxydans, which were differentially abundant from T2 to T5 (DSeq2, Fig. 3, Table 3) was decomposed in terms of host community (sea anemone, clownfish, water), experimental groups and time. At T0, Cellulophaga sp. counts were both low and variable across experimental groups and host communities. ASVs with log2-normalized fold-change over 1 and FDR corrected p-values < 0.0001 were kept (Table 2).
Tank system water: From T0 to T1, the abundance of Cellulophaga sp. dropped in all experimental groups to become undetectable except for FC, where only ASV 2 was significantly higher to PI (8.9 fold change). At T2, Cellulophaga sp. was still undetectable in anemone control group, and dropped under the detection threshold in clownfish control group. On the contrary, Cellulophaga sp. counts increased for the three ASVs both in PI and RI tank water to become statistically higher than in AC and FC groups (8.6 to 13.6 fold changes). At T3, the three ASVs were still undetectable in both control group water, whereas peaking in both PI and RI groups (9.7 to 13.7 fold changes). From T4 to T5, one and two weeks after clownfish retrieving from PI and RI tank systems, the three ASVs counts decreased gradually (5.2 to 13.7 fold changes at T5) in both PI and RI tank system water.
Sea anemone epithelium: From T0 to T1, Cellulophaga sp. counts dropped under the detection threshold in the three experimental groups hosting anemones (AC, PI, RI). At T2, Cellulophaga sp. counts increased for ASVs 49 and 177 in PI and RI anemones to become significantly higher than in control (7.7 and 9.5 fold changes). At T3, the counts of the three ASVs peaked in PI and RI groups, and were still significantly higher than in control (8 and 13 fold changes). From T4 to T5, one and two weeks after clownfish retrieving from PI and RI tank systems, Cellulophaga sp. counts decreased quickly: ASV 2 was no more differentially abundant, and ASVs 49 and 177 dropped from 5.8 to 10.2 fold changes at T4, and were no more significantly different from their control at T5.
Clownfish skin mucus: At T0, the three Cellulophaga sp. related ASVs counts were low and comparable between the three clownfish groups, which had shared the same tank system water for the three weeks acclimation period. From T0 to T1, Cellulophaga sp. counts remained low and comparable between the three clownfish groups, despite the transfer of PI and RI individuals to their respective PI and RI tank systems hosting anemones. At T2, Cellulophaga sp. counts of ASVs 49 and 177 increased in PI and RI to become significantly higher than in control (7.7 to 9.4 fold changes). At T3, the counts of the three ASVs peaked in PI and RI groups, and were still significantly higher than in control (10.7 to 13.3 fold changes). At T4, one week after PI and RI clownfish were reintroduced into the fish control water system, the counts of the three ASVs remained high and significantly higher than in control fish (7.1 to 12.9 fold change), despite sharing the same tank system water. At T5, two weeks after PI and RI clownfish reintroduction into the fish control water system, the three ASVs counts remained high only in PI clownfish (5 to 10.6 fold changes), whereas ASVs 2 and 49 dropped drastically in RI clownfish, both of them being no more significantly higher than in control fish.