Exploring mechanisms that affect coral cooperation: symbiont transmission mode, cell density and community composition

Coral collection and acclimation

Fragments from 60 unique coral colonies, ∼20 cm in diameter, were collected from reefs on the Central GBR from the 8–22 April 2015 under the Great Barrier Reef Marine Park Authority permits G12/35236.1 and G14/37318.1, prioritizing collection of focal species pairs by transmission mode from the same reef environment, including location and depth. A total of 10 corals of each species, Galaxea acrhelia and Galaxea astreata were collected from Davies Reef (18°49.816′, 147°37.888′, 11 April 2015), 10 corals each of A. millepora and M. aequituberculata were collected from Pelorus Island (18°33.358′, 146°30.276′, 18 April 2015) and 10 corals each of Goniopora columna and P. lobata were collected from Pandora Reef (18°48.778′, 146°25.593′, 21 April 2015) from depths of <10 m across all reefs, and from the same depth within each site. Corals were transported to the Australian Institute of Marine Science Sea Simulator facility and placed in shaded outdoor holding tanks with 0.2 μM filtered flow-through seawater (FSW, 27 °C, 150 μmol quanta m²/s). Each individual coral colony was further cut into six replicates using a diamond blade band saw and these fragments were mounted on aragonite plugs using either super-glue or marine epoxy. On 1 May 2015, fragments were moved into indoor experimental aquaria, which consisted of 12 50-L treatment tanks fitted with 3.5 W Turbelle nanostream 6,015 pumps (Tunze, Penzberg, Germany) with filtered flow-through seawater (FSW, ∼25 l/h) at 27 °C with Hydra 52 lights (Aqua Illumination, C2 Development, Inc., Ames, IA, USA) on a 12:12 light/dark cycle mimicking natural irradiance patterns, peaking at 130–160 μmol quanta m²/s at ‘midday’, which actually occurred at 10:30 GMT +10. Tanks, coral racks and plugs were cleaned daily. Beginning on 7 May, at 20 mins before dark (16:30 GMT +10) every day, corals were fed Artemia naupli at a density of 1–1.5 naupli/5 ml and rotifers at a density of 1–3/ml.

Experimental conditions and physiological trait measurements

Beginning on 28 May 2015, effective quantum yield (EQY) of Symbiodinium photosystem II was measured daily for all experimental fragments using a pulse amplitude modulated fluorometer (diving-PAM, Walz) fitted with a plastic fibre optic cable (Fig. S1). Measurements were made using factory settings with a measuring intensity of 12 and a gain of 5 and taken at peak light intensity (between 09:30 and 11:30 GMT +10). EQY values were used to guide the timing of the final sampling time point (Fig. S2), where a decline reflects an impact on the photosynthetic condition of the Symbiodinium (Ralph, Larkum & Kuhl, 2005) because we aimed to target the onset of the coral bleaching response rather than the end-point.

On 31 May 2015, temperatures in the heat treatment tanks were increased at a rate of 1 °C per day until temperatures reached 31 °C (day 4, Fig. S3). Sample time points occurred on day 2 (29 °C, 1 June), day 4 (31 °C, 3 June) and day 17 (31 °C, 16 June, Fig. S3). On each sampling day, one replicate fragment of each colony (n = 10) and species (n = 6) from each temperature treatment (n = 2; n = 120 total per sampling day) were used to measure net photosynthesis following the two-point method originally described and validated for A. millepora in Strahl et al. (2015). It is important to note that this method was not validated for additional species or under the experimental conditions used in the present study, nor were species specific photosynthesis-irradiance curves quantified to determine an appropriate saturating irradiance. Briefly, corals were incubated in enclosed 600 ml acrylic chambers at their respective treatment temperatures and light levels for 1.5 h. Chambers were placed onto custom built tables with rotating magnets, which served to power stir bars within each chamber to facilitate water mixing. For each run, four chambers without corals were used as blanks to account for potential changes in oxygen content due to the metabolic activity of other microorganisms in the seawater. For net photosynthesis measures, the O₂ concentration of the seawater in each chamber was measured at the end of the run using a hand-held dissolved oxygen meter (HQ30d, equipped with LDO101 IntelliCAL oxygen probe; Hach, Loveland, CO, USA). Values from blank chambers were subtracted from measures made in coral chambers and the subsequent rate of net photosynthesis was related to coral surface area, calculated in μg O₂/cm²/min.

