Synergistic Effects of Thermal Stress and Estuarine Discharge on Transcriptomic Variation of Montastraea cavernosa Corals in Southeast Florida

Coral reefs at the northern extent of Florida’s coral reef tract are exposed to many localized anthropogenic influences including controlled freshwater discharges, runoff, upwelling, and seasonal environmental variability. To better understand coral responses to sublethal stressors in nearshore environments, we conducted complementary experiments to assess the impacts of estuarine runoff and temperature stress on local populations of the scleractinian coral species, Montastraea cavernosa, using Tag-Seq global gene expression profiling. In an in situ time series experiment, fate-tracked colonies were sampled during periods of relatively low and high estuarine discharge over 4 years to investigate temporal trends in transcriptional patterns and to identify if coral stress indicators were regulated through time. There was significant transcriptomic variation through time, but patterns did not appear to be attributed to distance from nearby estuarine tidal flux. In an ex situ factorial experiment, clonal replicates of coral genotypes were exposed to temperature (25°C and 30°C) and water (offshore and estuarine discharge, representing typical oceanic conditions and episodic discharge conditions, respectively) treatments to quantify the potential individual and synergistic effects of sublethal stress on coral and algal gene expression. Comparative analyses suggested that corals and their algal symbionts were more responsive to thermal stress than to estuarine discharge, although there was evidence of a synergistic relationship between the two stressors. Strong genotype effects also demonstrated that transcriptomic responses to thermal stress were largely based on coral genotype, indicating the potential for stress resilience among certain members of coral populations from southeast Florida.


