Microbiomes of clownfish and their symbiotic host anemone converge before their first physical contact

Background One of the most charismatic, and yet not completely resolved example of mutualistic interaction is the partnership of clownfish and its symbiotic sea anemone. The mechanism explaining this tolerance currently relies on the molecular mimicry of clownfish epithelial mucus, which could serve as camouflage, preventing the anemone's nematocysts' discharge. Resident bacteria are known as key drivers of epithelial mucus chemical signature in vertebrates. A recent study has proposed a restructuration of the skin microbiota in a generalist clown fish when first contacting its symbiotic anemone. We explored a novel hypothesis by testing the effect of remote interaction on epithelial microbiota restructuration in both partners. Methods With metataxonomics, we investigated the epithelial microbiota dynamic of 18 pairs of percula clownfish (Amphiprion percula) and their symbiotic anemone Heteractis magnifica in remote interaction, physical interaction and control groups for both partners during a 4-week trial. Results The Physical and Remote Interaction groups’ results evidence gradual epithelial microbiota convergence between both partners when fish and anemone were placed in the same water system. This convergence occurred preceding any physical contact between partners, and was maintained during the 2-week interaction period in both contact groups. After the interaction period, community structure of both fish and anemone’s epthelial community structures maintained the interaction signature 2 weeks after fish–anemone pairs’ separation. Furthermore, the interaction signature persistence was observed both in the Physical and Remote Interaction groups, thus suggesting that water-mediated chemical communication between symbiotic partners was strong enough to shift the skin microbiota durably, even after the separation of fish–anemone pairs. Finally, our results suggest that fish–anemone convergent microbiota restructuration was increasingly associated with the parallel recruitment of three Flavobacteriaceae strains closely related to a tyrosinase-producing Cellulophaga tyrosinoxydans. Conclusions Our study shows that bacterial community restructuration, in the acclimation process, does not only rely on direct physical contact. Furthermore, our results challenge, for the first time, the traditional unidirectional chemical camouflage hypothesis, as we argue that convergence of the epithelial microbiota of both partners may play essential roles in establishing mutual acceptance. Video abstract Fish−anemone symbiotic relationship. Supplementary Information The online version contains supplementary material available at 10.1186/s40168-021-01058-1.


Abstract Background
One of the most charismatic, and yet not completely resolved example of mutualistic interaction is the partnership of clown sh and its symbiotic sea anemone. The mechanism explaining this tolerance currently relies on the molecular mimicry of clown sh epithelial mucus, which could serve as camou age, preventing the anemone's nematocysts' discharge. Resident bacteria are known as key drivers of epithelial mucus chemical signature in vertebrates. A recent study has proposed a restructuration of the skin microbiota in a generalist clown sh when rst contacting its symbiotic anemone. We explored a novel hypothesis by testing the effect of remote interaction on epithelial microbiota restructuration in both partners.

Methods
With metataxonomics, we investigated the epithelial microbiota dynamic of 18 pairs of percula clown sh (Amphiprion percula) and their symbiotic anemone Heteractis magni ca in remote interaction, physical interaction and control groups for both partners during a four weeks trial.

Results
Physical and Remote Interaction groups' results evidence gradual epithelial microbiota convergence between both partners when sh and anemone were placed in the same water system. This convergence occurred preceding any physical contact between partners, and was maintained during the two-weeks interaction period in both contact groups. After the interaction period, community structure of both sh and anemone's epthelial community structures maintained the interaction signature two weeks after shanemone pairs separation. Furthermore, the interaction signature persistence was observed both in Physical and Remote Interaction groups, thus suggesting that water-mediated chemical communication between symbiotic partners was strong enough to shift the skin microbiota durably, even after the separation of sh-anemone pairs. Finally, our results suggest that sh-anemone convergent microbiota restructuration was increasingly associated with the parallel recruitment of three Flavobacteriaceae strains closely related to a tyrosinase-producing Cellulophaga tyrosinoxydans.

Conclusions
Our study shows that bacterial community restructuration, in the acclimation process, does not only rely on direct physical contact. Furthermore, our results challenge, for the rst time, the traditional unidirectional chemical camou age hypothesis, as we argue that convergence of the epithelial microbiota of both partners may play essential roles in establishing mutual acceptance.