To measure the fraction of autotrophically derived carbon translocated to host animals, five colony fragments of each species from each treatment (n = 30 control, n = 30 heat) were placed into 18-L of FSW with a 5-W aquarium pump for circulation and ¹⁴C-bicarbonate (specific activity: 56 mCi/mmol) was added to a final concentration of 0.28 μCi/ml. Corals were incubated for 5 h in experimental tanks which experienced the normal experimental irradiance profile from 09:00 to 14:00 GMT +10, rinsed with flow-through FSW for 1 h to remove remaining unfixed ¹⁴C then snap frozen in liquid nitrogen.

Tissue was removed from snap frozen coral skeletons using an air gun and homogenized for 60 s using a Pro250 homogenizer (Perth Scientific Equipment, Perth, WA, Australia). A 300 μl aliquot of the tissue homogenate was fixed with 5% formalin in FSW and used to quantify Symbiodinium cell density. The average cell number was obtained from four replicate haemocytometer counts of a 1 mm³ area and cell density was related to host protein content (as assessed below) and expressed as cells/mg host protein. Although surface area is commonly used as the normalization factor for Symbiodinium cell density, it has been recognized that areal abundance does not account for differences in host biomass (reviewed in Cunning & Baker (2014)). We found that normalization to soluble host protein more accurately reflected this difference in biomass, given that tissue thickness is significantly greater in Goniopora columna and not adequately accounted for by skeletal surface area normalization alone (Fig. S4). We determined the degree of bleaching, or the reduction in symbiont density in response to sustained thermal stress, by calculating the difference in symbiont cell densities between heat-treated and paired control samples following the full 17 days of experimental treatment and used this value as a proxy of holobiont fitness.

An additional 1 ml aliquot of holobiont homogenate was frozen at −20 °C. The remaining homogenate was centrifuged for 2 min at 3,500 rcf to separate host and symbiont fractions and two ml of the host tissue slurry was frozen at −20 °C. Total soluble protein was quantified for host tissue samples in duplicate using a colorimetric assay (Bio-Rad Protein Assay Kit II; Bio-Rad Laboratories, Inc., Hercules, California, USA) following the manufacturer’s instructions. Coral skeletons were rinsed with 5% bleach then dried at room temperature. Skeletal surface area was quantified using the single wax dipping method (Veal et al., 2010) and skeletal volume (used to standardize respirometry chamber volumes) was determined by calculating water displacement in a graduated cylinder.

Sample radioactivity was determined using a liquid scintillation counter (Tri-Carb 2810TR v2.12; Perkin Elmer, Waltham, MA, USA). Triplicate host and duplicate holobiont 300 μl tissue homogenate aliquots were mixed with 3.5 ml Ultima Gold XR liquid scintillation cocktail (Perkin Elmer, Waltham, MA, USA). Samples were temperature and light adapted for 1 h and then counted for 1.5 min using the default parameters. Counts per minute were converted to disintegrations per minute (DPM) using a standard curve derived from a ¹⁴C Ultima Gold Quench Standards Assay (Perkin Elmer, Waltham, MA, USA). Technical replicates were averaged. Host DPM values were divided by holobiont DPM values to yield the fraction of autotrophically derived carbon shared between partners, which we defined as our metric of host-symbiont cooperation. The degree of symbiont parasitism was calculated as the difference in this proportion of photosynthetically fixed carbon translocated to hosts between heat-treated and paired control samples following the full 17 days of experimental treatment, sensu Baker et al. (2018).

To identify major Symbiodinium clades hosted by focal species, DNA was extracted from symbiont fractions of the tissue homogenate for each replicate colony of each species (n = 60, all control tank fragments) using Wayne’s method (Wilson et al., 2002). A restriction digest of the large subunit (LSU) region of Symbiodinium rRNA (Baker & Rowan, 1997; Palstra, 2000) consistently revealed single bands indicating the dominance of a single Symbiodinium clade for four of the six species (A. millepora, M. aequituberculata, Goniopora columna and P. lobata) whereas communities in Galaxea astreata and Galaxea acrhelia appeared more variable (Fig. S5). Therefore, DNA was pooled in equal proportions by species (n = 6 samples) to identify general species-specific communities using amplicon sequencing of the ITS2 region of Symbiodinium rRNA. Additional DNA samples for each of the Galaxea astreata and Galaxea acrhelia colony replicates (n = 20 samples) were also submitted for sequencing at the Genome Sequencing and Analysis Facility at the University of Texas at Austin. Given the high diversity subsequently observed in pooled samples of P. lobata and M. aequituberculata, additional colony replicate samples of each (n = 3 and n = 5, respectively) were later submitted for sequencing at Oregon State University’s Center for Genome Research and Biocomputing to compare symbiont community composition from individual samples.

