Clownfish and sea anemone rearing
Eighteen H. magnifica and 18 captive bred naïve A. percula juveniles from a tropical fish 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 flow to avoid any chemical contact prior to the experiment. H. magnifica 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. Clownfish 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 filter, and containing six biological replicates (i.e. six tanks). Two control groups: anemone control (AC), fish control (FC), and two test groups: physical interaction (i.e. fish and anemone in the same tank, PI), and remote interaction (i.e. fish and anemone in different tanks, all being connected to the same water flow, RI) (Fig. 1). To minimize bacterioplankton taxonomic drift due to independent water flow 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 clownfish control groups, we did not add the water mix (control anemone and control fish) in the control tanks. Following the acclimation period (Fig. 1a), the interaction period between clownfish and anemone for physical and remote interaction groups lasted two weeks (Fig. 1b), after which clownfish and anemones were separated for a two weeks resilience period (Fig. 1c). To mitigate possible effect of fish transfer from fish control to physical and remote interaction units, all the fish specimens were manipulated, including the six control fish 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 first 15 min of physical interaction between clownfish and anemone (PI), 15 min after transfer of fish 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 fish-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 fishes was immediately sampled by gently rubbing a sterile cotton swab on ≈50% of the total surface (upper half) of the right side of each fish outside of the water as in [23]. The same area was sampled on each fish 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 filtered 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 clownfish 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 amplified in a two-step dual-indexed PCR approach specifically designed for Illumina instruments by the Plateforme d’Analyses Génomiques (IBIS, Université Laval, Quebec city, Canada). The first PCR was performed with 16S region specific 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 flowcell and library specific barcodes). Please note that primers used in this work contain Illumina specific 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 amplification:
Forward specific primer (PCR #1):
ACACTCTTTCCCTACACGACGCTCTTCCGATCT-(347F)GGAGGCAGCAGTRRGGAAT,
Reverse specific primer (PCR #1):
GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT (803R)CTACCRGGGTATCTAATCC,
Forward specific primer (PCR #2):
AATGATACGGCGACCACCGAGATCTACAC[index1]ACACTCTTTCCCTACACGAC
Reverse generic primer (PCR #2):
CAAGCAGAAGACGGCATACGAGAT[index2]GTGACTGGAGTTCAGACGTGT.
Polymerase Chain Reactions (PCRs) were conducted with Q5 High-Fidelity DNA Polymerase from New-England Biolabs, in 25 uL reactions. The PCR program was: (a) 30 s 98°C; (b) 10 s 98°C; (c) 30 s 64°C; (d) 20 s 72°C; (e) 2 min at 72°C; 35 amplification cycles total. Amplified DNA was purified according to the manufacturer's instructions with AMPure beads (Beckman Coulter Genomics) to eliminate primers, dimers, proteins and phenols. Amplicon libraries were then sequenced on Illumina MiSeq (San Diego, CA, USA), including control samples.
Reads de-noising and ASV identification
Bioinformatics’ processing was undertaken using dada2 as reported in [24]. 4,770,388 raw reads were quality filtered 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 profiling 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 five 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, q-values).
Alpha diversity metrics
Simpson index was calculated using the R “vegan” package (v2.5.6). For the Faith PD index, the sequences were first 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 finally 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 significant 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: first, 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 influence 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 profile changes. LOWESS is a non-parametric and computationally demanding but robust and flexible regression method that uses locally weighted polynomials to fit a smoothed curve to scatterplot data [33]. In addition to smoothed curves, a confidence 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 significant 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-to-one (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 significant (Supplementary data), therefore confirming 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 first 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 clownfish 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 fish-anemone epithelial microbiota convergence. Then, de novo ASV abundances of clownfish and anemones were combined, PI and RI groups were combined as an interaction group, and clownfish and anemone controls were combined as a control group. ASVs with log2-normalized fold-change 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 identified during the first 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 clownfish), 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.