Uncovering Spatio-Temporal and Treatment-Derived Differences in the Molecular Physiology of a Model Coral-Dinoflagellate Mutualism with Multivariate Statistical Approaches

Multivariate statistical approaches (MSA), such as principal components analysis and multidimensional scaling, seek to uncover meaningful patterns within datasets by considering multiple response variables in a concerted fashion. Although these techniques are readily used by ecologists to visualize and explain differences between study sites, they could theoretically be employed to differentiate organisms within an experimental framework while simultaneously identifying response variables that drive documented experimental differences. Therefore, MSA were used herein to attempt to understand the response of the common, Indo-Pacific reef coral Seriatopora hystrix to temperature changes using data from laboratory-based temperature challenge studies performed in Southern Taiwan. Gene expression and physiological data partitioned experimental specimens by time of sampling, treatment temperature, and site of origin upon employing MSA, signifying that S. hystrix and its dinoflagellate endosymbionts display physiological and molecular signatures that are characteristic of sampling time, site of colony origin, and/or temperature regime. These findings promote the utility of MSA for documenting biologically meaningful shifts in the physiological and/or sub-cellular response of marine invertebrates exposed to environmental change.


Introduction
In recent years, a concerted effort has been made to better understand the basic biology of reef-building corals [1][2][3][4], as well as their response to changing environments [5][6][7]; the latter topic is especially pertinent given the extent of the anthropogenic pressures currently facing the high biodiversity ecosystems constructed by these cnidarian-dinoflagellate (genus Symbiodinium) endosymbioses [8,9].The impact of changing environments on corals is undoubtedly complex, and many species have been shown to acclimate to abiotic regimes previously hypothesized to compromise the integrity of these calcium carbonate-accreting mutualisms.As an example, although most scleractinian coral-Symbiodinium associations readily dissociate when exposed to even small changes in their aquatic milieu (particularly with respect to temperature [10,11]), those of Southern Taiwan have proven to be markedly resilient to an array of laboratory-simulated environmental challenges (Table 1).Corals can acclimate to OA on a multi-month timescale.[22] For instance, the Indo-Pacific reef-builder Seriatopora hystrix showed no mRNA-level molecular chaperone response when exposed for 2 days to 30 • C [12], a temperature hypothesized to ultimately elicit bleaching in this species based on observations made in Japan [23] and elsewhere.Instead, the expression of only 2 genes out of the 14 targeted (6 from Symbiodinium and 8 from the coral host), the cytoskeleton genes β-actin (actb) and α-tubulin (tuba), were determined by real-time PCR (qPCR) to be affected by temperature [13].Mayfield et al. [12] hypothesized that such a lack of a molecular chaperone response, in particular, in either compartment of this holobiont ("host + symbiont") may have been due to mRNA "front-loading" (sensu [24]) in samples of this "S.hystrix short-term temperature experiment" (SHSTTE).Briefly, corals of Southern Taiwan inhabit environments characterized by episodic upwelling, whereby temperatures may change by up to 9 • C in a matter of hours [25].Therefore, they could be predicted to express high levels of heat shock proteins (HSPs) and other stress genes/proteins even during ambient conditions in order to have the molecular machinery requisite for a temperature change-induced stress response at the time temperatures begin to fluctuate due to upwelling.
As an unexplored counter-hypothesis, it is plausible that concerted, biologically meaningful changes in expression of multiple gene mRNAs and other molecular physiological response variables were simply overlooked due to having used univariate statistical approaches only.Multivariate statistical approaches (MSA), such as principal components analysis (PCA) and multidimensional scaling (MDS), can uncover treatment-derived and spatio-temporal differences not revealed by univariate statistics-based approaches employing standard ANOVA models by instead looking at the relationships or correlations between various combinations of response variables simultaneously.Specifically, MSA can differentiate samples and treatments by integrating data across multiple parameters and so can partition samples within an experimental dataspace in a holistic manner.Another advantage of MSA, such as multivariate ANOVA (MANOVA), is that such techniques are more statistically conservative when analyzing datasets featuring a large number of response variables; by assessing all parameters (e.g., 17 in the SHSTTE) in an integrated, single-step model, the chances of making a type I error are substantially reduced.
Given these merits, as well as the fact that MSA have been utilized extensively in recent years to gain insight into the ecophysiology of reef-building corals (e.g., [26][27][28][29][30][31]), MSA were used herein to ascertain whether the S. hystrix-Symbiodinium holobiont was truly unresponsive to a short-term exposure to a temperature treatment hypothesized to elicit stress (sensu [32]).As a comparison, the dataset of the "S.hystrix variable temperature study" (SHVTS), which was also conducted at Taiwan's National Museum of Marine Biology and Aquarium (NMMBA), was re-explored, as corals of this study showed clear physiological and gene expression differences across both temperature regimes (stable vs. variable) and site of origin [14].Regarding the latter factor, unlike the SHSTTE, in which all corals were from an upwelling site within Nanwan Bay (Taiwan's southernmost embayment), Houbihu, half of those corals of the SHVTS were from a non-upwelling site, Houwan, which abuts NMMBA and is characterized by low coral cover and poor water quality due to agricultural runoff [33].It was predicted that MSA could be used to conclusively demonstrate the lack of a gene expression effect on high temperature samples of the SHSTTE and, similarly, further verify both site of origin and temperature treatment differences in the molecular physiology of samples of the SHVTS.Furthermore, it was hypothesized that MSA could determine whether corals exposed to variable temperature in the SHVTS showed a greater diversity in their physiological response than those exposed to stable temperature; this phenomenon was hypothesized to occur based on the fact that only half of the exposed experimental nubbins in the SHVTS (those of Houbihu) actually experience variable temperature profiles in situ [14], whereas all colonies from which experimental nubbins were made routinely experience stable temperature regimes.