Amplicon sequencing analysis

For sample libraries prepared and sequenced at UT Austin’s GSAF, the ITS2 primers of Pochon et al. (2001) were used. For colony replicate samples sequenced at OSU’s CGRB (individual P. lobata and M. aequituberculata samples), the ITS2 primers of LaJeunesse (2002) were used to prepare libraries. All sequencing libraries were subsequently analysed together. Prior to analysis, raw read data was filtered to remove reads which contained Illumina sequencing adapters or did not begin with the correct ITS2 amplicon primer sequence using the BBMAP package ver. 37.75 (http://sourceforge.net/projects/bbmap/). The DADA2 pipeline (Callahan et al., 2016) implemented in R Development Core Team (2017) was then used to infer sequence variants. Read data was analysed according to the following tutorial (https://benjjneb.github.io/dada2/tutorial.html) with the following modifications: ITS2 primers were trimmed from the beginning of each read and forward and reverse reads were truncated at 210 and 160 bp, respectively, to remove low quality bases at the end of reads. Additional paired end reads were discarded if they exhibited more than one expected error or when a quality score of 2 or less was encountered. In addition, when inferring sequencing variants using the dada() command, the BAND_SIZE flag was set to 32 as is recommended for ITS data (https://benjjneb.github.io/dada2/tutorial.html). The distribution of raw read data per sample and reads lost during quality filtering and processing steps can be found in Table S1. Post-clustering curation of identified amplicon sequence variants (ASVs) was accomplished with the LULU algorithm, using the default parameters (Frøslev et al., 2017). The MCMC.OTU package (Green et al., 2014) was then used to remove sample outliers with low counts overall (z-score < −2.5) and remaining ASVs of abundance less than 0.001 (Quigley et al., 2014) prior to statistical analysis. In total, 12 ASVs were identified across samples that satisfied these criteria. These high confidence sequence variants were taxonomically classified through a blast search against the GeoSymbio ITS2 database https://sites.google.com/site/geosymbio/downloads (Franklin et al., 2012), and the best match was recorded. In cases where variants matched equally well to multiple references, all top hits were reported (Table S2). Resequencing of the individual P. lobata and M. aequituberculata samples indicated that the P. lobata pooled sample was contaminated at some stage of the sequencing process, as individual samples did not match the host pool (Fig. S6). In addition, two Galaxea samples (Gacr1 and Gast4) appeared mislabeled, based on the pattern of symbiont community differences among species (Fig. S6) and we removed these samples prior to further analysis. To account for differences between individual samples and species pools, we calculated the normalized mean abundance of each variant across individual samples within a species and replaced the pooled sequence sample with this in silico pooled value, when available.

Statistical analyses

All statistical analyses were performed in R version 3.4.2 (R Development Core Team, 2017). A series of linear mixed models were used to determine the effect of species, temperature treatment and sampling day on physiological metrics. We used a conservative approach to evaluate the effect of transmission mode on physiological metrics. Rather than modeling fixed effects of transmission mode directly, we modeled the fixed effects of species (levels: Amil, Maeq, Gcol, Plob, Gast, Gacr), sampling day (levels: day 2, day 4, day 17), temperature treatment (levels: heat, control) and all possible interactions on Symbiodinium cell density, net photosynthesis, and host-symbiont resource sharing using the lme command of the nlme package (Pinheiro et al., 2013), including source coral colony identity nested within reef of origin as a scalar random factor. All models were assessed for normality of residuals and homoscedasticity. Symbiont densities required a log-transformation to meet normality assumptions. To assess goodness-of-fit, we used the function rsquared.glmm to calculate the conditional R² value for each of our mixed models (Johnson, 2014). Significance of fixed factors within models was evaluated using Wald tests and Tukey’s post hoc tests were used to evaluate significance among levels within factors and interactions when warranted. In these models, a significant effect of transmission mode would have been detected first as a significant difference among species, but significant differences in the hypothesized direction between each pair of vertical and HTs in the subsequent Tukey’s tests were also required to satisfy reporting an effect of transmission mode overall.

To determine the impact of specific Symbiodinium community characteristics on responses among species, the DESeq package (Anders & Huber, 2010) was used to construct a series of generalized linear models to evaluate differences in in the abundance of each ASV by species with respect to aspects of host-symbiont cooperation: mean initial symbiont density (Day 2, heat samples); the mean fraction of carbon shared by symbionts under ambient conditions (Day 2, 4, 17 control samples), the mean difference in the fraction of carbon shared by symbionts during bleaching (Day17, heat vs. control samples) and the mean difference in symbiont density during bleaching (Day17, heat vs. control samples). Models were run for 30 iterations and for models that did not converge, p-values were converted to NAs (the standard notation for missing data), prior to applying a multiple test correction (Benjamini & Hochberg, 1995).

Predictive relationships between Symbiodinium cell density, the degree of bleaching, cooperation, and symbiont parasitism, in addition to Symbiodinium community diversity (quantified using the inverse Simpson index, (Magurran, 2004)) and composition (the proportion of Clade D) were explored at the species level using a series of linear models. When a significant difference was detected among species for physiological trait comparisons, the model was re-run within each species and a multiple test correction using the method of Benjamini & Hochberg (1995) was applied. Analyses of symbiont community composition were run both using the full six species averages and excluding the species for which individual sample sequencing replicates were not available (A. millepora and Goniopora columna). As models did not substantially differ, we report results for the full six species.

Physiological metrics of thermal tolerance and cooperation

Significant fixed effects of species, temperature treatment, sampling time and the time*treatment interaction were detected for the Symbiodinium cell density model (R_GLMM² = 0.53, Table S3). The fixed difference among species was driven by the low symbiont density on average in P. lobata, which differed significantly from densities in A. millepora, M. aequituberculata, Galaxea astreata and Galaxea acrhelia (Tukey’s HSD < 0.05, Fig. 1). No differences were detected in post hoc tests between focal species pairs by symbiont transmission mode. Across species, cell densities did not differ between control and heat-treated corals on days 2 and 4, but were reduced in heat treated corals on day 17, by 5.4 × 10⁵ cells/mg host protein on average (Tukey’s HSD < 0.001, Fig. 1).

The model for changes in the rate of net photosynthesis identified significant effects of species, sampling time and the temperature treatment*time interaction (R_GLMM² = 0.53, Table S4). Average net O₂ production rate was highest in M. aequituberculata, differing significantly from A. millepora, Goniopora columna, Galaxea acrhelia and P. lobata Tukey’s HSD < 0.05, Fig. 1), and lowest in Goniopora columna, significantly more so in comparison to M. aequituberculata, P. lobata, Galaxea astreata and Galaxea acrhelia (Tukey’s HSD < 0.001, Fig. 1). Net O₂ production rate was also lower on average in A. millepora than in M. aequituberculata and Galaxea astreata (Tukey’s HSD < 0.01, Fig. 1). Considering differences between focal species pairs by symbiont transmission mode, the rate of net photosynthesis was significantly higher in the vertically transmitting M. aequituberculata and P. lobata than in their horizontally transmitting counterparts, A. millepora and Goniopora columna. However, it is important to note that absolute rates of oxygen production among focal species are much lower than previously reported values for corals (Anthony et al., 2008). Corals were acclimated to a common light environment, but data on species and time-point specific photosynthesis-irradiance curves are lacking. Consequently, the differences among species may be influenced by variation in species-specific photobiology while differences within species may be influenced by shifting photophysiological performance over the experimental time-course.

