The role of host promiscuity in the invasion process of a seaweed holobiont

Sample collection

Algae were sampled from August 27^th to September 21^st (2017) from seven populations also collected for Bonthond et al. [28], including three native populations; Akkeshi (Japan), Soukanzan (Japan), Rongcheng (China); and four non-native populations; Pleudihen-sur-Rance (France), Nordstrand (Germany), Cape Charles Beach (Viriginia) and Tomales Bay (California, Fig. 1, Table S1). Individuals fixed to hard substratum (see [30]) were sampled at least a meter apart from one another and stored in separate plastic bags. As A. vermiculophyllum has a complex, haplodiplontic life-cycle only diploids were included in the experiment. Life-cycle stages were identified in the field with a dissecting microscope or post-hoc by microsatellite genotyping [31]. After transport in coolers and storage at 4 °C in the lab, bags with algae were shipped to Germany, arriving within 4–6 days after collection. In the climate room (15 °C), individuals were transferred to separate transparent aquaria with transparent lids, containing 1.75 L artificial seawater (ASW) prepared from tap water and 24 gL⁻¹ artificial sea salt without CaCO₃ (high CaCO₃ concentrations increase disease risk, Weinberger data unpublished) and exposed to 12 h of light per day (86.0 µmol m⁻²s⁻¹ at the water surface). Aquaria were moderately aerated with aeration stones. Per population, four diploid individuals were acclimated over 31–32 days to climate room conditions prior to starting the experiment. Water was exchanged weekly with new ASW enriched with 2 mL Provasoli-Enrichment Solution (PES; [32]). At the start of the experiment, wet weight was recorded and individuals were divided into two parts of ~10 g each and placed into two plastic tanks with 1.75 L water and 2 mL PES (Fig. 1).

Experimental setup

To rigorously disturb the microbial community, one of each of the pairs of aquaria containing the same algal individual was treated with a combination of antibiotics, aiming to increase the effectivity (10 mgL⁻¹ ampicillin, 10 mgL⁻¹ streptomycin, 10 mgL⁻¹ chloramphenicol) and the other (control) remained untreated. All experimental work was conducted with disposable gloves and sterilized equipment, to minimize contamination. After three days, the water was removed from all tanks (treated and control) and the wet weight was recorded for all algae. All individuals were rinsed with one 1.75 L volume ASW and re-incubated in 1.75 L ASW. Subsequently, both groups received new ASW with 2 mL PES weekly and individuals treated with antibiotics received also 2 mL inoculum. The inoculum was prepared from individuals of all 7 populations, following the procedure to remove epibiota as described in Bonthond et al. [28]. Briefly, apical fragments of 1 g were separated from the thallus and transferred to 50 mL tubes containing 15 ± 1 glass beads (3 mm) and 15 mL ASW and vortexed for 6 min to separate epibiota from the algal tissue. In total, 8 samples were prepared from one individual per population. The resulting suspensions were pooled and mixed with glycerol (20% final glycerol concentration), aliquoted in 50 mL tubes and stored at −20 °C. For each water exchange, a new aliquot was defrosted at room temperature and added to the water of treated algae. Wet weight was recorded weekly with water exchanges. Before weighing the individual on aluminum foil, it was dipped twice on a separate aluminum foil sheet, to reduce attached water in a systematic way. Endo- and epiphytic microbiota were sampled in the field (t_field, [28]), at the start of the experiment (t₀), after one week (t₁), two weeks (t₂), four weeks (t₄) and six weeks (t_6, Fig. 1). To equalize acclimation times across populations the experiment was stacked into five groups (Table S2). At each sampling moment, 0.5 or 1 g of tissue was separated from all individuals with sterilized forceps and epibiota were extracted similarly to the preparation of the inoculum. The resulting suspension was filtered through 0.2 µm pore size PCTA filters. Both the filters and the remaining tissue were preserved at −20 °C.

DNA extraction and amplicon sequencing

Tissue samples were defrosted, rinsed with absolute ethanol and DNA free water to remove hydro- and moderately lipophilic cells and molecules from the surface and cut to fragments with sterilized scissors. DNA was then extracted from these fragments (endobiota) and from preserved filters (epibiota) using the ZYMO Fecal/soil microbe kit (D6102; ZYMO-Research, Irvine, CA, USA), following the manufacturer’s protocol. Although this method to separate endo- and epibiota was shown to resolve distinct communities [28], tightly attached epiphytic cells may not be completely removed from the surface and detectable in endophytic samples as well. Two 16S-V4 amplicon libraries, over which the samples were divided in a balanced manner, were prepared as in Bonthond et al. [28], following the two-step PCR strategy from Gohl et al. [33], using the same set of 16S-V4 target primers and indexing primers. The libraries were sequenced on the Illumina MiSeq platform (2×300 PE) at the Max-Planck-Institute for Evolutionary Biology (Plön, Germany), including four negative DNA extraction controls and four negative and positive PCR controls (mock communities; ZYMO-D6311). The fastq files were de-multiplexed (0 mismatches). Relevant field samples from Bonthond et al. [28] were combined with the new dataset and assembled, quality filtered and classified altogether with Mothur v1.43.0 [34] using the SILVA-alignment release 132 [35]. Sequences were clustered within 3% dissimilarity into OTUs using the opticlust algorithm. Mitochondrial, chloroplast, eukaryotic and unclassified sequences were removed. To prepare the community matrix we discarded singleton OTUs (in the full dataset), samples with <1000 read counts and OTUs from which all sequences were removed after the previous step. De-multiplexed reads and corresponding metadata were deposited in the SRA database (accession: PRJNA612003).