The Experiments
The SHSTTE [12,13] and SHVTS [7,[13][14][15] are described in prior works.Briefly, whole, unfragmented S. hystrix colonies from Houbihu were exposed to either the control (27 • C) or elevated temperature (30 • C) for 48 h in the former (Table 1), and RNAs, DNAs, and proteins were extracted from triplicate colonies housed within each of triplicate tanks at each of the two temperature treatments at each of four sampling times (6, 12, 24, and 48 h; 18 samples/sampling time).A tank average was calculated across the three pseudo-replicates within the same tank sampled at each time, resulting in 24 data points that were analyzed by the MSA discussed below (n = 3 biological replicates/sampling time/treatment × 4 sampling times × 2 temperature treatments).All coral colonies were collected under Kenting National Park permit 0992900398 issued to Dr. Tung-Yung Fan (2009).Univariate repeated measures ANOVAs were used previously [12,13] to assess the effects of time, treatment temperature, and their interaction on the 17 response variables described below.
In the SHVTS, six coral colonies were sampled from each of the two sites of origin (Houbihu [upwelling site] and Houwan [non-upwelling site]), and 48 nubbins were generated from the 12 colonies, all of which were of the same genotype [13].Half of the nubbins from the six colonies of each site were randomly assigned to the stable temperature aquaria maintained at 26 • C, whereas the other half were placed into those aquaria of the variable temperature treatment, which fluctuated from 23 to 29 • C over a 6-h period.The 12 tanks (n = 3 for each interaction group) each contained four nubbins, two of which were sampled at time 0 (while all tanks were still at the acclimation temperature of 26 • C) and two of which were sampled after seven days of treatment exposure; only the later 24 samples are discussed herein.In the case of the physiological response variables (discussed below), a tank average was calculated across the two pseudo-replicates sampled at the same time, resulting in a total sample size of 12 only.The molecular-scale data (n = 17 parameters) were left un-pooled for MANOVA, but not for PCA.MDS was used with the 12 pooled samples after having incorporated all 23 response variables.These 23 parameters were previously assessed individually with two-way ANOVAs to determine the effects of site of origin, temperature regime, and their interaction [13][14][15].

Response Variables
The same 17 molecular-scale response variables were assessed in the samples of each experiment and included three biological composition parameters: (1) the Symbiodinium genome copy proportion (GCP; a proxy for cell density [34]); (2) the RNA/DNA ratio (a proxy for total transcription); and (3) the protein/DNA ratio (a proxy for total translation).Expression of 6 Symbiodinium and 8 host mRNAs was also quantified in each sample.The Symbiodinium genes spanned three cellular processes: photosynthesis, metabolism, and the stress response.The photosynthesis genes included ribulose-1,5-bisphosphate carboxylase/oxygenase large subunit (rbcL), photosystem I (subunit III; psI), and phosphoglycolate phosphatase (pgpase).The lone metabolism gene was nitrate transporter-2 (nrt2), and the two stress genes were ascorbate peroxidase (apx1) and hsp70.The host genes spanned the cellular processes of (1) cytoskeleton, (2) stress response, and (3) transport processes involved in osmoregulation.The former four were actb, tuba, tropomyosin (trp1), and ezrin.The three osmoregulation genes were transient receptor cation channel (trcc), organic anion transporter (oatp), and phospholipase α2 (cplap2).The lone stress gene was hsp70.The SHVTS featured six additional response variables for a total of 23 parameters assessed [13].These included the Symbiodinium RBCL protein and five physiological response variables: growth, chlorophyll a concentration (chl a; normalized two different ways (areal and per cell)), Symbiodinium density, and the maximum quantum yield of photosystem II (Fv/Fm).

MSA
Data normality and homogeneity of variance were first tested with JMP ® (version 12, Cary, NC, USA) using Shapiro-Wilk W and Levene's tests, respectively, and transformations (typically log or square root) were performed when these assumptions were not met (p < 0.05 for either test).Except for MDS, JMP was used for all statistical analyses.Prior to analysis, all data were converted to Z-scores to account for the various parameters having different units of measure (online supplemental data file 1).First, heat maps were created to portray the relative levels of the 17 molecular response variables only (Figure 1).PCA using a correlation matrix was then used to depict variation in two dimensions for each of the two experiments, and it was hypothesized that meaningful groupings of samples might be unveiled with this approach alone.PCA finds combinations of response variables that account for the greatest proportion of the variation in a dataset via their condensation into a single dimension.The second principal component (PC) is constrained by being orthogonal to the first, and only the first two PCs were considered herein (online supplemental data file 2).PCA was conducted with a variety of different combinations of samples and response variables (Table 2) to find the minimum number of parameters that could visibly partition samples by temperature or time in the SHSTTE and by temperature regime or site of origin in the SHVTS.
MANOVA was next used to determine whether the multivariate means of the experimental treatment groups could be quantitatively separated within a 2-dimensional dataspace.Briefly, while univariate ANOVA tests each dependent variable in isolation from all others, MANOVA tests multiple dependent variables simultaneously, generating a multivariate vector of univariate means (known as a "centroid") that depicts variation in multiple dimensions along canonical axes (CA).Not only can MANOVA determine differences between experimental treatments for multiple dependent variables, but it can also identify response variables that contribute most significantly to such differences.When sufficient replicates existed for the comparison of interest, Wilks' lambda was calculated, and p values < 0.05 were considered to represent significance.
For the SHSTTE, the effects of treatment temperature and time were tested by MANOVA (n = 17 parameters).The interaction effect was analyzed separately as a one-way ANOVA across the eight temperature × time interaction groups; briefly, a fully factorial, single-model MANOVA of temperature, time, and temperature × time could not be performed due to missing data in the control-24-h treatment (one replicate only).JMP instead used a pooled estimate of variance for the control-24-h sample to compute the test statistics in this one-way MANOVA, and Roy's max root was calculated in place of Wilk's lambda for this comparison.For the SHVTS, MANOVA tested the effects of temperature, site of origin, and their interaction (n = 23 parameters).MANOVA was also performed with subsets of response variables in each experiment: molecular parameters only (n = 17), the Symbiodinium molecular response only (n = 6-7 parameters), the host mRNAs only (n = 8 mRNAs), the physiological variables only (SHVTS only; n = 4 parameters), and photosynthesis parameters only (n = 3-6 parameters; Table 2).
For both experimental datasets, PRIMER (ver.5) was used to perform a MDS analysis using the Bray-Curtis similarity matrix on Z-score-transformed data.In an ordination-based MDS plot, the distance between two points is inversely proportional to their similarity; adjacent points are nearly identical, while those at opposite ends of the plot differ to the greatest extent.Analysis of similarity (ANOSIM) was then performed to determine the effects of time, temperature, and their interaction in the SHSTTE and the effects of site of origin, temperature treatment, and their interaction in the SHVTS.Global R distribution p-values were considered significant at an α of 0.05 based on 999 permutations.Finally, JMP's predictor screening function was used to rank the response variables in terms of their proportional contribution to the cumulative difference between sampling times and temperature in the SHSTTE and between site of origin and temperature regime in the SHVTS; only those parameters contributing >10% of the cumulative difference are discussed herein.
Table 2. Summary of comparisons and major findings.Dominant loading factors (principal components analysis (PCA)) and canonical correlations (multivariate analysis of variance (MANOVA)) were only included in the right-most column when the respective technique resulted in good partitioning of the data."Symbiodinium molecular response" included the genome copy proportion (GCP) + 6 mRNAs for all comparisons except for PCA of the S. hystrix variable temperature study (SHVTS), in which case the ribulose-1,5-bisphosphate carboxylase/oxygenase large subunit (RBCL) protein was also included.In certain cases, Wilks' lambda or comparable MANOVA-based statistics could not be calculated due to having a large number of response variables relative to observations.Only when the Symbiodinium data were included could temporal partitioning (i.e., "discrimination") be achieved in the S. hystrix short-term temperature experiment (SHSTTE); in contrast, data from the eight host coral genes were required to separate samples by temperature treatment (TT) and site of origin (SO) × TT in the SHVTS.PCA eigenvectors and eigenvalues can be found in online supplemental data file 2. * statistically significant observation.When p-values approached significance (α = 0.05), the exact values have been included.HBH = Houbihu.HWN = Houwan.stab = stable.var = variable.4).

Overview of the Dataset
Of the 24 samples in each of the SHSTTE and SHVTS, 4 and 5, respectively, were excluded from most MSA due to missing data.Some such samples are evident from heat maps of the SHSTTE (Figure 1a) and SHVTS (Figure 1b), and the associated data were generally lacking due to failed nucleic acid extractions.All data (as Z-scores) can be found in online supplemental data file 1.For details of the PCA (discussed below), please see online supplemental data file 2.

Overview of the Dataset
Of the 24 samples in each of the SHSTTE and SHVTS, 4 and 5, respectively, were excluded from most MSA due to missing data.Some such samples are evident from heat maps of the SHSTTE (Figure 1a) and SHVTS (Figure 1b), and the associated data were generally lacking due to failed nucleic acid extractions.All data (as Z-scores) can be found in online supplemental data file 1.For details of the PCA (discussed below), please see online supplemental data file 2. hystrix variable temperature study (SHVTS) datasets.The relative scale extends from dark blue (Zscore = −3) to dark red (Z-score = 5) in the SHSTTE (a) and SHVTS (b) heat maps."X"s denote missing data.Samples marked by asterisks (*) were excluded from most multivariate statistical approaches (MSA).Likewise, one sample from each of the Houbihu (HBH)-variable (var) and Houwan (HWN)var groups was excluded from the heat maps themselves due to the respective RNA extractions having failed.Please see the main text for full names of the target genes.C = control temperature treatment.H = high temperature treatment.stab = stable.

PCA
PCA was first performed on all 17 molecular response variables across 20 of the 24 samples (Figure 2a), and the first two eigenvectors captured only ~40% of the variation.However, corals sacrificed at the 6-h sampling time appeared to be partitioned away from the others.It was hypothesized that a reduction in the number of parameters could lead to eigenvectors collectively encompassing a greater percentage of the variation in the dataset.When looking only at the seven Symbiodinium parameters (GCP + 6 mRNAs), the first and second PC encompassed ~67% of the variation (Figure 2b), and the response variable contributing the loading score with the highest positive value in PC1 was the Symbiodinium GCP.The second PC was dominated by the

PCA
PCA was first performed on all 17 molecular response variables across 20 of the 24 samples (Figure 2a), and the first two eigenvectors captured only ~40% of the variation.However, corals sacrificed at the 6-h sampling time appeared to be partitioned away from the others.It was hypothesized that a reduction in the number of parameters could lead to eigenvectors collectively encompassing a greater percentage of the variation in the dataset.When looking only at the seven Symbiodinium parameters (GCP + 6 mRNAs), the first and second PC encompassed ~67% of the variation (Figure 2b), and the response variable contributing the loading score with the highest positive value in PC1 was the Symbiodinium GCP.The second PC was dominated by the photosynthesis genes (excluding rbcL), and one such gene, pgpase, was negatively correlated with the GCP (linear regression t-test p < 0.001; R 2 = 0.33).PCA of the Symbiodinium parameters only did not partition the data by treatment temperature or time, and the data of both temperatures and all four sampling times were inter-mixed.
When looking at the host coral mRNAs only (Figure 2c), the first two PCs explained a similar percentage of the variation (~65%) as did the first two Symbiodinium PCs (Figure 2b); furthermore, as when looking only at the Symbiodinium response variables, samples were not well separated by temperature or time.However, the 6-h data appear to be somewhat distinct from the others, with PC1 accounting for this separation.The dominant loading factors for PC1 were two cytoskeleton genes (ezrin and trp1) and two osmoregulation genes (trcc and cplap2).There was a significant degree of correlation between the expression levels of several cytoskeleton genes, as evidenced by the similar trajectory of their biplot axes, and three of the four cytoskeleton genes (excluding ezrin) contributed most significantly to PC2 in terms of eigenvector loading scores (online supplemental data file 2).
To determine whether extensive variation within treatments (i.e., a tank effect) accounted, in part, to the failure to document significant partitioning by temperature, data were pooled across triplicate tanks within each of the eight temperature × time interaction groups (Figure 2d).However, the first two PCs accounted for less than 60% of the variation, and the four data points of each temperature treatment were essentially mixed with those of the other.Despite this, the time 6-and 24-h data were somewhat distinct, as was evident when all data were considered (Figure 2a).photosynthesis genes (excluding rbcL), and one such gene, pgpase, was negatively correlated with the GCP (linear regression t-test p < 0.001; R 2 = 0.33).PCA of the Symbiodinium parameters only did not partition the data by treatment temperature or time, and the data of both temperatures and all four sampling times were inter-mixed.
When looking at the host coral mRNAs only (Figure 2c), the first two PCs explained a similar percentage of the variation (~65%) as did the first two Symbiodinium PCs (Figure 2b); furthermore, as when looking only at the Symbiodinium response variables, samples were not well separated by temperature or time.However, the 6-h data appear to be somewhat distinct from the others, with PC1 accounting for this separation.The dominant loading factors for PC1 were two cytoskeleton genes (ezrin and trp1) and two osmoregulation genes (trcc and cplap2).There was a significant degree of correlation between the expression levels of several cytoskeleton genes, as evidenced by the similar trajectory of their biplot axes, and three of the four cytoskeleton genes (excluding ezrin) contributed most significantly to PC2 in terms of eigenvector loading scores (online supplemental data file 2).
To determine whether extensive variation within treatments (i.e., a tank effect) accounted, in part, to the failure to document significant partitioning by temperature, data were pooled across triplicate tanks within each of the eight temperature × time interaction groups (Figure 2d).However, the first two PCs accounted for less than 60% of the variation, and the four data points of each temperature treatment were essentially mixed with those of the other.Despite this, the time 6-and 24-h data were somewhat distinct, as was evident when all data were considered (Figure 2a).)).All 24 samples were included in the PCA of the Symbiodinium response variables only (GCP + 6 genes; (b)).For the eight host genes (c), the same four samples as in (a) were omitted, and the ends of the vectors representing three cytoskeleton genes have been encircled to emphasize their correlation.Data were also analyzed as pooled across aquaria for all 17 parameters for each of the eight temperature × time interaction groups (d).The legend in (a) applies to all panels.