Over the course of the experiment, net photosynthesis rates were elevated in heat-treated corals on days 2 and 4 relative to their respective controls (day 2: 0.09 μg O₂/cm²/min; day 4: 0.06 μg O₂/cm²/min; Tukey’s HSD < 0.01, Fig. 1) and reduced in heat-treated corals on day 17, by 0.1 μg O₂/cm²/min (Tukey’s HSD < 0.001, Fig. 1).

The proportion of carbon photosynthetically fixed by symbionts then translocated to hosts was significantly different among species, temperature treatment, sampling time, the treatment*time interaction and the species*treatment interaction (R_GLMM² = 0.73, Table S5). Carbon sharing in ambient conditions was significantly lower in Goniopora columna in comparison to all other species (Tukey’s HSD < 0.05, Fig. 1) and highest in the vertically transmitting Galaxea acrhelia and P. lobata, significantly more so than in M. aequituberculata and their horizontally transmitting counterparts, Galaxea astreata and Goniopora columna (Tukey’s HSD < 0.01, Fig. 1). Over time, no significant differences were detected in carbon sharing between control and heat-treated corals on days 2 and 4, but proportional translocation was significantly lower in heat-treated corals on day 17 (Tukey’s HSD < 0.001, Fig. 1). Relative to controls, carbon sharing decreased slightly under heat treatment on average in all species save Goniopora columna, though the only significant difference was observed in Galaxea acrhelia (heat vs. control Tukey’s HSD < 0.01, Fig. 1). Although the species*treatment*time interaction was not significant in the final model, this was likely driven by the decrease in heat-treated corals on day 17, which is most evident in A. millepora, M. aequituberculata, Galaxea astreata and Galaxea acrhelia.