Page 3/24
The interaction of anemones and clown sh is a charismatic example of mutualistic partnership [1], in which the anemone protects the clown sh against predators [2], while the clown sh provides the anemone's endosymbiotic zooxanthellae algae with excreted nutrients (ammonia, sulfur, and phosphorus) [3]. This mutualism is contingent upon a protective mechanism for clown sh against the anemone's nematocyst discharge. Numerous studies and reviews that attempted to identify the protective mechanism in different clown sh species, have highlighted two main non-exclusive hypotheses: either clown shes bene ciate from an innate protective mechanism in their skin mucus, and/or they need to coat their body with anemone mucus [4,5] as a chemical camou age. Interestingly, during clown shanemone acclimation (i.e. prior to rst physical contact), the epithelial mucus immunological pro le of clown sh changes to mimic that of the anemone [6][7][8]. Thus, clown sh epithelial mucus is suspected to act as a "chemical camou age" preventing "not-self" recognition associated with nematocysts' discharge [4,5,9,10]. Most importantly, transferred amino acids from the anemone mucus to the clown sh skin mucus were measured after few hours of physical and remote interaction [11], thus suggesting that chemical modi cation of clown sh skin mucus starts before physical interaction with the anemone. Given that skin microbial communities are important drivers of the chemical and immunological pro les of vertebrates' epithelia, which modulate host-parasite interactions [12], including in sh [13,14], it is relevant to investigate to which extent clown sh anemone symbiosis translates into epithelial microbiome shifts in both partners. In addition, as clown sh are known to hover near/above their partner anemone, before rst contacting its tentacles, it is relevant to test whether sh and anemone microbiome shifts may precede their physical contact. To date, the structure of clown sh skin microbiome after contact with their symbiotic anemone has only been investigated partly, without anemone control groups during the contact phase, nor testing the remote interaction [15,16]. For instance, Pratte et al. (2018) [15] conducted an investigation targeting a generalist clown sh species (A. clarkii) known to have an innate mechanism of protection against nematocyst discharge. However, their experimental design included physical interaction groups that shared the same water ow than the sh control group [15], and thus, the control group was constantly in remote interaction with sh-anemone pairs. Therefore, the general mechanism underlying the microbial community changes observed in the initiation of the mutualist partnership is still unknown. First, it is essential to characterize the dynamic of anemone microbiome to assess its involvement in anemone-clown sh mutual acceptance. Then, the study of control groups with remote interaction (i.e. absence of physical contact) between both hosts has also never been done, and yet, this data is essential to determine if the putative restructuring of the transient epithelial microbiome of both partners actually participates in establishing mutual acceptance, or if it is merely an artifact of the physical contact between both hosts. In order to detect evolutionary relevant microbiota changes regarding mutualism, we focused on Amphiprion percula, a clown sh species exhibiting narrow host speci city, which is mainly associated with the Heteractis magni ca sea anemone. A. percula is a common and locally abundant clown sh species, inhabiting lagoons and seaward reefs from the western Paci c off Queensland and Melanesia (e.g. northern Great Barrier Reef, northern New Guinea, New Britain, Solomon Islands and Vanuatu) [17,18]. In nature, A. percula lives in small groups consisting in one breeding pair and few non-breeders, all of them being in close interaction with one host sea anemone. pelagic to demersal habitat: after a 10-15 days dispersive oceanic phase at the larval stage, and following the completion of metamorphosis two days later, a sedentary reef phase starts as young juveniles actively look for settling into a host sea anemone [19]. Clown sh lay their eggs on a rock in the close vicinity of sea anemones in such way it is hypothesized that eggs are exposed to anemone metabolites without been stung by nematocysts [20]. From hatching to the rst contact of the sedentary reef phase, clown sh species are sensitive to sea anemone nematocysts [21]. Therefore, we targeted young juvenile stage of so-called "naïve" individual (i.e. no prior contact with anemone) to model as much as possible the natural acclimation process step in controlled conditions. Our objective was to test two hypotheses: (1) is anemone-clown sh mutualistic partnership associated with a signi cant restructuration of epithelial microbiotas in both partners during acclimation? and (2) does this skin microbiota restructuration precedes physical contact between partners? Then, we aimed to characterize the microbial taxa driving the observed community dynamics. To achieve our goal, four experimental groups were compared: anemone control, sh control, physical interaction (i.e. sh and anemone in the same tank), and remote interaction (i.e. sh and anemone in different tanks, both connected to the same water ow). We tested the hypothesis that there would be a gradual convergence of the skin mucus microbiota structure of both symbiotic partners during the interaction period, which would not only rely on physical contact.