MANOVA
When using MANOVA to test the interaction of temperature and time (Figure 3a) with data from all 17 response variables, Roy's max root was statistically significant, and this is likely due to the wide separation of samples of times 6 and 24 h along CA1.This partitioning was driven by a negative relationship between Symbiodinium pgpase and apx1 mRNA expression (Table 2).Within the 6-and 24-h centroids, the high temperature samples were reasonably well separated from the controls; it should be noted, though, that only one sample comprised the control-24-h group due to missing data.www.mdpi.com/journal/jmse3.2.2.MANOVA When using MANOVA to test the interaction of temperature and time (Figure 3a) with data from all 17 response variables, Roy's max root was statistically significant, and this is likely due to the wide separation of samples of times 6 and 24 h along CA1.This partitioning was driven by a negative relationship between Symbiodinium pgpase and apx1 mRNA expression (Table 2).Within the 6-and 24-h centroids, the high temperature samples were reasonably well separated from the controls; it should be noted, though, that only one sample comprised the control-24-h group due to missing data.

Figure 3. Multivariate analysis of variance (MANOVA) and multidimensional scaling (MDS) analysis of the
Seriatopora hystrix short-term temperature experiment (SHSTTE).MANOVA was used to test for differences between the eight temperature × time groups (a), and 4 of the 24 total samples were omitted due to missing data points (see Figure 1.); therefore, Roy's max root, rather than Wilks' lambda, was calculated in (a).In (a), the control-12-h (C12) and high temperature-12-h (H12) centroids (black circles; 95% confidence intervals for these and all centroids presented in this manuscript) lie within the high temperature-48-h (H48; blue circle) centroid and have been labeled with red lines.Not all axes have been presented or labeled due to spatial constraints.MDS analysis of the same 20 samples with PRIMER (b).In this panel only, circles were drawn by hand and do not represent 95% confidence centroids.The legend for both panels lies in the lower-right corner of (b).* p < 0.05.TT = temperature treatment.