Relationships between symbiont density, degree of bleaching and carbon translocation

Symbiont cell density did not explain a significant proportion of the variance in carbon translocation under ambient conditions (control corals on days 2, 4 and 17) due to differences among species (F_1,5 = 25.57, P < 0.001, Fig. 2). Within species, carbon translocation increased with increasing symbiont density in Galaxea acrhelia (R² = 0.22) but this relationship became non-significant after applying a multiple test correction.

Symbiont cell density in heat-treated corals on day 2 predicted a small portion of the variance in bleaching intensity on day 17 (calculated as the difference in cell density between paired control and heat-treated coral fragments). Greater bleaching was associated with higher initial symbiont cell densities (β = −0.34, R² = 0.06, P = 0.04, Fig. 3A) and this relationship did not differ significantly among species. A similar, but non-significant trend was also observed between bleaching on day 17 and symbiont cell density on day 4 (β = −0.24, R² = 0.04, P = 0.07).

No significant relationship was detected between the degree of bleaching and the degree of parasitism (calculated as the difference in carbon translocation between heat and control treated samples of paired corals by source colony) at the end of the experiment (Fig. 3B). Nor was the degree of parasitism explained by differences in symbiont cell density of heat-treated corals on day 2, but relationships differed among species (F_1,5 = 6.88, P < 0.001). In P. lobata, a higher symbiont density on day 2 was associated with lower parasitism, or an increase in proportional carbon translocation in heat-treated corals on day 17, and this relationship remained marginally significant even after applying a multiple test correction (R² = 0.89, P = 0.06, Fig. 3C).

The role of Symbiodinium community composition

Across all individual and pooled samples, 12 ASVs of sufficient representation (>0.1%) abundance, sensu (Quigley et al., 2014) were identified and consisted of clade C and D-type Symbiodinium (Table S2 and Fig. S6). Symbiodinium community composition varied among species (Fig. 4) but no relationships were detected between the abundance of individual ASVs and symbiont density, degree of bleaching or carbon translocation. Nor did we detect any relationship between Symbiodinium community diversity overall and symbiont density on day 2, mean carbon translocation under ambient conditions or the degree of bleaching at the end of the experiment (Figs. 5A, 5B and 5D). A marginally significant negative relationship was detected between community diversity and the degree of parasitism at the end of the experiment, where a more diverse community was associated with a greater degree of symbiont parasitism (R² = 0.56, P = 0.054, Fig. 5C). We did not detect any relationships between the percent of Clade D in the symbiont community and symbiont density on day 2, mean carbon translocation under ambient conditions, or the degree of parasitism or bleaching at the end of the experiment (Figs. 5E–5H), but this was likely due to the fact that the abundance of D was very low in all species except Galaxea acrhelia, however, did exhibit the highest carbon translocation rates under ambient conditions and the greatest transition towards parasitism on day 17 (Figs. 5F and 5G).