Methods
Clown sh and sea anemone rearing Eighteen H. magni ca and 18 captive bred naïve A. percula juveniles from a tropical sh distributer (Reef Solution, Laval, Qc, Canada; no prior contact with anemone) were acclimated in four recirculated aquatic systems (RAS) during three weeks in 20 L tanks with separated water ow to avoid any chemical contact prior to the experiment. H. magni ca is known to be almost unable to survive more than a few months in captivity [22]. Therefore, particular care was taken for optimizing rearing conditions. Tanks of synthetic sea water (Reef Crystal, Instant Ocean) were heated to 27°C and illuminated 12h/24h with bright lightning provided by pairs of Fluval Sea Marine 2.0 LED Light Fixture 48" ramps, each providing 1350 Lumens and 15 000 K. Nitrates were maintained below 5-10 mg/L. In each 20 L tank, water was pulsed with a 180 L/h water pump, in addition to the intake water from the recirculated system. Clown sh were fed daily, and anemones were fed three days a week with mysis shrimps, directly pipetted on the oral disc. Food waste was removed daily by syphoning water. During the experiment, anemones did not show any obvious signs of stress. Furthermore, anemones survived six to nine extra months in several reef tanks after the end of the experiment.
There were four experimental groups, two controls and two tests. Each experimental group consisted in a single RAS system connected to a deep sand bed biological lter, and containing six biological replicates (i.e. six tanks). Two control groups: anemone control (AC), sh control (FC), and two test groups: physical interaction (i.e. sh and anemone in the same tank, PI), and remote interaction (i.e. sh and anemone in different tanks, all being connected to the same water ow, RI) (Fig. 1). To minimize bacterioplankton taxonomic drift due to independent water ow in each experimental group, 30% water changes were conducted each day in both interaction groups with a water mix from both control groups. To prevent any induction of "remote interaction" between anemone and clown sh control groups, we did not add the water mix (control anemone and control sh) in the control tanks. Following the acclimation period (Fig.  1a), the interaction period between clown sh and anemone for physical and remote interaction groups lasted two weeks (Fig. 1b), after which clown sh and anemones were separated for a two weeks resilience period (Fig. 1c). To mitigate possible effect of sh transfer from sh control to physical and remote interaction units, all the sh specimens were manipulated, including the six control sh individuals, which were moved from one tank to the other in the control tank system.

Host microbiota and water sampling
Seven sampling steps were as follow: T0, at the end of a three weeks acclimation period; T1, after the rst 15 min of physical interaction between clown sh and anemone (PI), 15 min after transfer of sh to the recirculation system containing anemones (RI), and immediately in both control groups; T2 and T3, respectively one and two weeks after initial interaction (T1); T4 and T5, respectively one and two weeks after sh-anemone pairs separation from physical and remote interaction groups (T3). Note that at T1, the six PI individuals experienced an extended RI (from 0 to 24h) before being in physical contact with their respective anemone. After capture, the skin mucus of all shes was immediately sampled by gently rubbing a sterile cotton swab on ≈50% of the total surface (upper half) of the right side of each sh outside of the water as in [23]. The same area was sampled on each sh to standardize the sampling zone. Anemone epithelium mucus was sampled in the same way, by temporarily lowering the tank water, therefore allowing gentle cotton swab rubbing on tentacles out of the water. To characterize the bacterioplankton community of each group, at every sampling time 2 L of tank water were collected in sterile Nalgene bottles and immediately ltered on 0.22 μm membranes using a peristaltic pump. We used four bacterioplankton replicates per experimental group.
DNA extraction, libraries preparation, and 16S amplicons sequencing DNA extraction of epithelial mucus from clown sh and sea anemone, as well as 0.22 μm membranes from water samples, was performed using the Qiagen® Blood and Tissue Kit according to the manufacturer's instructions. The fragment V3-V4 of the 16S rRNA was ampli ed in a two-step dualindexed PCR approach speci cally designed for Illumina instruments by the Plateforme d'Analyses Génomiques (IBIS, Université Laval, Quebec city, Canada). The rst PCR was performed with 16S region speci c primers which were tailed on the 5' end with part of the Illumina TruSeq adaptors. A second PCR was performed to attach remaining adaptor sequence (regions that anneal to the owcell and library speci c barcodes). Please note that primers used in this work contain Illumina speci c sequences protected by intellectual property (Oligonucleotide sequences © Illumina, Inc. All rights reserved. Derivative works created by Illumina customers are authorized for use with Illumina instruments and products only. All other uses are strictly prohibited.). The following oligonucleotide sequences were used for two rounds of ampli cation:

Reads de-noising and ASV identi cation
Bioinformatics' processing was undertaken using dada2 as reported in [24]. 4,770,388 raw reads were quality ltered with a truncation to 270 base pairs. Before proceeding further, we assessed the quality of the sequence data. We observed a better taxonomic accuracy with the forward reads covering the V3 region of the 16S rRNA gene. Recent studies have shown that the V3 represents the best choice for the pro ling of complex microbial communities [25], for increased and more accurate estimates of community richness and diversity [26]. Thus, 270 bp reads covering V3 were used for downstream analyses, as in [25,26,27]. The standard dada2 pipeline was used, with a maximum of 2 expected errors (maxEE) and with default parameters as per version 1.14.1, except for the "dada" step where all samples were pooled for ASV inference [24]. The lowest number of reads in a sample post-dada2 was 6 736 and the maximum was 54 137 with a median number of reads across all samples of 23 762.

Taxonomic annotation
Taxonomic annotation of amplicon sequence variants (ASV) was performed by using blastn matches NCBI "16S Microbial" database. As the NCBI database for 16S sequences is updated more frequently than other sources [28], it matched our requirements for exhaustive information about lesser-known taxa, while minimizing ambiguous annotations. Matches above 99% identity were assigned the reported taxonomic identity. Sequences with no matches above the identity threshold were assigned taxonomy using a lowest common ancestor method generated on the top 50 blastn matches obtained. This method is closely inspired from the LCA algorithm implemented in MEGAN [29].