MDS and Predictor Screening
Upon analyzing the SHSTTE data with PRIMER, samples were significantly affected/sorted by time (ANOSIM Global R p = 0.002), but not by temperature (p > 0.05).It should be noted, however, that the stress was quite high (0.23) in the MDS plot (Figure 3b) and so this temporal variation should be interpreted cautiously.Specifically, although the 12-and 48-h samples were intermixed, the 6-and 24-h times were well separated, as was also evidenced by MANOVA (Figure 3a), and, to some extent, PCA (Figure 2a,c).PCA, MANOVA, and MDS all appear to suggest, then, that time, and not temperature, was more important in accounting for variation in the SHSTTE dataset.Therefore, the predictor screening function of JMP was used to identify the response variables that explained the greatest proportion of the differences between sampling times (Figure 4a), and these were found to be the protein/DNA ratio (21% of the cumulative difference), Symbiodinium pgpase mRNA expression (14%), and the RNA/DNA ratio (13%).MANOVA (Figure 3a) also found Symbiodinium pgpase to contribute to the separation of samples by time, specifically by distinguishing those samples of the 24-h sampling time (Table 2).
Regarding temperature (Figure 4a), only three response variables contributed to 10% or more of the documented cumulative difference in the molecular physiology of the control and high temperature samples: host trp1 (17%), Symbiodinium apx1 (14%), and the RNA/DNA ratio (13%).However, the expression of neither gene, nor the RNA/DNA ratio, differed significantly between MANOVA was used to test for differences between the eight temperature × time groups (a), and 4 of the 24 total samples were omitted due to missing data points (see Figure 1.); therefore, Roy's max root, rather than Wilks' lambda, was calculated in (a).In (a), the control-12-h (C12) and high temperature-12-h (H12) centroids (black circles; 95% confidence intervals for these and all centroids presented in this manuscript) lie within the high temperature-48-h (H48; blue circle) centroid and have been labeled with red lines.Not all axes have been presented or labeled due to spatial constraints.MDS analysis of the same 20 samples with PRIMER (b).In this panel only, circles were drawn by hand and do not represent 95% confidence centroids.The legend for both panels lies in the lower-right corner of (b).* p < 0.05.TT = temperature treatment.