INTRODUCTION
Increases in coastal development and changes to watershed use in recent decades have had profound implications for nearshore coral reefs. Coral reef habitats at the northern extent of Florida's coral reef tract are of particular concern, residing at their latitudinal environmental limits while also increasingly exposed to global and local anthropogenic stressors including land-based sources of pollution (e.g., freshwater runoff, sedimentation), annual temperature flux, and warming sea temperatures (Lirman and Fong, 2007;Gregg, 2013). Within this region, coral reefs inhabit an ecotone transition between subtropical and warm temperate climate zones (Lugo et al., 1999). Reefs in southeast Florida face high levels of both seasonal and episodic environmental variation that combine with anthropogenic impacts to impose multiple stressors on scleractinian corals. Land-based influences (e.g., controlled freshwater discharge and agricultural runoff), as well as natural environmental variability (e.g., seasonal rainfall and upwelling events), impact coastal habitats in the region (Yentsch et al., 2002;Beal et al., 2012;Gregg, 2013).
Two distinct climatological seasons dominate weather patterns along nearshore reefs in southeast Florida: (1) a wet season characterized by increased freshwater discharge and runoff, nutrient over-enrichment, and migration of dark, lowsalinity water over coastal zones, and (2) a subsequent dry season during the winter with decreased rainfall and runoff (Beal et al., 2012). This pattern is also influenced by global environmental variations such as the El Niño Southern Oscillation and regional tropical storm intensity and duration (Lapointe et al., 2012). During periods of prolonged rainfall, coastal habitats including nearshore coral reefs are frequently subjected to discharge from a series of natural tributaries and artificial canals which are meant to contain stormwater runoff (Graves et al., 2004;Lapointe et al., 2012). This system controls water levels in Lake Okeechobee and manages a portion of flow from the surrounding watershed via periodic releases through control structures located on major canals (e.g., C-23, C-24, C-44; Figure 1). As a result, the majority of runoff water is drained from urban and agricultural lands into coastal zones, rather than into the Everglades to the south (Perry, 2004). Characteristics of discharge water during high levels of flow from the watershed have been documented to include increased levels of total suspended solids, water color (i.e., colored dissolved organic matter; CDOM), and nutrients such as nitrogen and phosphorous (Collier et al., 2008;Lapointe et al., 2017).
Restructuring and management of the Okeechobee and St. Lucie watersheds has drastically changed environmental stability through time, whereas water quality can rapidly fluctuate near St. Lucie Inlet at the southern terminus of the Indian River Lagoon (Pickering and Baker, 2015). This region also corresponds with the northernmost portion of Florida's coral reef tract, St. Lucie Reef (SLR), located within the St. Lucie Inlet Preserve State Park (Figure 1 inset). This nearshore reef environment begins directly outside of the St. Lucie Inlet and continues approximately 7 km to the south (Walker and Gilliam, 2013). St. Lucie Reef is the northern known limit of at least eight scleractinian and octocoral species along the reef tract (Reed, 1982;Beal et al., 2012). Near the inlet, the reef is composed primarily of wormrock (Phragmatopoma lapidosa), with an increase in diversity and abundance of scleractinian corals toward the southern end of the park. The observed gradient of scleractinian diversity is hypothesized to be a result of increasing estuarine influence closer to the inlet (Beal et al., 2012;Shatters, 2017). Despite highly variable conditions, SLR supports a diverse community of 24 species of scleractinian corals, 100 other species of invertebrates, 23 algal species, four sea turtle species, and approximately 250 fish species (Beal et al., 2012). The two dominant scleractinian coral species found at SLR are Pseudodiploria clivosa and Montastraea cavernosa, with the former occurring throughout, and the latter absent from the northern third of the reef nearest the inlet (Walker and Gilliam, 2013). Their persistence on the reef despite highly variable water quality parameters may be in part due to their prolific ability to shed sediments through excessive mucus production and sloughing, heterotrophic feeding, and skeletal morphology characteristics (Bak and Elgershuizen, 1976;Lasker, 1980;Vargas-Ángel et al., 2006Vargas-Ángel et al., , 2007. There have been ongoing efforts to monitor reef populations and seasonal environmental variability within the state park since 2006 (Beal et al., 2012;Walker and Gilliam, 2013;Gilliam et al., 2018). From 2010 to 2017, quarterly surveys have assessed survival of fate-tracked colonies of the dominant scleractinian species, P. clivosa and M. cavernosa, with occasional tissue sampling events to quantify coral health through physiological, histological, and population genetics approaches. These two species have shown relatively stable coral-algal assemblages over time, suggesting resilient Symbiodiniaceae associations despite environmental and seasonal changes (Klepac et al., 2015). Coral transcriptomic analyses using microarrays, however, have documented an increase in expression of stress response genes during the wet season (Beal et al., 2012). Histological analyses have indicated that these two coral species at SLR are not actively reproductive (Beal et al., 2012), despite evidence of historical population connectivity to other reefs ∼20-50 km to the south (Dodge et al., 2020). Coral skeletons appeared less dense in fatetracked colonies relative to conspecifics elsewhere along Florida's coral reef tract, with additional tissue and skeletal loss due to the boring sponges Cliona spp. (Beal et al., 2012). In spite of these stressors, coral mortality was low at SLR from 2010 until Hurricane Irma in 2017 (Walker, 2018) and the recent outbreak of stony coral tissue loss disease (SCTLD) along Florida's coral reef tract starting in 2014 (reaching the northern regions in 2017) (Walton et al., 2018). Since 2017, monitoring efforts have increased both temporally and by individual fate-tracking to include all remaining scleractinian species on the reef (Voss and Combs, 2018;Walker, 2018). Prior to widespread coral mortality of the two dominant species at SLR due to SCTLD, colonies were remarkably persistent in the face of hypervariable conditions and stress.
This study builds on a decade of monitoring efforts to describe, at the molecular level, how corals at SLR respond to environmental variability. Using complementary temporal in situ and manipulative ex situ designs, we sought to quantify transcriptomic variation among individuals exposed to different levels of environmental stress. Recent advancements in RNA sequencing (RNA-seq) and analysis pipelines have enhanced assessments of coral health through whole-transcriptome gene expression profiling, especially in response to stressors (Kenkel et al., 2011(Kenkel et al., , 2013b(Kenkel et al., , 2014Meyer et al., 2011;Wright et al., 2019). There is a growing body of literature examining transcriptomic responses to thermal stress in corals (Meyer et al., 2011;Kenkel et al., 2013bKenkel et al., , 2014Anderson et al., 2016;Davies et al., 2016;Dixon et al., 2020); though few studies to date have examined transcriptomic variability in corals exposed to estuarine and riverine runoff [Aguilar et al., 2017[Aguilar et al., , 2019Edge et al., 2013;Mayfield et al., 2013; see nontranscriptomic responses reviewed by Fabricius (2005)], or synergistic effects of multiple stressors therein (Wright et al., 2019). St. Lucie Reef provides an ecologically relevant scenario to investigate the potential synergistic impacts of thermal stress and estuarine discharge on coral health. Through comparative analyses between two experiments, we aimed to identify common responses to thermal and water quality stressors, in order to provide insight into the fundamental metabolic processes that may enhance stress tolerance and coral resilience in marginal reef environments.

In situ Experiment
The scleractinian coral Montastraea cavernosa was chosen as the study species for both experiments as it is widely distributed throughout the study area, has annotated transcriptomic and genomic references (Kitchen et al., 2015) 1 , and has temporally stable Symbiodiniaceae communities (Cladocopium spp.) documented for the fate-tracked colonies in the in situ experiment (Klepac et al., 2015). Fifteen M. cavernosa colonies were repeatedly sampled over nine time points spanning periods of relatively low and high estuarine discharge in 2013-2016, resulting in a total of 92 coral tissue fragments collected (Table 1).
Variable water visibility and sea state conditions occasionally  Figure 1) with increasing distance from the St. Lucie Inlet, which has been shown to approximate a gradient of decreasing exposure to estuarine discharge water (Beal et al., 2012;Shatters, 2017). Cumulative estuarine discharge flow rates measured at the canals C-44, C-23, C-24, and Ten Mile Creek (Figure 1) were obtained from the South Florida Water Management District's DBHYDRO database ( Supplementary  Figure 1), where daily flow rates were averaged for 2 weeks leading up to the corresponding sampling time points. Daily rainfall data were obtained from the same database at the C-44 canal control structure. Small tissue fragments (∼5 cm 2 ) were collected from the perimeter of each colony using a hammer and chisel, transported to the surface in separate zip-top bags containing ambient seawater, and preserved in TRIzol reagent within 30 min of collection. Preserved samples were then kept on ice until storage at −80 • C prior to RNA extraction and library preparation.

Ex situ Experiment
Full details of the experimental apparatus used in the ex situ experiment are described in Shatters (2017). In summary, a two-factor (temperature and water) factorial apparatus was constructed using twelve 150-L glass aquaria across six open-top raceways. Three tank replicates per each of the four experimental treatments were randomly distributed among the raceways, with continuous monitoring of temperature, pH, and ORP with Neptune Systems Apex controllers, daily monitoring of salinity with a refractometer, and a 9-h photoperiod using Hydro Grow Sol LED lights set to comparable PAR levels (∼102 µmol m −2 s −1 ) found at the coral collection site based on in situ data collection with Onset HOBO data loggers. Temperature treatment levels (control: 25 • C; elevated: 30 • C) were determined from mean temperatures recorded before and during discharge events, respectively, by Odyssey conductivity and temperature data loggers deployed from October 2013 to January 2015 at the St. Lucie Reef study sites (Shatters, 2017; Figure 1 inset). Raceways were temperature-controlled using immersion heaters and in-line chillers to reflect the two temperature treatments, with three raceways per temperature treatment. Offshore water (∼32 psu) was collected via a subsurface pump from the Florida Oceanographic Society at high tide. The discharge water parameters were similarly determined from salinity data obtained before and during discharge events, where the mean salinity of the estuarine discharge plume during outgoing tides was ∼28 psu. Discharge water was collected via a pump at 0.5 m depth in the inlet during low slack tide, 5 days after flow rates through the canals from Lake Okeechobee reached a flow rate of 99.02 m 3 s −1 (Shatters, 2017). The resulting water (19 psu) was then adjusted with offshore water to achieve the treatment salinity of 28 psu, resulting in a ∼17/83% mix between discharge and offshore water, respectively. Both water types were initially filtered (50, 10, then 1 µm) and UV-sterilized, and stored in opaque holding tanks.
Due to state park regulations restricting whole-colony collection, eleven Montastraea cavernosa colonies of 150-175 cm 2 in size with no tissue loss, paling, or other visible signs of stress were collected south of SLR at Breakers Reef (26.6908 N, −80.0181 W) at depths of 13-17 m (Figure 1). Temperature data loggers deployed in June 2015 confirmed that ambient temperature was similar between sites (∼28 • C) despite difference in depth. Following transport to the lab, colonies were cut with a tile saw into twelve 3 cm × 4 cm colony replicates [genotyped in Dodge et al. (2020)], and randomly positioned in each of the twelve treatment tanks (n = 33 per treatment, n = 132 total; Table 2). The ex situ experiment was conducted over 32 days, with a 7-day recovery/acclimatization period, 5-day ramp to experimental conditions (achieved with daily ∼15% water changes and <1 • C temperature adjustments), and 20-day experimental treatment. Coral fragments were fed ∼10 g of brine shrimp per tank every 3-4 days over the course of the experiment. At the end of the treatment period, coral fragments were split into subsamples for gene expression and symbiont/chlorophyll analyses, and preserved in TRIzol/−80 • C and dry/−20 • C, respectively.

Symbiont and Chlorophyll Metrics
Symbiont and chlorophyll metrics, including Symbiodiniaceae density (cells cm −2 ), areal chlorophylls a and c 2 (chl cm −2 ), cellular chlorophylls a and c 2 (chl cell −1 ), and chlorophyll a:c 2 ratio, were quantified from the frozen subsamples as described in Shatters (2017). All statistical analyses were conducted in the R statistical environment (R Core Team, 2019). Symbiont and chlorophyll metrics could not be normalized with data transformation, therefore non-parametric tests were used. Data were first assessed for homogeneity of variance among treatment groups using the betadisper function of the package vegan (Oksanen et al., 2015), utilizing Euclidean distance matrices and 9,999 permutations. Following identification of statistically equal variance (all p > 0.05), univariate PERMANOVAs were conducted using the adonis function in vegan to identify potential effects of experimental treatments (genotype, temperature, water) on symbiont and chlorophyll metrics. Pairwise comparisons were made among single-factor combinations of temperature and water treatments for each of the metrics using the Total number of coral genotypes denoted by n total , with the number of genotypic replicates as n reps . The total number of samples for each treatment group given as n collected , with the number of samples used in statistical analyses as n ex situ .
nctp function of the package nparcomp (Konietschke et al., 2015). Multivariate effects were assessed using a three-factor PERMANOVA and pairwise colony tests using the packages vegan and RVAideMemoire (Hervé, 2019), respectively. Both multivariate tests used the same model parameters as indicated before, with FDR p value correction for the pairwise comparisons.

Library Preparation and Sequencing
For both experiments, RNA extractions and Tag-Seq library preparation protocols were followed as in Studivan and Voss (2020). In short, total RNA was extracted using a modified phenol-chloroform protocol Sacchi, 1987, 2006), DNase I treated using an Ambion kit, and purified using a LiCl precipitation. Only samples with at least 1 µg of total RNA were used for library preparation according to the Tag-Seq protocol first described in Meyer et al. (2011) and modified as in Studivan (2020b). In summary, cDNA libraries were amplified over 22 cycles, dual-barcoded with individual 5 and 3 Illumina indices, and equally combined using qPCR into three pools encompassing samples from both experiments. A total of 275 sample pools were generated (in situ: n = 111; ex situ: n = 164), including libraries that were removed for subsequent statistical analyses. Pools were then sequenced at the University of Wisconsin Biotechnology Center on an Illumina HiSeq 2500 with V4 chemistry and 15% PhiX spike-in, generating 1 × 50 bp single-end reads. For all samples across both experiments, 1.30 billion reads were produced, with a mean number of 4.73 ± 0.16 (mean ± SE) million reads per sample.

Differential Expression Analyses
Raw sequences were first processed separately for each experiment to obtain cleaned, mapped, and quantified reads using pipelines described in a GitHub repository associated with a previous study (Studivan, 2020b;Studivan and Voss, 2020). Cleaned reads were mapped to a concatenated reference transcriptome composed of the coral host M. cavernosa (Kitchen et al., 2015; updated genome-derived transcriptome from https://matzlab.weebly.com/data--code.html), and its dominant algal symbiont identified at SLR (Klepac et al., 2015; Supplementary Data 1), Cladocopium spp. (Davies et al., 2018). Creation and annotation of the combined coral host and algal symbiont reference transcriptome is described in a GitHub repository associated with a previous study (Studivan, 2020a;Studivan and Voss, 2020). Subset datasets were created to remove additional samples not part of the respective in situ and ex situ experimental designs. For the in situ experiment, 49 samples were retained (seven colonies with complete sampling coverage across seven time points) to achieve a balanced design with sufficient replication (Table 1). For the ex situ experiment, 130 samples were retained as part of the design ( Table 2). Mapping efficiency to the sequenced samples from both experiments was 17.0 ± 0.3%, resulting in 24,182 isogroups for the in situ experiment (15,712 host and 8,470 symbiont), and 27,933 isogroups for the ex situ experiment (18,210 host and 9,723 symbiont). Low mapping efficiency has been previously shown to be an issue with the updated version of the M. cavernosa transcriptome , however, improved mapping with the original transcriptome was countered by relatively poor GO and KOG functional annotations. Statistical analyses were therefore conducted using counts generated from the updated coral transcriptome, following the separation of count data into host and symbiont datasets for separate statistical analyses of coral and algal gene expression for each experiment. Datasets were pre-filtered to remove low-count genes (cumulative count <10 across all samples in a dataset) to decrease processing time, and then tested for outliers using the package arrayQualityMetrics (Kauffmann et al., 2009). Samples violating the distances between sample arrays (S a ) criterion for each of the respective datasets were removed, whereas three host and one symbiont samples were removed from the in situ experiment (Supplementary Figure 2), and three host and five symbiont samples from the ex situ experiment (Supplementary Figure 3). Count data were normalized using the variance stabilized transformation in the package DESeq2 (Love et al., 2014) for data visualization and PERMANOVAs. PERMANOVAs were conducted using vegan with 9,999 permutations for both experiments to identify if model terms were significantly affecting gene expression patterns in either the coral host or algal symbionts. Multivariate differences among samples were visualized using PCoAs of Manhattan distance in vegan and heatmaps of individual gene expression patterns using the package pheatmaps (Kolde, 2019). Tests of differential expression for both experiments were conducted using DESeq2 to identify differentially expressed genes (DEGs) associated with factors (model terms) and among factor levels (contrast statements) using the "reduced" flag on the respective main models. For the in situ experiment, a two-factor main model of time and site was used, with pairwise comparisons between sites and time points. For the ex situ experiment, a three-factor model of coral genotype, temperature, and water was used with interaction terms between the main model factors and a blocking factor (to assess whether tank location had an impact on gene expression patterns). Additionally, pairwise comparisons were made between single-factor levels of temperature and water treatments.

Gene Enrichment and Weighted Gene Co-expression Network Analyses
Differentially expressed genes were exported for each of the respective comparisons within experiments and analyzed for larger expression patterns across gene pathways using the package GO-MWU (Wright et al., 2015) for gene ontology (GO) databases, and the package KOGMWU (Dixon et al., 2015;Matz, 2016;Strader et al., 2016) for eukaryotic orthologous groups (KOGs). Both GO and KOG databases contain universal functional group annotations of individual genes based on biological and biochemical functions (Ashburner et al., 2000;Tatusov et al., 2000Tatusov et al., , 2003 The Gene Ontology Consortium, 2017). GO-MWU identifies GO categories that are significantly enriched by up-regulated or down-regulated genes using Mann-Whitney U tests. The resulting outputs identify significantly enriched GO categories and their relative levels of expression across the three GO divisions (Molecular Function, Biological Processes, and Cellular Components) according to an FDRcorrected p value cutoff of 0.05, or a fold change cutoff of 1. The KOGMWU analysis is similar, although enrichment of functional pathways is compared across universal KOG categories to produce a heatmap (Nordberg et al., 2014).
Prior to weighted gene co-expression network analysis (WGCNA) of the ex situ experiment, host and symbiont datasets were further filtered to remove genes with counts <10 across >90% of samples. Correlational networks were then created blind to experimental design based on gene expression patterns using the package WGCNA (Zhang and Horvath, 2005;Langfelder and Horvath, 2008). In short, groups of genes with similar expression patterns (modules) were first identified across samples, then correlated to experimental factors (coral genotypes, temperature/water treatments, symbiont, and chlorophyll metrics) based on strength of the module-trait associations. The protocol used to conduct WGCNA, as well as scripts for all other analyses described above, can be found in a GitHub repository associated with this study (Studivan, 2021). Soft thresholds were adjusted to fit each particular dataset (11 for host, 4 for symbiont) to achieve a suitable scale-free topology index, with no merging of original modules. Moduleassociated genes were then exported from the host dataset for GO enrichment analysis.

Differential Expression
Pre-filtering of raw counts produced datasets of 8,497 coral host and 1,652 symbiont transcripts across 46 and 48 samples, respectively. Analysis of gene expression patterns across in situ experiment factors (site, time) using PERMANOVA yielded that only time was a significant factor for both host and symbiont datasets (host: R 2 = 0.24537, p < 0.0005; symbiont: R 2 = 0.27559, p < 0.0001; Table 3). DESeq2 identified 927 host and 50 symbiont differentially expressed genes (DEGs) associated with time, while only 60 host and 30 symbiont DEGs were identified for site, and zero DEGs for the interaction term (Supplementary Table 1). DEGs identified for each of the pairwise comparisons across sites and sampling time points can also be found in Supplementary Table 1. PCoAs and heatmaps of samples through time suggested that transcriptomic variation was primarily attributed to differences from the first two time points (October 2013 and September 2014) to all subsequent time points (Figure 2). There appeared to be limited overlap in the PCoAs between samples collected during periods of relatively high estuarine discharge and rainfall (indicated by cooler colors in Figures 2A,C and Supplementary Figure 1), however, similarities between sample dispersions from consecutive time points (e.g., October 2013 and September 2014) suggest that other environmental factors beyond discharge rates alone may have influenced gene expression patterns in both the coral host and their algal symbionts.

Gene Enrichment Analyses
Since site did not have a significant impact on gene expression, subsequent analyses were conducted among sampling time points, combining samples across sites. Analysis of gene ontology (GO) categories enriched for the coral host dataset over sampling time points identified two main groups of pairwise time comparisons. From the first two time points (October 2013 and September 2014) to subsequent time points, all but one GO category in the biological process (BP) division demonstrated reduced expression in the later versus earlier time point (Supplementary Figure 4). Of these, several categories were consistently represented in pairwise time comparisons, including: "RNA catabolic process, " "nucleartranscribed mRNA catabolic process, nonsense-mediated decay, " "protein localization to endoplasmic reticulum, " and "proteincontaining complex disassembly." Between time points after October 2015, the majority of GO categories were instead overexpressed in later versus earlier time points, including several of the aforementioned categories (Supplementary Figure 4). Analysis of eukaryotic orthologous groups (KOGs) identified only one KOG class in the host dataset with relatively strong variation through time, where "nuclear structure" was negatively enriched between the first two time points and subsequent time points, and represented by higher-expressed genes between 2015 and 2016 time points (Supplementary Figure 5). Other KOG classes, including "cell motility, " "nucleotide transport and metabolism, " and "extracellular structures" demonstrated mild variation in expression through time. For the symbiont dataset, however, the majority of the KOG classes contained under-expressed genes in later versus earlier time points, with the exception of "extracellular structures" that showed variation through time. This KOG class was more highly-

Symbiont and Chlorophyll Metrics
Coral genotype and temperature treatments were significant factors affecting all symbiont and chlorophyll metrics, with the exception of cellular chlorophylls a and c 2 for temperature (Supplementary Table 2). Only areal chlorophyll a was marginally affected (p < 0.0612) by the water treatment; all other metrics were non-significant. Multivariate pairwise differences among coral genotypes are found in Supplementary  combinations of temperature and water treatments were only found for comparisons between control and elevated temperature treatments for symbiont density and both areal chlorophyll metrics. Symbiont density and areal chlorophylls a and c 2 were significantly reduced in both elevated temperature treatments, with no significant pairwise differences between offshore and discharge water treatments, respectively (Supplementary Figure 6). The ratio of chlorophylls a:c 2 and cellular chlorophylls a and c 2 did not vary significantly among treatments.

Differential Expression
Pre-filtering of raw counts resulted in 11,248 coral host and 1,834 symbiont transcripts across 127 and 125 samples, respectively. PERMANOVA results for the ex situ experiment indicated that genotype and temperature were significant factors, while water was marginally significant for the host dataset, and there was a significant temperature:water interaction (genotype: R 2 = 0.14197, p < 0.0001; temperature: R 2 = 0.01098, p < 0.0197; water: R 2 = 0.00984, p < 0.0516; temperature:water: R 2 = 0.01011, p < 0.0452;  (Figure 3). Corals appeared to respond more strongly to estuarine discharge than did the symbionts across all treatments. DEGs associated with each of the remaining interaction terms and pairwise comparisons between treatments are found in Supplementary Table 1. PCoAs demonstrated that there was stronger genotypic variation in the coral host than in the symbiont dataset, as well as evidence of a synergistic effect of temperature and water treatments for both datasets (Figure 4). Each of the temperature and water treatments alone were relatively close to one another in the ordination space, with a somewhat distinct sample spread for the combined treatment (elevated discharge).

Gene Enrichment and Weighted Gene Co-expression Network Analyses
Gene ontology categories across all three divisions (molecular function, MF; biological process, BP; cellular component, CC) were largely positively-enriched in elevated versus control temperature and discharge versus offshore water, respectively ( Figure 5). The overall number of significantly enriched GO terms, as well as the number of DEGs within those GO categories, were higher for the temperature treatment compared to the water treatment, although there were several similarities between the two. GO categories "structural constituent of ribosome" (MF), "translational initiation" (BP), "translational elongation" (BP), "RNA catabolic process" (BP), "cellular component disassembly" (BP), and several ribosomal pathways (CC) were elevated in both temperature and water treatments relative to their respective controls. Few GO categories, on the other hand, were underexpressed in treatments ["lysase" (MF), "transferase" (MF), "carboxylic acid catabolic process" (BP), "cell projection part" (CC), and "basolateral plasma membrane" (CC)], although none of these were represented across both temperature and water treatments. Eukaryotic orthologous group analyses across main factors identified three main trends: (1) expression patterns for genotype and temperature factors were strongly correlated for both host and symbiont datasets, (2) some KOG classes had an inverse relationship between responses to temperature and water stress in both datasets, and (3) corals and symbionts responded differently to the same factors (Supplementary Figure 7). Indeed, correlation coefficients between genotype and temperature  responses were r = 0.99 with highly significant p values in both datasets, suggesting that variation in response to temperature stress was innately tied to coral genotype, and was largely captured by the latter term in statistical models. Several KOG classes, including "cytoskeleton, " "extracellular structures, " and "nuclear structure" were differentially expressed between water and temperature treatments for the coral host (Supplementary Figure 7A). For the symbiont dataset, "extracellular structures" and "coenzyme transport and metabolism" were differentially expressed between water and temperature treatments, while other KOG classes showed generally weaker enrichment in the water treatment compared to the temperature and genotype factors (Supplementary Figure 7D). When examining pairwise combinations of treatments against the control (control offshore) for the host dataset, only two KOG classes ("nuclear structure" and "cell motility") were positively enriched in control discharge versus both elevated temperature treatments (Figure 3A). In the symbiont dataset, the same was true of "cell motility, " while "extracellular structures" and "nucleotide transport and metabolism" classes were more-expressed in control discharge and elevated offshore treatments versus elevated discharge samples ( Figure 3D). Both of these observations may imply that the coral host responded more strongly to elevated temperatures than estuarine discharge, while algal symbionts demonstrated a synergistic response to elevated temperatures and estuarine discharge combined.
Further filtering of raw gene counts for WGCNA resulted in a dataset of 967 host and 43 symbiont genes. WGCNA clustered genes into five modules for the host dataset and a single module for symbionts; therefore only host data are presented (Supplementary Figure 8; five host genes were unassigned). All of the modules demonstrated significant correlations to individual coral genotypes, while turquoise, yellow, green, and blue modules were also significantly correlated with experimental treatments, suggesting genotypedriven transcriptomic responses to stress. None of the symbiont or chlorophyll metrics were significantly correlated to any of the modules. Variation in module correlation directions between elevated offshore and elevated discharge treatments (turquoise, yellow, and green modules) corroborated a potential synergistic effect of temperature stress and water treatments on coral gene expression patterns. GO-MWU analysis of host modules identified significantly enriched GO categories within blue and brown modules only, with enriched ribosomal and cell division GO categories in the blue module, and immune/stress response and cell component GO terms enriched in the brown module (Supplementary Figure 8).

In situ Experiment
Temporal analyses of gene expression across three sites at St. Lucie Reef revealed that time was the only significant factor driving transcriptomic patterns for both the coral host and its symbiotic algae. This suggests that despite observations of rapid changes in water quality parameters following discharge events along a gradient with increasing distance from the St. Lucie Inlet (Beal et al., 2012;Pickering and Baker, 2015;Shatters, 2017), corals responded similarly at the transcriptomic level among sites. There did not appear to be a strong relationship between estuarine discharge rates at the time of sampling to either coral or symbiont expression profiles. Taken in the larger context of seasonal variation including rainfall and discharge rates, however, there were observed shifts in expression patterns from the first two time points to all subsequent time points (Figure 2). The wet and dry seasons in 2013-2014 followed the typical pattern described earlier (Beal et al., 2012;Klepac et al., 2015), but the 2015 "wet" season was unusually dry (Shatters, 2017), and was followed by an abnormally rainy 2015-2016 "dry" season (Stockley et al., 2018). Likewise, the number and magnitude of discharge events varied over the course of the study (Supplementary Figure 1). In 2013-2014, there were several larger pulses of discharge (>50 m 3 s −1 ) during the wet season, with relatively low or no discharge events during the dry season. Beginning in 2015, conversely, there were lower (<40 m 3 s −1 ), but consistent, discharges throughout the year, followed by approximately 9 months of discharge at rates in excess of 50 m 3 s −1 , and in excess of 150 m 3 s −1 in 2016 ( Supplementary Figure 1; Stockley et al., 2018).
The apparent shifts in gene expression patterns at later versus earlier time points is perhaps the result of covarying environmental factors through time, namely compounding rainfall and discharge flow rates that introduce hyposalinity, nutrient, and elevated temperature stressors to the reef. In the coral host, this is best illustrated by relative expression through time of genes within the "nuclear structure" KOG class. Despite relatively low sequencing coverage associated with the algal symbionts, genes in the "extracellular structures" class were differentially expressed over time (Supplementary Figure 5). Many of the gene ontology (GO) categories that were under-expressed in comparisons between 2015-2016 and initial 2013-2014 time points were then over-expressed in 2016 versus 2015 time points, suggesting variable expression of the same gene pathways in response to different levels of estuarine discharge exposure and hyposalinity stress associated with increased rainfall. Few studies have examined transcriptomic responses to sublethal exposures of estuarine water in natural environments, where differentially expressed stress response genes were linked to increased precipitation and runoff events [Beal et al., 2012;Edge et al., 2013; but see also Wright et al. (2019)].
A recent metanalysis, however, has suggested that transcriptomic responses may be conserved across various stressors, with diverging responses dependent on stress intensity. Dixon et al. (2020) identified a type "A" response to highintensity bleaching, heat, hyposalinity, or immune challenge stressors where cell death, immune response, and protein degradation pathways were generally up-regulated, and a type "B" response to lower-intensity stress exposures with an inverse relationship to type "A" expression levels. During the time period of our study, no phenotypic changes occurred in the fate-tracked colonies (i.e., bleaching or tissue mortality) to suggest that corals were exposed to high-intensity stress. Expression of gene pathways associated with cell growth [e.g., KOG terms "nuclear structure" in the coral host and "extracellular structures" in the symbionts (Supplementary Figure 5); GO terms "protein localization to endoplasmic reticulum" and several mitosis-related pathways in the coral host (Supplementary Figure 4)] were generally more expressed during later sampling time points when discharge rates and rainfall were higher, compared to earlier time points with less consistent freshwater inputs to the reef (Supplementary Figure 1). Immune GO pathways including "regulation of immune effector process" and "regulation of production of molecular mediator of immune response" were conversely downregulated during the later sampling time points (Supplementary  Figure 4). These patterns would therefore suggest that corals at St. Lucie Reef were exhibiting a more type "B" response following sublethal exposures to estuarine discharge and hyposaline stress.
While individual sublethal stressors may not cause significant disruptions to the coral-algal symbiosis, consecutive or concurrent stress events that are common in nearshore environments may form a "negative feedback loop" ultimately leading to colony mortality. In particular, exposure to salinity stress and/or sedimentation in estuarine discharge water has been shown to increase susceptibility of corals to diseases (Pollock et al., 2014;Shore-Maggio et al., 2018). While it was perhaps inevitable that SCTLD would reach coral populations in the northern extent of the reef tract, repeated exposure to estuarine discharge may have exacerbated individual corals' susceptibility to SCTLD infection and contributed to the spread of the disease across additional reefs. Testing this hypothesis requires further manipulative experiments applying multiple stressors experienced by nearshore reefs to observations of SCTLD transmission and progression. Transcriptomic approaches may also increase the ability to recognize early signs of infection following exposure to multiple sublethal stressors. Beyond synergistic effects on coral/symbiont physiology and health, estuarine discharge may have additional ramifications for the survival of coral populations in southeast Florida. For example, exposure to hyposalinity stress can severely limit coral fertilization success (Hédouin et al., 2015). The frequent hyposaline conditions (e.g., 27-28 psu commonly measured at St. Lucie Reef but as low as 19 psu for brief periods) in this region due to estuarine outflow, Lake Okeechobee discharges, and watershed runoff may contribute to the rare occurrence of reproduction in St. Lucie Reef corals (Beal et al., 2012).

Ex situ Experiment
In the ex situ factorial experiment, there was clear evidence of transcriptomic responses to thermal stress in both the coral host, and to a lesser extent, their algal symbionts (Supplementary Table 1 and Supplementary Figure 7). It is possible that reduced representation of symbiont transcripts in part explains the limited variation among treatments and lack of significantly enriched GO pathways in the symbiont dataset. In the coral host, several gene pathways were elevated following exposure to the increased temperature treatments, including GO categories for metabolic and catabolic processes, protein translation and localization, and translational elongation, which are part of the metabolic and photosynthetic responses to thermal stress (Murata et al., 2007;Reyes-Bermudez et al., 2009;Seneca and Palumbi, 2015;Davies et al., 2016;Studivan and Voss, 2020). Ribosomal expression was also higher in elevated temperature corals, which is likely associated with an attempt by the corals to return to homeostasis following exposure to stress (Kültz, 2005;DeSalvo et al., 2008;Kenkel et al., 2013b). Several of the GO categories described above, including "protein localization to endoplasmic reticulum, " "nuclear-transcribed mRNA catabolic process, " and protein-containing complex and cellular component disassembly, were positively expressed in the final time points of the in situ experiment (Supplementary  Figure 4). The KOG classes "nuclear structure" and "extracellular structures" were also differentially enriched for temperature and water treatments (Supplementary Figure 7), reinforcing the notion that in situ corals at St. Lucie Reef were likely responding to sublethal stress during 2015 and 2016 time points similarly to the artificial stress treatments applied in the ex situ experiment. These similarities also demonstrated that corals in the ex situ experiment were demonstrating a type "B" transcriptomic response to sublethal environmental stress (Dixon et al., 2020).
While exposure to estuarine discharge water alone did not significantly alter coral or symbiont gene expression patterns in the ex situ experiment, many of the same GO categories that were positively enriched with thermal stress were also enriched with exposure to discharge water ( Figure 5). This observation concurs with the hypothesis of a conserved stress response among broad taxonomic groups including corals (Wang et al., 2009;Dixon et al., 2020), with variation of genes or expression levels within a general environmental stress response that may be stressor-specific. For example, hyposaline conditions can lead to differential expression of heat shock protein and oxidative stress genes, among others (Seveso et al., 2013;Aguilar et al., 2019). Our data suggest that responses were similar in terms of gene function, but that thermal stress caused stronger responses in both corals and their symbionts. There was also a synergistic effect between temperature and water treatments that could not be explained by elevated temperature or estuarine discharge alone (Figures 3, 4 and Supplementary Figure 6). Previous studies have observed synergistic effects of multiple stressors including temperature, salinity, and nutrient factors, with proteomic indicators of oxidative stress (Dias et al., 2019b), bleaching and reduced growth/condition (Coles and Jokiel, 1978), and reduced photosynthetic performance and/or mortality (Li et al., 2009;Faxneld et al., 2010;Dias et al., 2019a). Exposure to multiple stressors may cumulatively result in more severe impacts than with individual type "B" stressors alone, and can lead to a type "A" stress response (Dixon et al., 2020). In addition to direct effects of multiple stressors on coral health, thermal stress has been shown to impede several species' ability to cope with enhanced sedimentation (Bessell-Browne et al., 2017), and causes reductions of surface mucopolysaccharide layers (i.e., where mucus is produced) in P. clivosa (Pratte and Richardson, 2014). Such synergistic relationships could have important ramifications for the health of corals in Florida, as nearshore reef environments are often exposed to several stressors simultaneously with tidal fluctuations (Yentsch et al., 2002;Lirman and Fong, 2007;Beal et al., 2012;Gregg, 2013).
Synergistic effects of discharge water and thermal stress treatments were not apparent in the symbiont and chlorophyll metrics; significant differences among treatments were attributed only to elevated temperatures (Supplementary Figure 6). Previous studies examining physiological responses to hyposaline conditions identified reductions in symbiont densities and photosynthetic efficiency (Kerswell and Jones, 2003;Downs et al., 2009;Gardner et al., 2016), although this response may be dependent on the coral species. Siderastrea radians, a common shallow-water species found extensively in nearshore reef environments in Florida, was not affected by hyposalinity stress until its predicted lower tolerance limit of 10 psu (Chartrand et al., 2009;Lirman and Manzello, 2009). Such limits have not been established for M. cavernosa, but it has been hypothesized that salinity and light attenuation are major limiting factors in the distribution of nearshore coral populations (Yentsch et al., 2002;True, 2012), and that inshore environments experience seasonal stress and recovery phases (Haslun et al., 2016). While our ex situ experiment did not examine as many salinity treatments as the aforementioned studies, significant differences were observed between similar salinity levels as tested here. Also, areal chlorophyll metrics differed among treatments whereas cellular chlorophyll metrics did not, implying that decreases in areal chlorophylls a and c 2 were largely driven by reduced symbiont density in elevated temperature treatments, rather than a change in the number of symbiont cells within coral symbiosomes.
Variation in host responses to thermal stress was highly linked with coral genotype (Figure 3 and Supplementary  Figure 6), suggesting a strong genetic component to thermal stress susceptibility/resilience. A similar relationship has been identified with several coral species (Baums et al., 2013;Kenkel et al., 2013a;Bay and Palumbi, 2014), and host transcriptomic patterns may also be influenced by symbiont genotypes as well (DeSalvo et al., 2010;McGinley et al., 2012;Rocker et al., 2012;Cunning and Baker, 2020). In our study, it appears that M. cavernosa genotypes demonstrated variable susceptibility to thermal stress, as genotype explained an order of magnitude more transcriptomic variation in the coral dataset compared to temperature and water treatments (14.2% versus 1.1% and 1.0%, respectively; Table 3).
Additionally, there was a strong correlation of KOG expression between genotype and temperature factors (Supplementary  Figure 7). Furthermore, in WGCNA, many of the significant correlations observed between experimental factors and modules of related genes were attributed to individual coral genotypes rather than to temperature or water treatments (Supplementary Figure 8).
In the symbiont dataset, there was also an influence of coral genotype on transcriptomic patterns, despite relatively homogeneous symbiont assemblages across all samples in the ex situ experiment (Supplementary Data  1). The marginally significant p value associated with the genotype factor was likely the direct result of a lack of variation in symbiont among coral genotypes. However, 8.9% of the transcriptomic variation was explained by coral genotype, versus 1.9% explained by the temperature treatment ( Table 3). The strong links between coral genotypes and responses to thermal stress in both host and symbiont gene expression patterns may be indicative of resilience among some individuals in the sampled coral populations. Nearshore reef environments in southeast Florida (including Breakers Reef) experience seasonal and freshwater discharge-related environmental fluctuations, notably including changes in temperature at episodic and chronic temporal scales (Yentsch et al., 2002;Beal et al., 2012;Gregg, 2013). An interesting hypothesis is that these corals are particularly resilient to thermal stress given repeated exposure to temperature variation, resulting in a conserved transcriptomic response to sublethal stress events. Additional research is therefore warranted to further elucidate potential relationships between coral-algal genotypes and susceptibility to multiple stressors in this species across populations in Florida and beyond (Cunning and Baker, 2020;Smith et al., 2020).

Conclusion
Few studies to date have examined coral responses to multiple stressors (e.g., reduced salinity, higher nutrients, sedimentation, increased light attenuation) associated with estuarine discharge in coastal zones. Through complementary ex situ and in situ experiments, this study demonstrated that corals and their algal symbionts responded to thermal and estuarine discharge stressors, and that synergistic impacts from these stressors could influence their transcriptomic responses. Compared to a recent metanalysis of general environmental stress responses that generalized two characteristic transcriptomic patterns (Dixon et al., 2020), our results suggest that M. cavernosa displays a type "B" response following exposure to sublethal stress in both their natural environments and following manipulative treatments. Cell growth and protein production/modification pathways were observed to be up-regulated while immune pathways were down-regulated following temperature and hyposaline stress. Nearshore coral populations at the northern extent of Florida's coral reef tract are commonly exposed to multiple stressors in the form of variable temperature, salinity, and water quality (Yentsch et al., 2002;Lirman and Fong, 2007;Beal et al., 2012;Gregg, 2013) that may combine to negatively affect coral health. Continued examination of marginal populations found in this region may not only identify the molecular mechanisms responsible for coral resilience, but may also identify resilient genotypes for conservation and restoration initiatives. Given that modern coral reef ecosystems, especially those found in coastal environments, are more frequently exposed to multiple stressors, management actions are needed to mitigate potentially synergistic impacts of stressors on the health and resilience of coral reef organisms.

DATA AVAILABILITY STATEMENT
Raw Tag-Seq sequences are available in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under BioProject PRJNA698305, Accession numbers SAMN17712326 through SAMN17712600. Datasets generated from this study, including analysis scripts, are available in a GitHub repository (Studivan, 2021). Pipelines used for transcriptome annotation and bioinformatics followed protocols in previously-released GitHub repositories (Studivan, 2020a,b).

AUTHOR CONTRIBUTIONS
JV, JB, AS, and DD secured funding for this research, designed and performed the research, experiments, and initial data analyses, with field and lab assistance from MS. MS analyzed the data and wrote the manuscript, with significant contributions and revisions by all co-authors. All authors contributed to the article and approved the submitted version. provided critical feedback and support on the design of this study. Coral tissue samples were collected from the St. Lucie Inlet Preserve State Park under a permit from the Florida State Park Service (FL State Park permit number 06261715). Coral colonies were collected from Breakers Reef by permission of FWC. Paul Wills, Richard Baptiste, Peter Stock, Gary Luisi, Susan Laramore, and Jay Adams were integral in the development and construction of the ex situ experimental apparatus. Dennis Hanisak and Kristen Davis from the Indian River Lagoon Observatory at Florida Atlantic University's Harbor Branch Oceanographic Institute, and Brittany Biber from the Florida Oceanographic Society facilitated collection of estuarine discharge water and seawater. Sequencing services were provided by the University of Wisconsin Biotechnology Center, with exceptional sequencing support from Jeremy Niece and Joshua Hyman, and high-performance computing services were provided by Research Computing Services at Florida Atlantic University. This is contribution 2290 from Harbor Branch Oceanographic Institute at Florida Atlantic University.