ASV abundance tables normalization processing
From the initial ASV abundance table supplied by dada2, two normalized tables were made for downstream analyses. A relative abundance normalized table was constructed and used in the following analyses: ANOSIM, alpha diversity, beta diversity. This normalization was made to produce relative abundance tables by dividing amplicon counts by the total counts per sample. No rarefaction was made, as random rarefaction is known to discard valid data, to introduce random errors and to reduce the overall quality of the data [30]. Raw data distribution was not transformed prior to DESeq2 analysis since the DESeq2 pipeline expects raw data. This is a crucial requirement for the statistical model (DESeq2 internally performs normalization by multiplying raw counts by normalization factors). ANOSIM Initial overall testing for phylogenetic structure dynamics in bacterioplankton with respect to our experimental conditions was done using the ANOSIM method from the R "vegan" package (v2.5.6). All ve possible conditions (Control Fish/Control Anemone/Physical/Remote), and all 6 time points were used. The test was permuted 999 times and p-values were corrected with false discovery rate (FDR, qvalues).

Alpha diversity metrics
Simpson index was calculated using the R "vegan" package (v2.5.6). For the Faith PD index, the sequences were rst aligned using the R "DECIPHER" package (v2.10.2) using -16 to -18 for gap opening penalties and -1 to -2 for gap extension penalties. Aligned sequences were then used with the R "phangorn" package (v2.5.5) to generate a phylogenetic tree (using a JC69 substitution model and then a Neighbor-Joining construction method). The tree and the ASV table were nally used with the R "picante" (v1.8.2) package to calculate the Faith PD index. Bonferroni corrected pairwise Mann-Whitney U tests, hereafter named MW tests, were performed to assess statistically signi cant changes in alpha diversity between experimental groups.

Beta diversity metrics
GUniFrac distance comparison. Generalized UniFrac (GUniFrac) distance analyses were performed with the following sample description: Control/Interaction groups and six times (T0, T1, T2, T3, T4, T5). GUnifrac distance matrices were calculated using the same above-mentioned phylogenetic tree and the R "GUnifrac" package (v1.1) with an alpha parameter of 0.5, as recommended by the package authors. The distance matrix was then split to separate each condition (Control/Physical/Remote) in order to allow for the plotting and comparison of said conditions across all six time points (using ggplot2). GUniFrac distances were preferred over other distance metrics for two main reasons: rst, it is a measure of phylogenetic distance, which captures ecological changes since functional repertories are coupled with taxonomy in bacteria (reviewed in [31]). Second, GUniFrac distance uses both the phylogenetic divergence and an extra parameter α controlling the weight on abundant lineages in order to mitigate the in uence of highly abundant lineages over other lineages in UniFrac distance computing [32]. Therefore, GuniFrac is expected to capture more ecologically relevant changes in microbial communities relatively to other phylogeny-informed distances metrics. As our GUnifrac distances evolve over time, a longitudinal approach by using locally weighted scatterplot smoothing (LOWESS) curves was used to model the taxonomic pro le changes. LOWESS is a non-parametric and computationally demanding but robust and exible regression method that uses locally weighted polynomials to t a smoothed curve to scatterplot data [33]. In addition to smoothed curves, a con dence interval of 99% was computed using a t-based approximation method. As with alpha diversities above, Bonferroni corrected pairwise MW tests were performed to evaluate statistically signi cant changes in GUniFrac distance between experimental groups and time points.
Other Beta diversity metrics. GUniFrac distances were nonetheless compared with Mantel tests to commonly used non-phylogenetic indices such as Thetayc, Bray and Curtis and the Jaccard distance. The Thetayc index is a function of species proportions from both the shared and non-shared species. In addition, in the Thetayc index, the shared species proportions in each community are compared one-toone (instead of a sum of the abundances of all shared species in the Bray-Curtis index, which gives no indication on which species are shared, and if the abundances of the shared species are similar or not).
As a result, the Thetayc index places more weight on those shared species, which have similar species proportions in both communities. All Mantel tests were all signi cant (Supplementary data), therefore con rming the overall consistency of dissimilarity measures.