MDS and Predictor Screening
Upon analyzing the SHSTTE data with PRIMER, samples were significantly affected/sorted by time (ANOSIM Global R p = 0.002), but not by temperature (p > 0.05).It should be noted, however, that the stress was quite high (0.23) in the MDS plot (Figure 3b) and so this temporal variation should be interpreted cautiously.Specifically, although the 12-and 48-h samples were intermixed, the 6-and 24-h times were well separated, as was also evidenced by MANOVA (Figure 3a), and, to some extent, PCA (Figure 2a,c).PCA, MANOVA, and MDS all appear to suggest, then, that time, and not temperature, was more important in accounting for variation in the SHSTTE dataset.Therefore, the predictor screening function of JMP was used to identify the response variables that explained the greatest proportion of the differences between sampling times (Figure 4a), and these were found to be the protein/DNA ratio (21% of the cumulative difference), Symbiodinium pgpase mRNA expression (14%), and the RNA/DNA ratio (13%).MANOVA (Figure 3a) also found Symbiodinium pgpase to contribute to the separation of samples by time, specifically by distinguishing those samples of the 24-h sampling time (Table 2).
Regarding temperature (Figure 4a), only three response variables contributed to 10% or more of the documented cumulative difference in the molecular physiology of the control and high temperature samples: host trp1 (17%), Symbiodinium apx1 (14%), and the RNA/DNA ratio (13%).However, the expression of neither gene, nor the RNA/DNA ratio, differed significantly between temperatures despite the fact that the expression of trp1 was 4-fold less in high temperature samples.When looking at the interaction of time and temperature (Figure 4a), the RNA/DNA ratio and Symbiodinium apx1 were the only factors that contributed to greater than 10% of the total difference between the four interactions groups (22% and 12%, respectively); neither showed an interaction effect when analyzed by repeated measures ANOVA (p > 0.05; [13]).
temperatures despite the fact that the expression of trp1 was 4-fold less in high temperature samples.When looking at the interaction of time and temperature (Figure 4a), the RNA/DNA ratio and Symbiodinium apx1 were the only factors that contributed to greater than 10% of the total difference between the four interactions groups (22 and 12%, respectively); neither showed an interaction effect when analyzed by repeated measures ANOVA (p > 0.05; [13]).Predictor screening analysis of the Seriatopora hystrix short-term temperature experiment (SHSTTE) and the S. hystrix variable temperature study (SHVTS).The predictor screening function of JMP was used to determine the response variables that accounted for the greatest proportions of the cumulative differences between temperature, time, and their interaction in the SHSTTE (17 response variables; (a)) and between temperature regime, site of origin, and their interaction in the SHVTS (23 response variables; (b)).It should be noted that the relative scales differ between the proportion plots.

PCA
PCA, MANOVA, and MDS were all able to separate the samples by temperature regime in the SHVTS, and some approaches were able to separate the four site of origin × temperature treatment interaction groups (Table 2).First, PCA distinguished a variety of informative groupings (Figure 5a).When looking at all 23 response variables for data pooled across pseudo-replicates for each of the 12   Predictor screening analysis of the Seriatopora hystrix short-term temperature experiment (SHSTTE) and the S. hystrix variable temperature study (SHVTS).The predictor screening function of JMP was used to determine the response variables that accounted for the greatest proportions of the cumulative differences between temperature, time, and their interaction in the SHSTTE (17 response variables; (a)) and between temperature regime, site of origin, and their interaction in the SHVTS (23 response variables; (b)).It should be noted that the relative scales differ between the proportion plots.

PCA
PCA, MANOVA, and MDS were all able to separate the samples by temperature regime in the SHVTS, and some approaches were able to separate the four site of origin × temperature treatment interaction groups (Table 2).First, PCA distinguished a variety of informative groupings (Figure 5a).When looking at all 23 response variables for data pooled across pseudo-replicates for each of the 12 aquaria, there was a clear separation of the stable and variable temperature samples along both PCs.However, the total variation encompassed by these two PCs was less than 55%.PC1 was comprised of a mix of host and Symbiodinium genes (Table 2), while the second PC consisted mainly of Symbiodinium photosynthesis genes (online supplemental data file 2).There was some partitioning by site of colony origin within each temperature treatment, though still some overlap.Clearly, the effect of temperature on the molecular physiology of S. hystrix was greater than that due to site of origin upon assessment of PCA alone.
aquaria, there was a clear separation of the stable and variable temperature samples along both PCs.However, the total variation encompassed by these two PCs was less than 55%.PC1 was comprised of a mix of host and Symbiodinium genes (Table 2), while the second PC consisted mainly of Symbiodinium photosynthesis genes (online supplemental data file 2).There was some partitioning by site of colony origin within each temperature treatment, though still some overlap.Clearly, the effect of temperature on the molecular physiology of S. hystrix was greater than that due to site of origin upon assessment of PCA alone.
Figure 5. Principal components analysis (PCA) of the Seriatopora hystrix variable temperature study (SHVTS).Data from pseudo-replicate samples of the same tank were pooled, resulting in 12 data points in each plot.All 23 parameters (a).The Symbiodinium molecular response only (genome copy proportion (GCP) + 6 mRNAs + the ribulose-1,5-bisphosphate carboxylase/oxygenase (large subunit) protein (RBCL); (b)).One Houwan-stable icon is nearly masked by a Houbihu-stable one towards the left side of the plot.The host molecular response only (8 mRNAs; (c)).In (c), the unlabeled phospholipase α2 (cplap2) axis is below the ezrin one.Likewise, the tropomyosin (trp1) and β-actin (actb) axes are overlapping.The six photosynthesis parameters only (d).The legend for all panels is located in (c).
In order to reduce the complexity of the dataset, PCA was also conducted only with the eight Symbiodinium molecular response variables: the GCP, the six mRNAs, and the RBCL protein (Figure 5b).Samples of the two temperature regimes, stable and variable, could be distinguished by PC1 (51.9%), which was dominated by the photosynthesis genes in terms of highest positive eigenvector loading factor scores (Table 2 and online supplemental data file 2).The stable temperature samples tended to cluster together, with the six variable ones showing greater variability and spread throughout the dataspace.Furthermore, the Houwan variable temperature samples showed greater dispersal than did the Houbihu ones (Figure 5b).The Symbiodinium GCP was the most significant contributor to the second PC (18.7%).When looking at individual correlations of photosynthesis gene expression vs. Symbiodinium GCP (data not shown), all slopes were negative; however, unlike for the Figure 5. Principal components analysis (PCA) of the Seriatopora hystrix variable temperature study (SHVTS).Data from pseudo-replicate samples of the same tank were pooled, resulting in 12 data points in each plot.All 23 parameters (a).The Symbiodinium molecular response only (genome copy proportion (GCP) + 6 mRNAs + the ribulose-1,5-bisphosphate carboxylase/oxygenase (large subunit) protein (RBCL); (b)).One Houwan-stable icon is nearly masked by a Houbihu-stable one towards the left side of the plot.The host molecular response only (8 mRNAs; (c)).In (c), the unlabeled phospholipase α2 (cplap2) axis is below the ezrin one.Likewise, the tropomyosin (trp1) and β-actin (actb) axes are overlapping.The six photosynthesis parameters only (d).The legend for all panels is located in (c).
In order to reduce the complexity of the dataset, PCA was also conducted only with the eight Symbiodinium molecular response variables: the GCP, the six mRNAs, and the RBCL protein (Figure 5b).Samples of the two temperature regimes, stable and variable, could be distinguished by PC1 (51.9%), which was dominated by the photosynthesis genes in terms of highest positive eigenvector loading factor scores (Table 2 and online supplemental data file 2).The stable temperature samples tended to cluster together, with the six variable ones showing greater variability and spread throughout the dataspace.Furthermore, the Houwan variable temperature samples showed greater dispersal than did the Houbihu ones (Figure 5b).The Symbiodinium GCP was the most significant contributor to the second PC (18.7%).When looking at individual correlations of photosynthesis gene expression vs. Symbiodinium GCP (data not shown), all slopes were negative; however, unlike for the pgpase gene in the SHSTTE, these slopes were not significantly different from 0 (linear regression t-tests, p > 0.05).
When looking at the host coral molecular response only (Figure 5c), PCA also partitioned samples of the two temperature regimes.The first PC encompassed all three cytoskeleton genes and four osmoregulation genes and explained 62.5% of the variation (Table 2).The second PC consisted of hsp70 as the only positive loading score (online supplemental data file 2), and this PC encompassed 16.3% of the variation.The expression levels of the osmoregulation and cytoskeleton genes tended to correlate.As with the PCA of the Symbiodinium molecular response variables only (Figure 5b), the spread of the variable temperature regime samples was greater than that of the stable temperature ones.In contrast to the results of the PCA of the Symbiodinium molecular response only, though, the Houbihu variable temperature samples showed more spread than those of Houwan exposed to this same temperature profile.Finally, when looking only at the six photosynthesis parameters (Figure 5d), samples of the two temperature regimes were well separated along PC1 (55.2%), in which the psI gene and the RBCL protein contributed the highest positive loading scores (Table 2).Samples of the two sites of origin for the stable temperature treatment were separated along PC2, in which case chl a content (pg/cell) contributed most substantially to such partitioning (online supplemental data file 2).