Functional profiling

To predict functional groups we used Picrust2 software [36] with default settings. Using KO-numbers from the Kyoto encyclopedia of Genes and Genomes KEGG [37], we defined the following functional groups: autotrophy (RuBisCo; K01601), aerobic heterotrophy (COXIII; K02276), anaerobic heterotrophy (adenylyl-sulfate reductase; K00394, methane/ammonia monooxygenase; K10944 and fumarate reductase; K00244 combined) and diazotrophy (nifH; K02588).

Identification of core microbiota

Geography-independent OTUs identified in Bonthond et al. [28] were re-identified in the new dataset by cross-comparing all OTU sequences with those of the core microbiota from the previous study [28]. Sequences that were identical, or most similar, were reclassified to epiphytic, endophytic or algal core OTUs.

Statistical modeling

We used the relative growth rate (RGR) as a measure of performance. As tissue was removed from all individuals at each timepoint (t₀–t₆) RGR was obtained by dividing the gained wet weight by the weight after sampling at the previous timepoint and divided by the number of days in between sampling points (thus expressed in % growth d⁻¹). RGR was analyzed with a linear mixed model as a function of range, treatment, time (weeks) and all interaction terms. To model the relationship between RGR and time in a flexible way time was included as a third order polynomial, therewith considering temporal variation, without imposing shapes more complex than a third order polynomial. We included population–identity and individual-identity as random intercepts to represent the genetic population structure by A. vermiculophyllum at the local scale [24] and to account for non-independence within individuals.

For alpha-diversity, we used OTU-richness rarefied to 1000 reads per sample and the probability of interspecific encounter (PIE) as a measure of evenness, obtained with the package mobr [38]. To compare diversity in the field with the beginning of the experiment, field and control samples from the first timepoint were used, including substrate and time as fixed and population and individual-identity as random effects. Then, we fitted third order polynomial functions of time on the subset of the data including experiment samples (t_0–6). These models also included the predictors substrate, treatment and random intercepts population- and individual-identity. To meet normality, PIE was logit transformed. To account for possible effects resulting from differences in read counts across samples, the log of the sequencing depth (LSD) was included in all models as a continuous variable. Predicted functional groups were analyzed with the same model structures. To meet normality, responses were log (+1 when including zeros) or squared-root transformed.

To analyze community composition between treatments we used multivariate generalized linear models (mGLMs) from the R package mvabund [39] in a two-step approach. The community matrices were trimmed to the 95% most abundant OTUs and split by substrate to analyze epi- and endobiota separately. First, a mGLM was used to remove the effects of sequencing depth (by including LSD) and differences among populations (by including population-identity). The mGLMs assumed a negative binomial distribution with a log in the link function. Second, on the residuals we ran a mGLM in response to treatment, time (a third order polynomial) and the interaction. This model assumed a Gaussian distribution as the residuals from the first model were normally distributed in the link function. Multivariate statistics were obtained by resampling the univariate models with 500 bootstrap iterations. Compositional differences between control and treatment over time were visualized using non-metric multidimensional scaling (nMDS) on the rescaled residuals of the first model. Group centroids and corresponding 95% confidence regions were computed with the R package vegan [40]. Compositional changes at the univariate level (i.e., specific OTUs) were visualized with a heatmap including the fifty most abundant OTUs by treatment and timepoint.

Beta-diversity was analyzed with pairwise Bray–Curtis distances from the epi- and endophytic datasets that were adjusted for the sequencing depth using mGLMs as a function of LSD. Additionally, we ran these models on weighted UniFrac distances, for which representative sequences of all OTUs were aligned with MAFFT v7.221 [41], with Saccharomyces cerevisiae as outgroup, and clustered into a maximum-likelihood tree with RAxML v8.2.12 [42] with the GTR + G substitution model and a 1000 bootstrap iterations. We compared distances among individuals, calculated between samples within the same timepoint (i.e., t₀–t₀, t₁–t₁, etc.) and regressed those against a third order polynomial of time, a new factor population (levels: within- and between-populations), range (native and non-native), treatment (control and treated) and all interactions. The random intercepts population-combination and individual-combination were included to account for non-independence resulting from calculating distances by making different combinations with the same individuals. To characterize how community composition changed with respect to the composition observed in the field, we calculated Bray–Curtis and weighted UniFrac distances between individuals in the experiment (t₀–t₆) and individuals in the field (t_field) within the same population. These models (for endo- and epibiota separately) included a third order polynomial of time, the range, their interaction and the random intercept individual-combination.

Proportional changes in core OTU abundances were analyzed with a mixed linear model using the subset of our data containing only field samples and the final timepoints (t₆) and included the variables treatment, substrate, range and the interactions. Population- and individual-identity were included as random intercepts. All univariate analyses were conducted using the R package lme4 [43], calculating marginal and conditional R² values (variation explained by fixed effects and fixed plus random effects, respectively) with the r.squaredGLMM function [44]. Violations of model assumptions were verified visually with QQ-plots and residual-vs-fitted-plots for univariate and multivariate analyses.

Source: Ecology - nature.com