Detection of differentially abundant ASVs
Two differential abundance analyses were performed with the R "DESeq2" package (v1.22.2). The rst analysis aimed to identify high impact ASVs that responded with a strong statistical signal to physical/remote interactions. To do so, de novo ASV abundances (i.e. ASV determined by DADA2 without regard to taxonomic annotation) in clown sh and anemone were monitored during the whole experiment (from T0 to T5) using differential abundance analysis (DESeq2) to identify bacterial taxa that were mostly associated to sh-anemone epithelial microbiota convergence. Then, de novo ASV abundances of clown sh and anemones were combined, PI and RI groups were combined as an interaction group, and clown sh and anemone controls were combined as a control group. ASVs with log2-normalized foldchange over 1 between the ASV abundance in interaction versus control groups, and with a Bonferroni corrected p-value < 0.05, were kept for further analysis (Table S3). The second analysis aimed to validate the differential abundance of the ASVs identi ed during the rst analysis (three ASVs related to Cellulophaga tyrosinoxydans), but on every possible time/group/sample type contrast possible, including all four experimental groups (RI, PI, CF, CA), sample type (bacterioplankton, sea anemone and clown sh), and time points (T0, T1, T2, T3, T4, T5). Thresholds used were an FDR adjusted p-value of 0.0001 and a fold-change of 1 between compared groups.

Results
Bacterioplankton did not exhibit any time or treatment-speci c pattern.
We rst assessed the patterns observed in the bacterioplankton community, a factor known to covariate with sh skin mucus communities [31,34,35]. ANOSIM tests performed on GUniFrac distances (Table 1), and Kruskal-Wallis tests performed on Simpson index (Supplementary data, Figure S1, Table S1) showed that phylogenetic structure and alpha diversity of bacterioplankton did not exhibit any time or treatmentspeci c pattern. This result suggests that bacterioplankton was not signi cantly associated to the microbial community restructuring observed in clown sh and anemones from physical (PI) and remote (RI) interaction groups from T1 to T3, as detailed below.
Gradual convergence of the epithelial microbiome of clown sh and anemones.
GUniFrac dissimilarity analyses. Analysis of the anemones' epithelial microbiota (Figure 2a) shows that prior to contact with clown sh, GUniFrac dissimilarity between test and control groups was minimum (0.49 ± 0.02) and not signi cantly different between RI and PI (MW test p = 0.85) (T0: after three weeks of acclimation). At T1, 15 minutes after the clown sh test individuals were transferred from the sh 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 sh tank), and after the rst 15 min of physical contact between physical interaction (PI) clown sh individuals with their respective anemone (i.e. six biological replicates of physical interaction), dissimilarity between test (remote and physical interaction) and control anemones was signi cantly 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 signi cant difference between time nor group (T1 to T3). From T4 (one week after PI / RI clown sh individuals were retrieved), to T5 (two weeks after PI / RI clown sh individuals were retrieved), dissimilarity between test and control anemones reached a plateau for RI anemones, whereas increased again for PI anemones to become signi cantly 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 clown sh skin microbiota (Fig. 2b), a similar pattern to that of the anemones occurred from T0 to T3: GUniFrac dissimilarity at T0 between PI / RI clown sh test and control groups was minimum (0.040 ± 0.01 for RI and 0.43 ± 0.01 for PI) and signi cantly different between RI and PI (MW test, Bonferroni corrected p = 0.035) prior to sh contact with their respective anemone. At T1, after the rst 15 min of physical contact with their anemone for PI clown sh test individuals and 15 min after RI clown sh test individuals were placed into their respective tank systems, dissimilarity between PI / RI test and control clown sh was signi cantly 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 clown sh 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 signi cantly 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 clown sh individuals were retrieved and moved back to the control clown sh water system), to T5 (two weeks after), dissimilarity between PI / RI test and control clown sh groups remained stable (mean from T3 to T5: 0.69 for RI and 0.67 for PI) and not signi cantly 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 sh and anemone microbiota ( Fig. 2c; Table S2), 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 signi cantly 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 signi cantly lower than dissimilarity between control sh 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, signi cantly below the stable dissimilarity (0.79 ± 0.01) observed between sh and anemone control groups (MW test, Bonferroni corrected p < 6.3E-13). At T4, one week after sh-anemone pairs' separation, the dissimilarity values of PI and RI were still signi cantly lower than in control sh 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 sh-anemone pairs' separation, the dissimilarity values converged to that of control for both RI and PI test groups, although being still signi cantly 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 reconvergence of the distances between sh and anemones of experimental groups towards the values observed for the control group is likely explained by the stability of the dissimilarity between PI / RI test and control clown sh 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. S3) was never signi cantly different between the test and control groups. Thetayc dissimilarity between sh 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 sh-anemone epithelial microbiota convergence peaked after two weeks of interaction.
Tables detailing the most differentially abundant ASVs at each time point in each group have been provided in the Supplementary Materials (Table S3). At T0, there were only 5 differentially abundant ASVs between interaction and control sh-anemone pairs. At T1, after the rst 15 minutes of clown shanemone interaction, differentially abundant ASVs increased to 10. 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 sh-anemone pairs, the number of differentially abundant ASVs remained at 30. At T5, two weeks after separation of interaction sh-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 sh-anemone pairs: fold change ranging from 6 to 11 at T4, and from 3 to 8 at T5, relatively to sh-anemone control group. Therefore, these three ASVs related to Cellulophaga tyrosinoxydans were further analyzed in terms of abundance dynamics in water, sea anemone and clown sh mucus.
Differential abundance analysis (DESeq2) onCellulophaga sp. in water, sea anemone and clown sh. 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, clown sh, 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 for subsequent analyses (Table 2, Fig. 3).
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 signi cantly 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 clown sh 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 clown sh 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 signi cantly 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 signi cantly higher than in control (8 and 13 fold changes). From T4 to T5, one and two weeks after clown sh 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 signi cantly different from their control at T5.
Clown sh skin mucus: At T0, the three Cellulophaga sp. related ASVs counts were low and comparable between the three clown sh 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 clown sh 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 signi cantly 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 signi cantly higher than in control (10.7 to 13.3 fold changes). At T4, one week after PI and RI clown sh were reintroduced into the sh control water system, the counts of the three ASVs remained high and signi cantly higher than in control sh (7.1 to 12.9 fold change), despite sharing the same tank system water. At T5, two weeks after PI and RI clown sh reintroduction into the sh control water system, the three ASVs counts remained high only in PI clown sh (5 to 10.6 fold changes), whereas ASVs 2 and 49 dropped drastically in RI clown sh, both of them being no more signi cantly higher than in control sh.