MANOVA
MANOVA was performed with a variety of combinations of response variables to determine which best modeled differences between the four site of origin × temperature treatments groups (Table 2).First, MANOVA was performed with four of the five physiological response variables alone (areal chl a was excluded in place of chl a/cell; Figure 6a).Although Wilks' lambda for the interaction effect was significant, only two of the four groups appear well-separated: Houbihu-variable and Houwan-stable.Houbihu-stable and Houwan-variable samples were inter-mixed.Fv/Fm was the most significant factor contributing to the separation of the former two groups (Table 2 and online supplemental data file 2).When looking at all 17 molecular response variables with 20 of the 24 samples, all four interaction groups were well separated (Figure 6b), though Wilks' lambda could not be computed due to the large number of response variables relative to the number of biological replicates (~1:1).A negative relationship between host actb and Symbiodinium psI + protein/DNA drove the partitioning of the four interactions groups (Table 2), though such separation was not considered statistically significant by Roy's max root (p > 0.05).Figure 6.Multivariate analysis of variance (MANOVA) and multidimensional scaling (MDS) analysis of the Seriatopora hystrix variable temperature study (SHVTS).MANOVA was used to model differences between the four site of origin (SO) × temperature treatment (TT) interaction groups with four of the five physiological parameters (excluding areal chlorophyll a (chla a); (a)), all 17 molecular response parameters (b), seven of the eight Symbiodinium (Sym) molecular response variables (excluding the ribulose-1,5-bisphosphate carboxylase/oxygenase (large subunit) protein (RBCL); (c)), all eight host coral mRNAs (d), and all 23 response variables (e).The sample sizes displayed in the individual panels reflect the number of samples and not the number of response variables.Wilks' lambda could not be computed in (b) and (e) due to the large number of response variables relative to the number of samples.In (b), not all axes have been labeled.In (c), the Sym genome copy proportion (GCP) axis falls between those of the mRNAs nitrate transporter-2 (nrt2) and rbcL.In (e), only the dominant axes/response variables have been shown due to spatial constraints in the panel itself.In (a-d) Wilks' lambda (a, c, d) and Roy's max root (b) values test the interaction of SO and TT.In (f), PRIMER's MDS function was used to portray the SHVTS dataspace, and the circles were handdrawn to encompass the four SO × TT groups.All other circles represent 95% confidence centroids.The legend in (e) corresponds to all panels.* Global R p < 0.01.
When analyzing seven of the eight Symbiodinium molecular response variables only (RBCL protein expression was excluded since it was only assessed in 12 of the 24 samples; Figure 6c), the four interaction groups were not found to differ significantly by MANOVA.In contrast, when looking at the eight host genes alone (Figure 6d), significant discrimination was achieved.However, Figure 6.Multivariate analysis of variance (MANOVA) and multidimensional scaling (MDS) analysis of the Seriatopora hystrix variable temperature study (SHVTS).MANOVA was used to model differences between the four site of origin (SO) × temperature treatment (TT) interaction groups with four of the five physiological parameters (excluding areal chlorophyll a (chla a); (a)), all 17 molecular response parameters (b), seven of the eight Symbiodinium (Sym) molecular response variables (excluding the ribulose-1,5-bisphosphate carboxylase/oxygenase (large subunit) protein (RBCL); (c)), all eight host coral mRNAs (d), and all 23 response variables (e).The sample sizes displayed in the individual panels reflect the number of samples and not the number of response variables.Wilks' lambda could not be computed in (b,e) due to the large number of response variables relative to the number of samples.In (b), not all axes have been labeled.In (c), the Sym genome copy proportion (GCP) axis falls between those of the mRNAs nitrate transporter-2 (nrt2) and rbcL.In (e), only the dominant axes/response variables have been shown due to spatial constraints in the panel itself.In (a-d) Wilks' lambda (a, c, d) and Roy's max root (b) values test the interaction of SO and TT.In (f), PRIMER's MDS function was used to portray the SHVTS dataspace, and the circles were hand-drawn to encompass the four SO × TT groups.All other circles represent 95% confidence centroids.The legend in (e) corresponds to all panels.* Global R p < 0.01.
When analyzing seven of the eight Symbiodinium molecular response variables only (RBCL protein expression was excluded since it was only assessed in 12 of the 24 samples; Figure 6c), the four interaction groups were not found to differ significantly by MANOVA.In contrast, when looking at the eight host genes alone (Figure 6d), significant discrimination was achieved.However, the two sites of origin were only well separated in the variable temperature dataset and were inter-mixed for the stable temperature regime samples.When looking at all 23 response variables, the four site of origin × temperature treatment interaction groups were well separated (Figure 6e), though neither Wilks' lambda nor Roy's max root could be calculated due to the large number of response variables relative to observations.