Discussion
Remote interaction between naïve clown sh and sea anemone triggered epithelial microbiota convergence. Our results from remote and physical interaction groups revealed that prior to the rst physical contact, both clown sh and anemone epithelial microbiotas converged from T1. This convergence was not accompanied with a shift of the bacterioplankton pro le according to alpha diversity dynamics between T0 and T1 (Supplementary results). Contrastingly, an increase of the Simpson index in both PI and RI sh mucus (Fig. S1) may suggest that the interaction between symbiotic partners involved a quick restructuration of clown sh microbiota in terms of richness and evenness.
Furthermore, after the rst 15 minutes of physical interaction, PI anemones microbiota exhibited a decrease of the Simpson index (Fig. S1, Supplementary data). Interestingly, this restructuration in terms of richness and evenness in RI anemones was not detected at T1 but at T2. This apparent delay is most likely explained by the fact that the six PI anemones were basically sampled along a 24h period from T1, corresponding to the time needed for each PI sh to freely initiate physical contact with their anemone (see methods). Therefore, both PI sh and anemones experienced extended remote interaction compared to the RI group.
It is interesting to note that the restructuration in terms of richness and evenness in RI anemones detected at T2 coincided with the signi cant increase of phylogenetic diversity, but not evenness in RI sh, (Supplementary data, Figure S1). This co-occurence of shifts in both RI sh and anemone coincides with their dissimilarity increase from T1 to T2, which precedes the maximum convergence between RI and PI at T3 (GUniFrac distances, Fig. 2c). Moreover, the gradual occurrence and persistence of convergent restructuration between sh and anemone epithelial microbiota in the remote interaction (RI) group from T1 to T3 (Fig. 2c) suggests that microbial community restructuration in both partners does not only rely on physical interaction and occurs before such direct contact. Therefore, the convergent restructuration between symbiotic partners likely starts as soon as the clown sh skin mucus is exposed to sea anemone chemical compounds that are released into the surrounding water, as observed by Schlichter (1975Schlichter ( , 1976) [9,10]. In addition, this clown sh / anemone epithelial microbiota convergence has to be paralleled with Schlichter (1972) [36], which quanti ed the amount of tritium labelled amino acids transferred from the anemone to the naïve clown sh in physical interaction. The pattern of tritium activity from labelled anemones corresponded to that of clown sh adapted to those anemones. Thus, this experiment suggested that the sh skin mucus composition changed during acclimation to resemble that of the anemone [36]. Interestingly, the naïve control clown sh in this experiment was in remote interaction, as it was isolated form the anemone using a perforated plexiglass sheet. To this respect, a residual tritium activity was measured in these naïve control sh, likely due to sh mucus coating with water released anemone metabolites. In a more recent survey focusing on clown sh / anemone epithelial microbiota, Roux et al. (2019) [16] observed in a closely related clown sh species (A. ocellaris) that taxonomical composition of sh skin microbiota in physical contact with sea anemone (H. magni ca) was closer to those of the anemone when compared to control clown sh. However, as there was neither control anemone nor replication of experimental groups, their observations need further validation.
Separation between symbiotic partners revealed contrasting microbiota resilience dynamics. After separation of symbiotic partners (T4-T5), clown sh and sea anemone interaction group microbiota exhibited a contrasting response according to their respective interaction group. RI anemone and clown sh microbiota divergence with that of their respective controls stabilized from T4 to T5. Contrastingly, PI anemone and sh microbiota changed in opposite ways with that of their respective controls: PI anemone divergence was still increasing from T4 to T5, while PI sh remained stable from T3 to T4, despite being moved to the control sh water system since one week, then decreasing from T4 to T5. Then, when considering sh-anemone microbiota dissimilarity, both RI and PI partially converged to that of their controls. In addition, control sh-anemone microbiota dissimilarity signi cantly decreased between T3, T4 and T5, therefore also converging towards that of both contact groups, although remaining signi cantly higher (Fig. 2c, Table S2). Interestingly, the convergence of RI with the control was delayed of one week comparatively to PI as it occurred between T4 and T5 instead of between T3 and T5 for PI group. This observation may suggest that the delayed convergence between RI partners relatively to PI partners translated in a delayed convergence with control during the resilience period (T4 to T5).
Taken together, two observations show that neither environmental water nor time, two factors that are known to drive bacterial community shifts [31,34,35], did play any major rule in reshaping the clown sh microbiota in interaction groups: i the lack of full convergence between ex-interaction PI and RI sh with control sh during the resilience period, despite sharing the same water system during two weeks, and ii the maximum dissimilarity index measured between contact and control sh at T3, which remained stable during the resilience period.
Furthermore, the sh skin microbiota signature of C. tyrosinoxydans related ASVs in physical/remote interaction with sea anemone, remained detectable two weeks after clown sh-anemone pair's separation (Fig. 1c, 2c; Table 2). This observation can be paralleled with what was observed in previous pioneering experiments in controlled conditions from Schlichter (1972) [36] focusing on anemone mucus proteins and antigens in A. clarkii: those molecules persisted in clown sh skin mucus after clown sh-anemone pairs' separation, either after physical or remote interaction [8,9]. In some cases, it is possible that remote interaction is su cient to elicit mutual acceptance between clown sh and the sea anemone as observed for A. xanthurus: 24-48 h of 5-7 cm remote interaction with the sea anemone Antheopsis sp in the eld conferred a protection against nematocysts of many other anemone species [10]. Cellulophaga sp.are the main mucus symbionts involved in sea anemone clown sh contact. The microbiota convergence observed between sh and anemone during the contact period starting at T1 was followed from T2 by the gradual parallel recruitment of three initially rare Flavobacteriaceae symbionts closely related to Cellulophaga tyrosinoxydans. This parallel recruitment of Cellulophaga sp. peaked at T3, two weeks after sh-anemone pair contact and fade out with contrasting dynamics in sh and anemones from their separation (T4-T5). In anemones, two out of three ASVs were still signi cantly more abundant than in controls at T4, but no more at T5, whereas being signi cantly more abundant in water. Contrastingly in sh, all the three ASVs were still signi cantly more abundant than in controls at T4, as well as at T5 for PI clown sh only, despite sharing the same FC water. Therefore, these results suggest that physical interaction during two weeks (T1-T3) exerted a more sustainable imprinting in sh skin microbiota than remote interaction, where two out of three C.tyrosinoxydans related ASVs counts converged to that of control sh. Given that transfer of anemone mucus proteins and antigens to clown sh skin mucus and their persistence after clown sh-anemone pairs' separation was documented [8], the persistence of C. tyrosinoxydans strains two weeks after clown sh-anemone pairs' separation at least demonstrates their relationship with the biochemical imprinting of the clown sh-anemone mutualistic association, and possibly suggests that those ASVs are tightly involved in the remote communication between clown sh and its anemone host. Interestingly, Flavobacteriaceae were observed to be the most enriched bacterial family in clown sh hosting anemones in natural conditions, relatively to non-hosting anemones, [37]. This work is relevant to our ndings as they included the same anemone species used in our experiment, Heteractis magni ca. Unfortunately, clown sh epithelial microbiota was not characterized, so that the involvement of Flavobacteriaceae in clown sh anemone microbiota convergence is still to be con rmed in nature. Regarding clown sh microbiota, in the study of Roux et al. (2019) [16] on the H. magni ca / A. ocellaris system, Flavobacteriaceae appeared to be enriched in clown sh in physical contact with their anemone, one week after the start of the contact phase. However, the signi cance of this increase was not reported, and sh-anemone microbiota convergence was mainly explained by three other bacterial families (Pseudoalteromonadaceae, Saprospiraceae and Haliangaceae). In Pratte et al. (2018) [15] study on the E. quadricolor/ A. clarkii system, Flavobacteriaceae appeared to be enriched in clown sh in physical contact with their anemone, one and two weeks after the start of the contact phase. However, the overall differentiation between physical contact and "remote interaction" control was just below or above the 5% threshold with no correction for multiple testing. This result suggests that differentiation between experimental groups was erased by the fact they shared the same water ow, thus compromising the comparison with control. Therefore, additional experiments in eld and laboratory conditions are needed to assess to which extent Flavobacteriaceae are involved in clown sh anemone microbiota convergence, and ultimately to establish their potential contribution in the mutual acceptance of both partners.
A complex network of inter-kingdom interactions One of our most salient results is the epithelium microbiota convergence between symbiotic partners, which could be linked to an important driver in the evolution of symbioses: the use of a shared signaling pathway as a common 'language' [38] between symbiotic partners. Furthermore, the persistence of the clown sh-anemone interaction signal in sh microbiota after symbiotic partner separation, and the potential direct or indirect link with the biochemical imprinting, needs further discussion. Overall, results from other vertebrate model species (e.g. mouse, human, zebra sh, stickleback) showed that there are multiple ways in which the dialogue between two eukaryotic hosts can involve microbial communities: host-microbiota interactions, microbiota-microbiota interactions, and host-host interactions [reviewed in 39,40].
(1) Host-microbiota interactions. A complex bidirectional communication is taking place between microorganisms and their host via chemical signals. The interaction between the symbiotic microbiota and the eukaryotic host brain was extensively documented [41][42][43][44][45]. For instance, a few studies have shown that prokaryote responsiveness to eukaryotic catecholamine hormones is widespread, and bacteria associated with animal surface epithelia are especially stress hormone responsive [45][46][47][48][49]. Thus, the catecholamines potentially produced by the clown sh brain during the initial phase of mutual acceptance could affect the taxonomic and functional pro le of the epithelial microbiome of both eukaryotic partners (clown sh and anemone). In a similar manner, symbiotic bacteria can also in uence host physiology. In our experiment, we observed that three ASVs related to Cellulophaga tyrosinoxydans, a tyrosinase producer, were especially associated to the convergence of microbiomes during the interaction period. These taxa might play signi cant roles in the chemical signaling convergence. For instance, melanin synthetized by bacterial tyrosinases are immunologically active compounds, known to bind diverse chemicals, and have many pharmaceutical applications including host skin protection against radiation, antioxidants, antiviral agents, or immunogens [49]. The metabolites repertory of the three ASVs related to C. tyrosinoxydans merits further investigations to highlight its functional importance in clown sh-anemone mutual acceptance.
(2) Microbiota-microbiota interactions. An important communication system between bacterial communities relies on quorum sensing, which regulates a wide variety of bacterial population densitydependent processes [50]. An increasing population density results in increased concentrations of bacterial autoinducers released by bacterial cells, which then elicit speci c responses. In our experiment, it is possible that the proximity of clown sh and anemones in physical and remote interaction groups induced the detection of an autoinducer shared in the epithelial microbiome of both eukaryotic host species. Then, density-dependent quorum sensing processes could have modulated the phylogenetic structure of bacterial communities from both hosts, potentially leading to their convergence (as observed in this study).
(3) Host-host interactions. Numerous studies have shown that symbiotic interactions between eukaryotic hosts can perturb the host-associated microbial habitats, which translates in a restructuring of the microbial communities involved [51]. For instance, in epithelial surface membranes, the initial colonization of a host by a parasite markedly alters the host's epithelial barrier by affecting mucus production and composition, tight junctions, and epithelial cell turnover [39]. Analogously to parasite colonization, the colonization of a sea anemone by a clown sh may modify the physical characteristics of the epithelial mucus habitat, and thus indirectly perturb the microbiomes of both symbiotic partners.
However, as observed with the remote interaction (i.e. no physical contact), our study suggests that most of the epithelial microbiota restructuring relies on remote biochemical communication.

Conclusions
Whether the restructuration of the microbiome of both hosts when in physical/remote interaction is a prerequisite for mutual acceptance (i.e. the "camou age hypothesis" [5]), or whether this community remodeling is an indirect consequence of an unknown remote biochemical detection system between partners is yet to be completely resolved. Here, we provide salient insights supporting the multilayered model of microbiome structuring from Shapira (2016) [52], which proposes that the variable environmentally modulated exible microbial pool, termed as transient microbiota, allows adaptation of holobionts to changing environments at the scale of microbial generation time. We propose that a remote interaction between both symbiotic partners triggers a gradual convergence of the transient epithelial microbiome of both partners, which may participate in establishing mutual acceptance. As such, our results suggest that the protective mechanism of A. percula clown sh, which has a narrow spectra of host sea anemone species, not only involves the coating of sh skin with components of the host anemone mucus [4,5], but also implicates a chemical dialog between symbiotic partners prior to the rst physical contact. Then, the present work provides additional insights regarding the hypothesis that epithelial mucus immunological pro le of clown sh changes to mimic that of the anemone prior to rst contact [6][7][8]: the convergence of the transient epithelial microbiota of both partners in remote interaction reveals that the chemical dialog might also trigger a restructuring of the sea anemone epithelial mucus.
As such, our results may challenge the traditional unidirectional chemical camou age hypothesis. In other words, both symbiotic partners each take a step towards the other to establish their mutualistic relationship. We hope that our study serves as a foundation for the design of other mechanistic studies to unravel this complex inter-kingdom interaction network between clown sh, anemones, and their bacterial symbionts.   a) Acclimation, three weeks until T0. b) Interaction, two weeks (from T1 to T3). c) Resilience, two weeks (from T4 to T5).