MDS and Predictor Screening
MDS can readily depict relationships between samples even when the number of response variables is large relative to the number of observations, and MDS was able to distinguish the four interaction groups of the SHVTS (Figure 6f).ANOSIM found site of origin, temperature treatment, and interaction effects to be significant, and, within the stable temperature regime samples, data points of the two sites of origin were well separated.However, there was some overlap between the Houbihu-variable and Houwan-variable samples due to one Houbihu-variable data point falling closer to those points of the latter group.JMP's predictor screening function (Figure 4b) found the five physiological parameters and Symbiodinium psI to contribute most significantly to the cumulative difference between the four interaction groups (Table 2).psI also contributed significantly to the cumulative difference between temperature regimes (Figure 4b), in conjunction with Fv/Fm and host hsp70.

Discussion
Corals of Southern Taiwan have proven to be resilient to a number of laboratory-induced environmental challenges (Table 1), and MSA confirmed this univariate ANOVA-based observation for samples of the SHSTTE.This leaves at least two hypotheses remaining: (1) the corals were truly unresponsive to a short-term exposure to 30 • C or (2) non-temperature-sensitive parameters were chosen.Given the well-documented photoinhibition that occurs when Symbiodinium are exposed to elevated temperatures (e.g., [35]), it seems likely that at least several photosynthesis genes should have undergone changes in mRNA expression.However, dinoflagellates exhibit a significant degree of post-transcriptional activity [36], suggesting that a proteomic analysis of thermally-challenged coral-Symbiodinium holobionts may ultimately be required to conclusively understand the sub-cellular basis of Symbiodinium acclimation to elevated or variable temperatures.
Corals of the four sampling times possessed unique gene expression + biological composition signatures.Temporal changes in the molecular phenotype of a coral have been documented previously [37] and may be related to the complex, dual-compartmental metabolism displayed by organisms, such as reef-building corals, that have acquired the capacity for photosynthesis via symbiosis [38,39].Specifically, coral metabolism is temporally variable due to the periodic nature of light-driven photosynthesis [38][39][40][41], and metabolic hysteresis driven by dinoflagellate photosynthesis as a function of the light cycle likely contributed to the temporal variation observed in the SHSTTE.Circadian rhythm may also have accounted, in part, for the separation of the 6-and 24-h samples in the SHSTTE [42].The former were collected at 13:45, and, while stable, artificial light was used in this experiment, it is possible that the temporal gene expression signatures were driven by an entrained response to high light levels that would normally be experienced at such times.The 12, 24, and 48-h sampling times corresponded to 19:30, 7:30, and 7:30, respectively, times at which light levels would be relatively low in situ.However, the experimental corals had been reared under non-fluctuating, artificial light (12:12-h light-dark cycle) for nearly one month at the time of sampling (including the pre-experimental acclimation period), and circadian rhythm can be abolished within two days of changing the light regime in endosymbiotic anthozoans [39].Therefore, other factors besides metabolic hysteresis due to photosynthesis and circadian rhythm may have accounted for the unique molecular phenotype of corals sampled after 6 h.
All MSA documented temperature regime, and oftentimes site of origin, differences in the SHVTS, and the host coral parameters were more likely to partition samples by temperature in the SHVTS; this contrasts with the SHSTTE, in which the Symbiodinium response variables were relatively more important in the separation of samples across the dataspace.However, the predominant experimental factor leading to sample partitioning in the SHSTTE was time, rather than temperature.The fact that the Symbiodinium response variables more significantly contributed to temporal variation in the SHSTTE, while host coral parameters led to a relatively greater partitioning of samples by temperature in the SHVTS, emphasizes that notion put forth by Mayfield and Gates [18] and others that it is important to consider both compartments of the coral-Symbiodinium endosymbiosis when attempting to gauge the molecular physiological response of the composite holobiont to environmental change.
When performing PCA on the Symbiodinium molecular response only, the Houwan variable temperature samples showed greater dispersal than did the Houbihu ones.This could be because the former corals had never before experienced such variable temperature profiles in situ; as such, the variability in their response to fluctuating temperatures could be hypothesized to be greater than that of corals of Houbihu, which do routinely experience upwelling.Likewise, when looking at the MDS plot of the SHVTS dataset, the spread of the variable temperature samples was greater than that of the stable ones, and a similar explanation could account for this observation.Although variability in the physiological response to an environmental change has indeed been predicted to be important in gauging marine animal health [43], it is unlikely that such enhanced phenotypic variability in response to variable temperature exposure documented herein is indicative of stress, as corals of both sites of origin ultimately grew at normal rates under variable temperature conditions [13].However, the unique molecular physiological signatures of the four interaction groups suggest that the sub-cellular and physiological acclimation responses to fluctuating temperatures differed between corals of the two sites or origin.Future work will attempt to more thoroughly characterize the molecular constituents responsible for such variable temperature acclimation for corals of Houwan and Houbihu, as doing so could help elucidate the role of environmental history in coral acclimation to future environmental shifts.

Conclusions
MSA were successfully used to define time-specific phenotypes in the SHSTTE and molecular physiologies with fidelity to both site of origin and temperature regime in the SHVTS.MDS proved to be an especially well-suited means of displaying a molecular phenotype that integrates a number of different response variables and macromolecules as, unlike MANOVA, the MDS-associated ANOSIM can still be performed when the number of response variables is large relative to the number of samples.Rather than the absolute expression level of a gene mRNA characterizing a sample group, relationships between multiple response variables and genes instead better distinguished corals sampled at different times in the SHSTTE.In the SHVTS, multiple groupings of response variables (e.g., gene expression, physiological, and biological composition parameters) partitioned samples by temperature regime, and the molecular physiological phenotypes differed significantly between corals of these two temperature profiles.Unlike the SHSTTE, in which the Symbiodinium response was a greater contributor to the overall variation of the dataset, host coral response variables better partitioned data points of the two temperature regimes in the SHVTS.Finally, corals exposed to variable temperatures showed a greater range in their molecular physiological response relative to those exposed to stable temperature, as had been hypothesized.Future works with these same samples will exploit next generation transcriptome sequencing and mass spectrometry-based protein sequencing to attempt to gain further insight as to how these corals acclimate to highly variable temperature regimes in situ, as the underlying molecular pathways are likely to be those that these corals will use to combat the rising temperatures that will characterize their environments in the coming years.

Figure 1 .
Figure 1.Heat maps of the Seriatopora hystrix short-term temperature experiment (SHSTTE) and S. hystrix variable temperature study (SHVTS) datasets.The relative scale extends from dark blue (Zscore = −3) to dark red (Z-score = 5) in the SHSTTE (a) and SHVTS (b) heat maps."X"s denote missing data.Samples marked by asterisks (*) were excluded from most multivariate statistical approaches (MSA).Likewise, one sample from each of the Houbihu (HBH)-variable (var) and Houwan (HWN)var groups was excluded from the heat maps themselves due to the respective RNA extractions having failed.Please see the main text for full names of the target genes.C = control temperature treatment.H = high temperature treatment.stab = stable.

Figure 2 .
Figure2.Principal components analysis (PCA) of the Seriatopora hystrix short-term temperature experiment (SHSTTE) dataset.All 17 response variables (a), including 3 biological composition parameters (RNA/DNA ratio, protein/DNA ratio, and the Symbiodinium genome copy proportion (GCP)), expression of 6 Symbiodinium mRNAs, and expression of 8 coral mRNAs, were assessed across 20 of the 24 samples (4 samples were omitted due to missing data (see Figure1.)).All 24 samples were included in the PCA of the Symbiodinium response variables only (GCP + 6 genes; (b)).For the eight host genes (c), the same four samples as in (a) were omitted, and the ends of the vectors representing three cytoskeleton genes have been encircled to emphasize their correlation.Data were also analyzed as pooled across aquaria for all 17 parameters for each of the eight temperature × time interaction groups (d).The legend in (a) applies to all panels.

Figure 2 .
Figure2.Principal components analysis (PCA) of the Seriatopora hystrix short-term temperature experiment (SHSTTE) dataset.All 17 response variables (a), including 3 biological composition parameters (RNA/DNA ratio, protein/DNA ratio, and the Symbiodinium genome copy proportion (GCP)), expression of 6 Symbiodinium mRNAs, and expression of 8 coral mRNAs, were assessed across 20 of the 24 samples (4 samples were omitted due to missing data (see Figure1.)).All 24 samples were included in the PCA of the Symbiodinium response variables only (GCP + 6 genes; (b)).For the eight host genes (c), the same four samples as in (a) were omitted, and the ends of the vectors representing three cytoskeleton genes have been encircled to emphasize their correlation.Data were also analyzed as pooled across aquaria for all 17 parameters for each of the eight temperature × time interaction groups (d).The legend in (a) applies to all panels.

Figure 3 .
Figure3.Multivariate analysis of variance (MANOVA) and multidimensional scaling (MDS) analysis of the Seriatopora hystrix short-term temperature experiment (SHSTTE).MANOVA was used to test for differences between the eight temperature × time groups (a), and 4 of the 24 total samples were omitted due to missing data points (see Figure1.); therefore, Roy's max root, rather than Wilks' lambda, was calculated in (a).In (a), the control-12-h (C12) and high temperature-12-h (H12) centroids (black circles; 95% confidence intervals for these and all centroids presented in this manuscript) lie within the high temperature-48-h (H48; blue circle) centroid and have been labeled with red lines.Not all axes have been presented or labeled due to spatial constraints.MDS analysis of the same 20 samples with PRIMER (b).In this panel only, circles were drawn by hand and do not represent 95% confidence centroids.The legend for both panels lies in the lower-right corner of (b).* p < 0.05.TT = temperature treatment.

Figure 4 .
Figure 4. Predictor screening analysis of the Seriatopora hystrix short-term temperature experiment (SHSTTE) and the S. hystrix variable temperature study (SHVTS).The predictor screening function of JMP was used to determine the response variables that accounted for the greatest proportions of the cumulative differences between temperature, time, and their interaction in the SHSTTE (17 response variables; (a)) and between temperature regime, site of origin, and their interaction in the SHVTS (23 response variables; (b)).It should be noted that the relative scales differ between the proportion plots.

Figure 4 .
Figure 4. Predictor screening analysis of the Seriatopora hystrix short-term temperature experiment (SHSTTE) and the S. hystrix variable temperature study (SHVTS).The predictor screening function of JMP was used to determine the response variables that accounted for the greatest proportions of the cumulative differences between temperature, time, and their interaction in the SHSTTE (17 response variables; (a)) and between temperature regime, site of origin, and their interaction in the SHVTS (23 response variables; (b)).It should be noted that the relative scales differ between the proportion plots.

Table 1 .
Summary of eight environmental challenge studies performed at Taiwan's National Museum of Marine Biology and Aquarium between 2009 and 2014.In general, only corals exposed to 31.5 • C for 2-4 weeks were found to bleach, die, and/or, more generally, display a significantly different phenotype.pCO 2 = carbon dioxide partial pressure.ppm = parts per million.NA = not applicable.var.= variable.OA = ocean acidification.
a Data not depicted graphically; b could not compute Wilks' lambda or Roy's max root; c determined by JMP's predictor screening function (Figure