Variation in Coral Thermotolerance Across a Pollution Gradient Erodes as Coral Symbionts Shift to More Heat-Tolerant Genera

Phenotypic plasticity is one mechanism whereby species may cope with stressful environmental changes associated with climate change. Reef building corals present a good model for studying phenotypic plasticity because they have experienced rapid climate-driven declines in recent decades (within a single generation of many corals), often with differential survival among individuals during heat stress. Underlying differences in thermotolerance may be driven by differences in baseline levels of environmental stress, including pollution stress. To examine this possibility, acute heat stress experiments were conducted on Acropora hyacinthus from 10 sites around Tutuila, American Samoa with differing nutrient pollution impact. A threshold-based heat stress assay was conducted in 2014 and a ramp-hold based assay was conducted in 2019. Bleaching responses were measured by assessing color paling. Endosymbiont community composition was assessed at each site using quantitative PCR. RNA sequencing was used to compare differences in coral gene expression patterns prior to and during heat stress in 2019. In 2014, thermotolerance varied among sites, with polluted sites holding more thermotolerant corals. These differences in thermotolerance correlated with differences in symbiont communities, with higher proportions of heat-tolerant Durusdinium found in more polluted sites. By 2019, thermotolerance varied less among sites, with no clear trend by pollution level. This coincided with a shift toward Durusdinium across all sites, reducing symbiont community differences seen in 2014. While pollution and symbiont community no longer could explain variation in thermotolerance by 2019, gene expression patterns at baseline levels could be used to predict thermotolerance thresholds. These patterns suggest that the mechanisms underlying thermotolerance shifted between 2014 and 2019, though it is possible trends may have also been affected by methodological differences between heat stress assays. This study documents a shift in symbiont community over time and captures potential implications of that shift, including how it affects variation in thermotolerance among neighboring reefs. This work also highlights how gene expression patterns could help identify heat-tolerant corals in a future where most corals are dominated by Durusdinium and symbiont-driven thermotolerance has reached an upper limit.

Phenotypic plasticity is one mechanism whereby species may cope with stressful environmental changes associated with climate change. Reef building corals present a good model for studying phenotypic plasticity because they have experienced rapid climate-driven declines in recent decades (within a single generation of many corals), often with differential survival among individuals during heat stress. Underlying differences in thermotolerance may be driven by differences in baseline levels of environmental stress, including pollution stress. To examine this possibility, acute heat stress experiments were conducted on Acropora hyacinthus from 10 sites around Tutuila, American Samoa with differing nutrient pollution impact. A threshold-based heat stress assay was conducted in 2014 and a ramp-hold based assay was conducted in 2019. Bleaching responses were measured by assessing color paling. Endosymbiont community composition was assessed at each site using quantitative PCR. RNA sequencing was used to compare differences in coral gene expression patterns prior to and during heat stress in 2019. In 2014, thermotolerance varied among sites, with polluted sites holding more thermotolerant corals. These differences in thermotolerance correlated with differences in symbiont communities, with higher proportions of heattolerant Durusdinium found in more polluted sites. By 2019, thermotolerance varied less among sites, with no clear trend by pollution level. This coincided with a shift toward Durusdinium across all sites, reducing symbiont community differences seen in 2014. While pollution and symbiont community no longer could explain variation in thermotolerance by 2019, gene expression patterns at baseline levels could be used to predict thermotolerance thresholds. These patterns suggest that the mechanisms underlying thermotolerance shifted between 2014 and 2019, though it is possible trends may have also been affected by methodological differences between heat stress assays. This study documents a shift in symbiont community over time and captures potential INTRODUCTION Anthropogenic climate change presents a bleak future for many species with major declines in global biodiversity projected (Bellard et al., 2012). Rising global temperatures and other changing environmental conditions are predicted to push many species past their physiological tolerance limits (Somero, 2010). An important, yet often overlooked component of projecting the effects of climate change on species persistence is estimating both a taxa's thermal threshold for physiological stress, and its capacity for acclimatization (Seebacher et al., 2015).
Reef building corals present a good model for studying withingeneration acclimatization processes (i.e., phenotypic plasticity) because they are sessile and long-lived, and have experienced rapid climate-driven declines in the past two decades, often with differential survival among individuals (Marshall and Baird, 2000;West and Salm, 2003;Hughes et al., 2017). These differences may be driven by microenvironmental differences, genetics, or phenotypic plasticity. Since corals have high specieslevel and individual-level differences in thermotolerance (e.g., Fuller et al., 2020;Loya et al., 2001), they are good candidates for studies investigating plastic responses to climate change.
Acclimatization processes can occur at multiple levels within a coral "holobiont." A coral holobiont is the coral animal plus all of the associated microorganisms and symbiotic algae that live on and within the coral tissue (Rohwer et al., 2002). Plasticity can occur within the coral animal itself (e.g., gene expression shifts to increase heat shock proteins or antioxidants during thermal stress; Dixon et al., 2020), or within the members of the coral holobiont community that contribute to coral thermotolerance (e.g., endosymbiont shifts toward more heattolerant species or "symbiont shuffling"; Berkelmans and van Oppen, 2006). Symbiont shuffling may increase thermotolerance when there is a change in the proportions of different genera within the Symbiodiniaceae family, often seen as an increase in the proportion of heat-tolerant Durusdinium (Berkelmans and van Oppen, 2006;Stat and Gates, 2011;Ladner et al., 2012;LaJeunesse et al., 2018). Gene expression shifts and symbiont shuffling may increase thermotolerance, but there are limits to that increase. Once symbiont shuffling has occurred and corals host entirely one species/strain of symbiont, they may have maximized their ability to improve their symbiont-driven thermotolerance (Howells et al., 2020). While other plasticity processes exist (e.g., transgenerational plasticity; Putnam and Gates, 2015), gene expression shifts and symbiont shuffling are well-studied mechanisms by which corals are known to adjust their thermotolerance within the lifetime of an individual (Berkelmans and van Oppen, 2006;Bellantuono et al., 2012) and are the focus of this study.
While gene expression shifts and symbiont shuffling have been examined in response to temperature stress, fewer studies have investigated how these processes respond when temperature stress interacts with local stressors such as pollution. Landbased pollution can carry nutrients and sediments into marine environments, affecting corals and their symbionts (Silbiger et al., 2018). There have been multiple accounts of elevated nutrients lowering thermotolerance (Wooldridge, 2009;Wiedenmann et al., 2013;Donovan et al., 2020) or conversely improving thermotolerance (Béraud et al., 2013;Zhou et al., 2017;Riegl et al., 2019;Becker et al., 2021). These discrepancies may be partially explained by differences in nutrient type or ratios among nutrient concentrations (Morris et al., 2019;Burkepile et al., 2020). While the effects of nutrient pollution have been investigated in lab-controlled studies, there are fewer instances of field-based assessments of pollution and thermotolerance (Wooldridge, 2009;Becker and Silbiger, 2020;Becker et al., 2021). Recent fieldbased studies of how pollution impacts thermotolerance have not explored underlying mechanisms such as symbiont shifts and coral host gene expression (Becker and Silbiger, 2020;Becker et al., 2021).
Long-term pollution stress could increase or decrease coral thermotolerance to acute heat stress. When corals are exposed to multiple stressors, many studies have shown synergistic effects, whereby the cumulative effect of two stressors is greater than each stressor independently (Kersting et al., 2015;Towle et al., 2016;Ellis et al., 2019). However, other studies have shown that multiple stressors can produce antagonistic effects in corals, whereby the cumulative effect of two stressors is less than that of each stressor alone (Zhou et al., 2017;Marangoni et al., 2019;Darling et al., 2020). Since the cellular stress response (CSR) is similar even for different types of environmental stresses (e.g., Dixon et al., 2020;Kültz, 2020), it is possible that if mild long-term pollution stress triggers macromolecular damage and increases constitutive expression of stress response genes, this could lead to higher thermotolerance during heat stress. This may be due to the "frontloading" of stress response genes in polluted waters, whereby the baseline gene expression more closely resembles CSR gene expression, better preparing the coral to tolerate acute stress (Barshis et al., 2013;Thomas et al., 2018). It is also possible that long-term pollution stress has induced symbiont community shifts in favor of more stress-tolerant symbionts, leading to higher thermotolerance (LaJeunesse et al., 2010). This concept of increased tolerance to one stressor due to exposure to a different stressor is termed "cross-tolerance" FIGURE 1 | Topographic map of Tutuila, American Samoa of Tutuila, American Samoa with study site locations labeled. Green represents low pollution, yellow represents moderate pollution, and red represents high pollution. Pollution designations were determined using human population and nutrient loading data (DiDonato, 2004;Sudek and Lawrence, 2017;Shuler et al., 2019;Shuler and Comeros-Raynal, 2020). Unfilled circles were sampled only in 2014 and filled in circles were sampled in 2014 and 2019. and has been demonstrated in many species (Li and Hahn, 1978;Sabehat et al., 1998;Ely et al., 2014;Gunderson et al., 2016).
In this study, we investigate thermotolerance of a common reef-building coral, Acropora hyacinthus, on reefs of differing pollution levels around the island of Tutuila, American Samoa in 2 years, 5 years apart (2014 and 2019). During this 5 year period, mass bleaching events occurred in 4 of the 5 years on American Samoan coral reefs (Supplementary Figure 1, Witze, 2015;Sudek and Lawrence, 2017;Morikawa and Palumbi, 2019). We also measure two underlying mechanisms of plasticity that may affect thermotolerance: endosymbiont community composition and gene expression shifts. We investigate the role of pollution in triggering these plastic mechanisms and their effect on coral thermotolerance.

Site Selection and Field Collections
Ten sites of differing pollution levels were chosen to represent a gradient of pollution around Tutuila, American Samoa (Figure 1). Sites were categorized as either "low, " "moderate, " or "high" pollution based on human population and nutrient loading data, following established designations (DiDonato, 2004;Sudek and Lawrence, 2017;Shuler et al., 2019;Shuler and Comeros-Raynal, 2020). In March and July 2014, ten colonies of Acropora hyacinthus were sampled at each of ten sites between 07:00 and 11:00 (Figure 1). Eight fragments of approximately 2 cm 3 were collected haphazardly from each colony, with two replicates processed for genetic sampling (stored in RNA preservation buffer). The remaining six fragments were transported back to our experimental facility (in Vaitogi in March, and Coconut Point in July), and held at 28 • C overnight, with the assay to begin at sunrise the following day. Colonies were sampled at least 10 m apart to reduce the likelihood of sampling clones. All fragments were collected on snorkel and were sampled between 0 and 2 m in depth.
In August 2019, eight colonies of A. hyacinthus were sampled at a subset of five of the original 2014 sites between 07:00 and 11:00 (Figure 1). Seventeen fragments of approximately 2 cm 3 were collected haphazardly using stainless steel coral cutters from each colony: four replicates for each of four temperature treatments and one field control which was never placed in the temperature stress assay. Colonies were sampled at least 10 m apart to reduce the likelihood of sampling clones. All fragments were collected on snorkel and were sampled between 0 and 2 m in depth.

Threshold-Based Heat Stress Assay
Temperature dependent mortality in 2014 was measured using a 3-day threshold-based heat stress assay consisting of a ramp tank and a control. Replicate coral branches (n = 60 per tank, with three replicates of each of 10 colonies from two sites) were allowed to acclimate to tank conditions at 28 • C for 16 h. The temperature stress assay began at 06:00, in which the control tank was held at 28 • C for the duration of the stress, and the heated tank was brought to 30 • C over the course of the first 30 min, and then continuously ramped at 2 • C per 12 h during "daylight hours" (06:00-18:00), held overnight at the last ramp value, and then ramped between 06:00 and 18:00 for the following two successive days (Supplementary Figure 2). During daylight hours (06:00-18:00), full spectrum lights were measured using a planar Apogee PAR sensor and maintained at 500 µmol m −2 s −1 . Partial water changes were performed two times per day over the course of the assay. Water temperatures were controlled using a custom-built Arduino controller linked to aquarium heaters (Finnex HMA-200S Titanium Aquarium Heater) and Nova Tec IceProbe thermoelectric chillers (Nova Tec Products, San Rafael, CA). Water temperatures were measured every minute using HOBO UA-002-64 temperature and light loggers. Temperatures were kept to within 0.5 • C of the desired temperature.
Coral color paling was measured every 4 h between 06:00 and 22:00 using the CoralWatch R Coral Health Chart (Siebeck et al., 2006) and was recorded over the course of the assay until fragments experienced mortality, as measured by tissue sloughing. A fragment's thermal threshold of "degree heating days" was calculated as the number of degrees above the local bleaching threshold of 29 • C, as reported by NOAA's Coral Reef Watch, multiplied by the number of elapsed days above that threshold until mortality was reached. Data from colonies with controls that experienced mortality throughout the course of the experiment were discarded.

Ramp-Hold Heat Stress Assay
Thermal resistance to temperature stress in 2019 was measured using a standardized short-term acute heat stress assay and the Coral Bleaching Autonomous Stress System (CBASS; Voolstra et al., 2020). This approach has been shown to determine differences in coral thermotolerance similarly to a classic longterm heat stress assay Evensen et al., 2021). We chose to use this 1-day ramp-hold assay in 2019 rather than repeat the 2014-style assay due to the logistical advantage of a 1-day assay and to be more comparable to recent studies using CBASS Cunning et al., 2021;Evensen et al., 2021;Savary et al., 2021). This assay consists of four replicate tanks to test three experimental temperature treatments and one control, for a total of eight tanks. The temperature stress assay was conducted at 13:00 for each study site and continued until the following morning at 06:00 (Supplementary Figure 2). Replicate coral branches (n = 16 per tank, with two replicates of each of 8 colonies) were allowed to acclimate to tank conditions at 28 • C for 1 h. At 14:00, corals were exposed to ramp-hold: control (28 • C), 33 • C, 34 • C or 35 • C during a 2 h ramp, 3 h hold, 1 h ramp down to 28 • C, and overnight recovery at 28 • C. Light levels were measured using an Apogee underwater quantum meter twice during the assay and maintained between 210 and 250 µmol m −2 s −1 . To mimic natural field conditions, lights were turned off at 19:00 and turned on in the morning at 06:00 (Roleadro LED Aquarium Light). Partial water changes (∼2 L) using water from the sampling site were performed 4-5 times over the course of the assay. Water temperatures were controlled using a custom-built Arduino controller linked to aquarium heaters (Finnex HMA-200S Titanium Aquarium Heater) and custom-built cooling loops hooked to a Hamilton Technology Aqua Euro Max Aquarium Chiller. Water temperatures were measured every minute using HOBO UA-002-64 temperature and light loggers. Temperatures were kept to within 0.5 • C of the desired temperature.
Thermotolerance was measured through changes in visual color paling using the CoralWatch R Coral Health Chart (Siebeck et al., 2006) with the same observer for all trials and colorimetric analysis using an Olympus TG-5 taken by the same photographer in the same location for all measurements (sensu Winters et al., 2009). A subset of coral fragments (n = 4 per temperature treatment) were collected during heat stress (at 19:00) for RNAseq analysis. Samples were preserved in an RNA preservation buffer and stored at -20 • C until analysis, though only the 28 and 35 • C treatments were analyzed.
CoralWatch R Color Card and red pixel intensity measurements were used to generate a single metric for thermotolerance for each colony. Raw CoralWatch and red intensity values were normalized in an open interval from 0 to 1 for each temperature treatment from 28 to 35 • C. Logistic curves were fit to the data across temperatures and the midpoint of the curves was used as an indication of temperature at which bleaching occurred, with a maximum of 36.5 • C if curves did not reach the midpoint by 35 • C (using the nlme function in R; Pinheiro et al., 2021). The mean of the two midpoints (CoralWatch and red intensity) was used to generate a two-variable mean threshold, in units of • C. This two-variable mean was used to determine the most and least thermotolerant corals. The highest and lowest 10, 20, and 30% of the two-variable means were used to classify the 10, 20, and 30% least and most thermotolerant corals. All data and code are available at: https://github.com/melissanaugle/SCLERA_ Tutuila_Thermotolerance. All analyses were conducted in R version 4.1.0 (R Core Team, 2021).

Quantification of Symbiont Communities
For symbiont analysis, ten replicates per site were used in 2014 and eight replicates per site were used in 2019. Total genomic DNA was extracted from preserved fragments that had never entered the heat stress assay using a Qiagen DNeasy R Blood and Tissue Kit. Samples were prepared by removing excess buffer and homogenizing in a TissueLyser LT for 5 min at 50 Hz. Total DNA was measured using a NanoDrop spectrophotometer and a Qubit fluorometer. All DNA quantifications met the following criteria: > 2 ng/µl (on Qubit), 260:280 > 1.8, 260:230 > 1.59. Total DNA was prepared for qPCR using methods described in Cunning and Baker (2013) to quantify symbiont communities at each study site. Samples were run in duplicate in 2014 and triplicate in 2019, with a no-template control on a Biorad CFX96 Touch Real-Time PCR Detection System. Reaction volumes were 10 µl, with 5 µl Taqman Genotyping Master Mix and 1 µl genomic DNA template. qPCR analysis uses speciesspecific tags to identify Cladocopium and Durusdinium, which comprise the majority of the Symbiodiniaceae in American Samoan Acropora hyacinthus (Ladner et al., 2012). Since Durusdinium are more thermally tolerant than Cladocopium, ratios between the two genera provide a metric to understand coral thermotolerance contributed by the symbiont community (Cunning and Baker, 2013).
Ratios of Cladocopium to Durusdinium cycle threshold (Ct) values were calculated using results from qPCR. Baseline thresholds were chosen for each run to remove background noise. Ct values were recorded for samples that amplified past the threshold in fewer than 40 cycles. Ct values were averaged across duplicates in 2014 and triplicates in 2019 and cell numbers of Cladocopium and Durusdinium were calculated using the following formula: 2 (40−Ct)/ cell copy number (where cell copy number = 9 for Cladocopium and 1 for Durusdinium (Cunning and Baker, 2013). Log base 2 Cladocopium to Durusdinium ratios were calculated using the following formula: log2(cell number Cladocopium/cell number Durusdinium). All data and code are available at: https://github.com/melissanaugle/SCLERA_ Tutuila_Thermotolerance. All analyses were conducted in R version 4.1.0.

RNA Sequencing
Due to budget constraints, only the control and highest (28 • C and 35 • C) temperature treatments were used for mRNA sequencing and gene expression analysis. RNA extraction, cDNA library construction and sequencing were performed in two separate batches: once for Coconut Point, Faga'alu, and Faga'tele in February-March 2020 and again for Cannery and Vatia in March-April 2021. Total RNA was extracted using a Qiagen RNeasy R Plus Mini Kit, following the manufacturer's protocol. Samples were homogenized using a TissueLyser LT for 10 min at 50 Hz in 2020 and for 3 min at 50 Hz in 2021. RNA was assessed using a NanoDrop TM (Thermo Fisher Scientific), Qubit fluorometer and on a BioAnalyzer. 38 cDNA libraries were constructed from 300 ng of total RNA using the NEBNext R Ultra II RNA Library Prep Kit for Illumina R (E7530) with the NEBNext R Poly(A) mRNA Magnetic Isolation Module (E7490). Paired-end libraries were sequenced on a single HiSeq4000 150 bp lane at NovoGene in Davis, CA. Raw data obtained from sequencing was submitted to the NCBI Sequence Read Archive (SRA) database (Bioproject accession: PRJNA762371; SRA accession: SRP339664).

Differential Gene Expression Analysis
Low quality reads and adapter sequences were discarded using Trimmomatic (Bolger et al., 2014). Trimmomatic parameters were set to remove reads below 25 bp long, leading and trailing bases below quality "5, " and reads that did not meet quality standards for a sliding window where in a four base sliding window, the average quality per base drops below a 5. Sequences were also trimmed of adapter sequences including standard Illumina adapters and polyT sequences. Quality of trimmed reads was assessed using FastQC. 1 Sequenced reads were aligned to the reference A. hyacinthus transcriptome described in Barshis et al. (2013) using bowtie2 and counted using 1 http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ RSEM in UNIX (Li and Dewey, 2011;Langmead and Salzberg, 2012). Out of 33,496 possible contigs, 18,387 were discarded due to low expression, following methods described in Chen et al. (2016), and 15,109 contigs were retained and used for differential gene expression analysis. Gene expression statistical analyses were conducted using edgeR with a FDR of 0.001 and a fourfold change cut-off. A weighted gene co-expression network analysis (WGCNA) was conducted to compare coregulated gene networks and their association with pollution level, symbiont community, and thermotolerance (Langfelder and Horvath, 2008). WGCNA identifies co-expressed gene modules using hierarchical clustering of expression data and relates those modules to sample traits. WGCNA was performed on baseline gene expression from corals at control conditions, and separately on gene expression from corals under heat stress. Gene annotations were acquired from Barshis et al. (2013), and significant WGCNA modules (p < 0.05) were analyzed for enrichment of Biological Processes (BP) Gene Ontology (GO) terms using the GO_MWU R script (Wright et al., 2015). The GO_MWU R script uses a Mann-Whitney U-test and tests the kME (module membership score or eigengene-base connectivity) in among-module genes compared to other genes in the transcriptome outside the module to test if genes in the module of interest are significantly enriched (Wright et al., 2015;Huerta-Cepas et al., 2017).

Linear Mixed-Effects Models of Thermotolerance
Linear mixed-effects models were used to investigate the role of multiple variables in setting thermotolerance in 2014 and We built and present 5 mixed models, each using colonylevel features (e.g., tolerance, symbionts), site-level features (e.g., pollution) or temporal data (e.g., season, year). Each model held Site as a random factor. We specifically model: the correlation of pollution level and season vs. thermal tolerance in 2014, the correlation of pollution and Cladicopium:Durusdinium ratio vs. thermal tolerance in 2019; the correlation of year and pollution level vs. Cladicopium:Durusdinium ratio; the correlation of pollution level, season, and Cladicopium:Durusdinium ratio vs. thermal tolerance in 2014, and the correlation of three eigengene expression modules vs. thermal tolerance in 2019.
To explore the potential for non-linearity, models were first built using a generalized additive model approach using the mgcv package in R (Wood, 2011), but due to minimal evidence of nonlinearity, subsequent linear mixed-effects models were built using the lme4 package in R (Bates et al., 2015). Model assumptions were checked using the ggplot2 package in R (Wickham, 2016). Separate models were built for each year when investigating thermotolerance due to the differences in experimental design and thermotolerance metric between 2014 and 2019. To determine the importance of individual variables within each final model, we generated model sets from the final model using all-subsets regression (dredge in the package MuMIn), and then small sample size corrected Akaike information criterion (AICc) weight across each parameter (using importance in the MuMIn package) to generate a relative importance metric, ranging from 0 to 1 (Kamil, 2020).
In 2014, season, log2 ratio of Cladocopium:Durusdinium, and pollution level were included as explanatory variables of thermotolerance, measured as mortality in degree days, with site was included as a random factor. Model selection was based on AIC.
In 2019, log2 ratio of Cladocopium:Durusdinium, pollution level, and eigengene expression of significant WGCNA modules at baseline conditions (pre-heat stress) were included as explanatory variables of thermotolerance, measured as a two-variable mean value that represented thermotolerance. All data and code are available at: https://github.com/ melissanaugle/SCLERA_Tutuila_Thermotolerance. All analyses were conducted in R version 4.1.0.

Thermotolerance Patterns in 2014 and 2019
We compared colony-level coral thermotolerance in March and July 2014 among sites, expressed as time to mortality in degree heating days (Figures 2A,B). Thermotolerance varied among sites, with the most thermotolerant corals found at Cannery (high pollution) in March 2014 and Cannery (high pollution), Coconut Point (high pollution), Faga'alu (moderate pollution) and Faga'itua (moderate pollution) in July 2014. The least thermotolerant corals were found at Faga'alu (moderate pollution) and Faga'itua (moderate pollution) in March 2014 and at Amouli (moderate pollution), Aoa (low pollution) and Vatia (moderate pollution) in July 2014. In 2019, thermotolerance We show that thermotolerance was higher both in the hotter season (i.e., March) and at sites ranked as high pollution. Marginal R 2 = 0.520. was measured as a two-variable mean and was less variable among sites, with the most thermotolerant corals found at Vatia (moderate pollution) and the least thermotolerant corals found at Faga'tele (low pollution; Supplementary Figure 3 and Figure 2C).
Sites were pooled into their pollution level categories to examine trends between pollution and thermotolerance. In linear mixed models for 2014, holding site as a random effect, we show that thermotolerance was higher both in the hotter season (i.e., March) and at sites ranked as "High" for pollution level (Table 1 and Figures 2A,B; N = 115, Season p << 0.001, Pollution p < 0.05). Pollution level held no explanatory power in 2019 ( Table 2).

Shifts in Symbiont Community Between 2014 and 2019
In both 2014 and 2019, the log2 ratio of Cladocopium to Durusdinium favored Cladocopium at lower pollution levels and favored Durusdinium at higher pollution levels (Figure 3 and Table 3; N = 154, Pollution p < 0.05). Yet, 2014 and 2019 differed in their log2 ratio of Cladocopium to Durusdinium, with 2019 showing a stronger preference for Durusdinium compared to that seen in 2014 (Figure 3 and Table 3; N = 154, Year p < 0.05). Model estimates of Cladocopium to Durusdinium ratios controlling for pollution level showed a 16.7-fold decrease in 2019 compared to 2014 (i.e., mean log2 ratio reduced by 4.07 between years), suggesting a substantial and significant shift toward Durusdinium. Log2 Cladocopium to Durusdinium ratios also varied by site, with most sites across time points including corals hosting a combination of Cladocopium and Durusdinium, though both high pollution sites in 2019 hosted 100% Durusdinium (Supplementary Figure 4).

Gene Expression Responses
Differential gene expression analysis was performed on coral fragments from 2019 from the five sample sites and from two treatments (heat stress at 35 • C and control at 28 • C).
Gene expression profiles were driven strongly by temperature treatment, with 6,020 genes differentially expressed between heat and control treatments across all sites (edgeR, FDR < 0.001; Supplementary Figure 5). A weighted gene co-expression network analysis (WGCNA) was conducted to investigate how gene networks (called modules) from corals at control conditions (i.e., baseline expression patterns) correlated to pollution level, presence or absence of Cladocopium, and thermotolerance (Figure 4). One module, module C9, correlated to pollution level. No modules correlated to symbiont type. However, multiple gene modules showed correlations with thermotolerance under control conditions. Two modules (C5 and C8) correlated with the top 10% most thermotolerant coral colonies and two modules (C6 and C7) correlated with the top 20-30% thermotolerant coral colonies. Two additional modules (C3 and C16) correlated with the bottom 10-20% thermotolerant coral colonies.
Gene Ontology (GO) analysis was performed on the four modules relating to the most thermotolerant corals: modules C5, C6, C7, and C8. These modules contained 496, 267, 406, 117 genes, respectively. Of these four modules, only C8 showed any significant enrichment for Biological Processes (BP) GO terms. Module C8 was enriched for genes involved primarily in immune response, cytokine production, and multi-organism processes (p < 0.05, Supplementary Table 1).
Gene Ontology (GO) analysis was performed on modules C3 and C16 to determine which gene expression patterns at control conditions correlated with poor performance under heat stress. Module C3 contained 69 genes and module C16 contained 280 genes. Module C3 contained overrepresented GO categories relating to ion transport (p < 0.05, Supplementary Table 1). Module C16 contained many overrepresented GO categories, including negative regulation of protein catabolic process,  negative regulation of proteasomal protein catabolic process, regulation of ubiquitin-dependent protein catabolic process, regulation of protein catabolic process, intrinsic apoptotic signaling pathway by p53 class mediator, tail-anchored membrane protein insertion into endoplasmic reticulum (ER) membrane, establishment of protein localization to membrane, negative regulation of neuron differentiation, intrinsic apoptotic signaling pathway in response to ER stress, and regulation of cellular protein catabolic process (p < 0.05, Supplementary Table 1). These overrepresented categories were contained within broader, parent categories of protein catabolic process regulation, protein localization, apoptotic signaling, negative regulation of nervous system development, and mRNA processing.
A WGCNA was also performed on gene expression on corals under heat stress at 35 • C to determine how gene expression under heat stress varies by pollution level, presence or absence of Cladocopium, and thermotolerance (Supplementary Figure 6). Module H5 related to high pollution; modules H3 and H4 related to low pollution. Four modules (H3, H4, H13, and H14) showed positive expression in corals hosting both Cladocopium and Durusdinium and negative expression in corals hosting only Durusdinium. Gene Ontology (GO) analysis for the four modules correlated to symbiont community type showed enrichment for calcium ion transport, DNA repair, RNA processing, and developmental processes, among other GO terms (Supplementary Table 2). Four modules (H7, H8, H13, and H14) correlated with the 10-20% most thermotolerant corals. GO analysis for the four modules correlated to high thermotolerance showed enrichment for RNA processing, gene silencing by RNA, RNA splicing, and membrane docking, among other GO terms (Supplementary Table 2). Two modules (H12 and H15) correlated with the 10-20% least thermotolerant corals. GO analysis for the two modules correlated to low thermotolerance showed enrichment for MAPK activity and vesicle-mediated transport, which were upregulated in the least thermotolerant corals (Supplementary Table 2).

Predictors of Thermotolerance in 2014 and 2019
Linear mixed-effects models were used to identify the factors contributing to thermotolerance in 2014 and 2019. In 2014, we investigated the role of symbiont community, season, and pollution level on setting thermotolerance. Log2 Cladocopium to Durusdinium ratio, season and pollution level were all retained in a linear mixed-effects model as explanatory variables of thermotolerance, as measured in mortality in degree days ( Table 4). Although Cladocopium to Durusdinium ratio is correlated with pollution level, pollution level adds additional explanatory power to the model, and its inclusion significantly improves the model fit and survives model simplification using likelihood ratio tests. Taken together, the model including symbiont genotype, season, and pollution level explained 46.9% of variation in thermotolerance in 2014 (marginal R 2 = 0.469; Table 4, N = 81). Using a model set approach to quantify the relative importance of each variable, we found that season was the strongest predictor of thermotolerance (sum of AICc weights = 1.00), followed closely by pollution level (0.92), and log2 Cladocopium to Durusdinium ratio (0.12).
Applying the same model structure to the data from 2019, however, results in a fit with almost no explanatory power. In 2019, we investigated the role of symbiont community, pollution level, and baseline gene expression on setting thermotolerance. All 2019 sampling occurred in the same season, so that parameter was invariant in 2019. Yet, neither log2 Cladocopium to Durusdinium ratio nor pollution category was significantly correlated with the thermal tolerance phenotype in 2019, with the model explaining only 3.71% of variation in thermotolerance (marginal R 2 = 0.0371%, Table 2; N = 37). Gene expression, however, showed substantive and significant power to predict thermotolerance phenotypes, even for a modest sample size (i.e., N = 18). After model simplification, eigengene expression values of three of the six significant WGCNA modules were retained in the model: modules C3, C5, and C8 (Table 5), though models retaining module C7 were competitive. While including module C7 and other modules improved model predictive

DISCUSSION
In this study we measured coral thermotolerance and symbiont communities across a gradient of pollution in 2014 and 2019. In 2014, we found variation in thermotolerance by pollution level, largely driven by symbiont community differences. By 2019, we found no variation in thermotolerance by pollution level, likely due to shifts in symbiont communities over time toward Durusdinium, which reduced variation in symbiont communities among sites. Instead, thermotolerance in 2019 was best predicted by gene expression of a select number of modules during our control (i.e., baseline) treatment. Our results suggest that baseline gene expression may act as an important indicator of thermotolerance variation and may be especially important after shifts to more ecologically homogeneous symbiont communities, as we saw with Durusdinium.

Thermotolerance Varied by Pollution Level in 2014, Driven in Part by Symbiont Community
In 2014, we found differences among sites in their degreeday mortality as an indicator of thermotolerance, with more polluted sites holding more thermotolerant corals. Variation in thermotolerance was best explained by season, followed by variation in pollution level and symbiont community. Thermotolerance was higher in March, at the end of austral summer, compared to in July, in the middle of austral winter. This finding is consistent with past work showing that thermotolerance fluctuates by season and tends to be highest during warmer months (Jurriaans and Hoogenboom, 2020). In addition to season, symbiont community was important in driving differences in thermotolerance in 2014. Symbiont communities varied between the two time points in 2014, with a higher proportion of heat-tolerant Durusdinium in March compared to July (Figure 3). It is well established that symbiont communities within corals fluctuate by season and corals tend to host a higher proportion of thermotolerant symbionts in warmer months (Fitt et al., 2000;Chen et al., 2005;Thornhill et al., 2006).

Symbiont Communities Shifted Toward Durusdinium by 2019
At all sites over both years, the proportion of Cladocopium was greater at sites with lower pollution levels. This pattern of higher proportions of Durusdinium at high pollution sites aligns with previous work showing that corals living in more variable or stressful regions tend to favor Durusdinium (Fabricius et al., 2004;Oliver and Palumbi, 2009;Carballo-Bolaños et al., 2019). Some work has also linked increased proportions of Durusdinium to areas with higher pollution level or human impact (LaJeunesse et al., 2010). Our findings support the idea that corals exposed to chronic pollution stress, similarly to heat stress, favor Durusdinium. In 2014, this pattern was most apparent, with stronger differences in the Cladocopium to Durusdinium ratio between high and low pollution sites. By 2019, coral fragments at all sites shifted significantly toward Durusdinium, with sites considered moderately or highly polluted being at or nearly 100% Durusdinium. This shift over time may have occurred due to symbiont shifts after bleaching, such as after any of the four mass bleaching events that occurred between 2014 and 2019 (Witze, 2015; Sudek and Lawrence, 2017; Morikawa and Palumbi, 2019; Figure 5). However, it should be noted that colonies sampled in 2019 were not the same colonies as those sampled in 2014, so the significant shift toward Durusdinium we describe should not be interpreted as single colony "shuffling, " but rather a landscape-scale phenomenon. Since Durusdinium outcompete other symbiont species in environmentally stressed corals, they are predicted to continue to overtake coral symbiont communities over time, especially after continual bleaching events (Stat and Gates, 2011;Logan et al., 2021). Shifts to Durusdinium typically increase coral thermotolerance, but are accompanied by tradeoffs, including to growth rate (Little et al., 2004). While growth rate was not included in this study, it is possible that corals with high levels of Durusdinium (e.g., corals in high pollution sites and/or most corals in 2019) also had lower growth rates or other physiological differences not measured here. In 2019, we observed an interesting exception to the established idea that corals hosting more Durusdinium are more thermotolerant. Though the majority of 2019 colonies hosted 100% Durusdinium, the two most thermotolerant coral colonies both hosted a small proportion of Cladocopium. The qPCR method used in this study to distinguish between Symbiodiniciae genera cannot distinguish among different species of Cladocopium. Therefore, it is possible that we have captured multiple different Cladocopium species, some of which have been shown to promote higher thermotolerance (Hoadley et al., 2021). Also interestingly, even corals that hosted 100% Durusdinium varied in their thermotolerance, indicating that factors other than symbiont community can play a significant role in setting thermotolerance. Previous work has also shown that corals with similar symbiont communities may exhibit variation in thermotolerance due to differences in their environment (Oliver and Palumbi, 2011).
The shift in symbiont communities toward Durusdinium from 2014 to 2019 is likely due to severe marine heatwaves that occurred in the interim (Supplementary Figure 1). As symbiont communities become dominated by heat-tolerant Durusdinium, this facet of coral's adaptive capacity to future warming may become effectively exhausted. Global projections suggest that anthropogenic warming will favor shifts to heattolerant symbiont types under future climate change scenarios, and associated thermotolerance gains may reach an upper limit (Logan et al., 2021). This suggests that other types of acclimatization, such as changes in host gene and protein expression, may play an especially important role in determining thermotolerance when symbiont-driven thermotolerance is maximized for species that can undergo symbiont shifts.

Thermotolerance Did Not Vary by Pollution in 2019, and Symbiont Community Variation Was Minimal
In 2014, we found differences among sites in their degreeday mortalities as an indicator of thermotolerance; yet by FIGURE 5 | Hypothesized mechanisms underlying thermotolerance in A. hyacinthus in Tutuila, American Samoa over time with green representing Cladocopium and blue representing Durusdinium. In 2014, thermotolerance varied by site, with that variation driven by symbiont community differences. Likely due to the four mass bleaching events that took place between 2014 and 2019, symbiont communities became dominated by Durusdinium, thereby evening out differences in thermotolerance among sites. From 2019, we expect that mechanisms other than symbiont community variation are underlying variation in thermotolerance.
2019, pollution was no longer significant in our linear mixed model (p = 0.87, Table 2). Thermotolerance differences in 2014 were related to variation in symbiont community among sites. By 2019, these differences were reduced and corals at all sites hosted almost entirely Durusdinium (with the exception of one individual). Since symbiont community was a strong predictor of thermotolerance in 2014, we attribute the minimal variation in thermotolerance in 2019 to the minimal variation in symbiont community. In a linear mixed model, we found that log2 Cladocopium to Durusdinium ratio was no longer predictive of thermotolerance in 2019 (p = 0.24, Table 2). The minimal predictive power of symbiont community in 2019 may be attributed to shifts toward Durusdinium. Previous work has shown that symbiont fidelity may promote higher thermotolerance compared to symbiont flexibility (i.e., symbiont shifts), especially as co-evolution occurs between symbionts and coral hosts (Howells et al., 2020). Our findings suggest that Tutuila corals are shifting toward Durusdinium dominance, and may continue to shift toward even higher proportions of Durusdinium after future bleaching events.
Nevertheless, it is possible that differences in thermotolerance variation between 2014 and 2019 may be attributed to the differences in the heat stress assays. In 2014, a threshold-based assay (3 days) was used to determine differences in mortality while in 2019, a ramp-hold assay (6 h) was performed to determine differences in bleaching intensity. However, prior comparisons of the acute ramp-hold assay to longer-term assays (i.e., 2-3 weeks) showed similar response differences Evensen et al., 2021). Regardless, for this reason, we did not directly compare differences in thermotolerance over time, but focused on how differences in thermotolerance varied by other factors, including site. Yet, it is possible that variation among sites may still be affected by the difference in the assay between years. For example, the 2014 threshold-based assay tests a more sustained acute heat stress (3 days) whereby corals may exhibit greater variation in their ability to acclimate (perhaps due in part to pollution level) while the 2019 ramp-hold assay tests shorter acute heat stress (6 h) that may not allow for as much variation in acclimation by pollution level.

Gene Expression Correlates With Variation in Thermotolerance in 2019
While thermotolerance differences among sites were minimal in 2019 even after controlling for corals that host predominantly Durusdinium, there was variation in thermotolerance within sites. Since pollution level and symbiont community did not explain these within-site thermotolerance differences, we investigated baseline gene expression levels and gene expression responses to heat stress. Gene expression at control conditions was linked to thermotolerance during heat stress in six gene modules (Figure 4). Interestingly, those gene expression patterns did not appear to be dictated by sample site, meaning that regardless of site, some corals express genes that correlate with future heat stress tolerance. In a linear mixed-effects model of thermotolerance, three of those six modules explained 54.1% of the variation in thermotolerance (Table 5).
At control conditions, gene expression modules were correlated to thermotolerance as measured in the heat stress assay to examine the possibility of "front-loading" (sensu Barshis et al., 2013). We compared the genes in the four modules where baseline expression correlated to high thermotolerance (e.g., modules C5, C6, C7 and C8; Figure 4) with frontloaded genes identified in Barshis et al. (2013). The four modules correlating to high thermotolerance in A. hyacinthus matched three of the 135 genes identified in Barshis et al. (2013), on the nearby American Samoa island of Ofu. These three genes were annotated as: a large repetitive protein, a non-collagenous (NC) domain protein, and a protein kinase family protein (Barshis et al., 2013). We examined the Biological Process (BP) Gene Ontology (GO) terms for these three shared genes, though no BP annotations were available. However, one gene was annotated for oxidoreductase activity as a Molecular Function (MF) GO term. We note that the analysis of gene expression differed between the two studies, which may have led to fewer similarities in frontloading genes than expected. We also compared these four modules (e.g., modules C5, C6, C7, and C8; Figure 4) to three significant bleaching-associated modules from Rose et al. (2016). We found 131 genes shared with the 2459 genes in module 1, 38 genes shared with the 442 genes in module 10, and 10 genes shared with the 277 genes in module 12 (Rose et al., 2016). We examined the Biological Process GO terms in these 179 shared genes with bleaching-related genes from Rose et al. (2016) and found five occurrences of "apoptosis, " "immune response" and "purine nucleotide biosynthetic process, " four occurrences of "homophilic cell adhesion" and three occurrences of "defense response to virus, " "DNA replication, " "innate immune response, " "negative regulation of viral reproduction, " "proteolysis, " and "regulation of transcription, DNA-dependent." The higher number of shared genes in the Rose et al. (2016) comparison vs. the Barshis et al. (2013) comparison may be due to the greater number of genes to compare against or due to the more comparable WGCNA study design. Additionally, both the Barshis et al. (2013) and the Rose et al. (2016) studies examined A. hyacinthus from Ofu, an island in American Samoa hosting remarkably thermotolerant corals. Regardless, these findings suggest that gene modules predictive of thermotolerance may vary regionally or over time, and that functional categories may be more important as predictors than individual genes.
In the least thermotolerant corals, baseline gene expression included upregulation of apoptosis and ion transport Gene Ontology (GO) categories, which are characteristic of the cellular stress response (CSR; Kültz, 2005). These patterns indicate that corals with baseline gene expression patterns characteristic of the CSR unexpectedly perform worse during heat stress. Expression of apoptotic or programmed cell death related genes could be indicative of severe or chronic stress. This idea has been proposed previously: constitutive expression of stress response genes may not benefit organisms if (1) overexpression of these genes is costly or (2) these genes drive tradeoffs in the stress response (Rivera et al., 2021).
Taken together, our results indicate that cytokine production and upregulation of immune response genes at baseline conditions may benefit A. hyacinthus during heat stress while upregulation of apoptosis-related genes may hinder thermotolerance. The most thermotolerant corals may use other gene pathways to protect against macromolecular damage without triggering apoptosis (Rivera et al., 2021). These results align with a finding that disease-tolerant corals upregulated cytokine-related pathways under stress while disease-susceptible corals upregulated apoptotic-related pathways (Fuess et al., 2017). These differences in baseline gene expression may be due to variables that were not measured in this study, including environmental, ecological, or evolutionary variation. It is also possible that the differences in baseline gene expression are not plastic, but rather signal differences in heritable variation in thermotolerance. Since genetic differences were not explored in this study, we cannot determine whether these gene expression differences occur due to plastic processes or genetic variation, though past work on A. hyacinthus in American Samoa has found that both adaptation and acclimatization shape thermotolerance (Oliver and Palumbi, 2011;Palumbi et al., 2014;Barshis et al., 2018;Thomas et al., 2018).
The linear mixed model that best predicted thermotolerance in 2019 contained eigengene expression of three gene modules ( Table 5). These modules encompassed 681 genes whose baseline expression could be used to predict 54.1% of variation in future thermotolerance. This can be compared to the 46.9% of thermotolerance variation in 2014 that could be explained by season, log2 Cladocopium to Durusdinium ratio, and pollution level. However, if we had baseline gene expression data available for 2014, it is possible that more variation in thermotolerance could have been explained by the model.
Baseline expression of these three gene modules were shown to predict a substantial amount of variation in future thermotolerance. This presents a useful tool for reef managers or restoration specialists, who may test baseline expression to identify corals with high thermotolerance. However, additional research with higher sample sizes and broader geographic areas would need to be conducted before useful baseline gene expression lists could be generated. Our gene list was useful in predicting thermotolerance of our samples, yet, when compared to published gene lists (e.g., Barshis et al., 2013;Rose et al., 2016), we find few genes overlapping. This suggests that baseline gene expression data may be context dependent and may be most useful when applied to the geographic region where they were generated. Further research into baseline gene expression and its ability to predict thermotolerance may give insight into its potential as a conservation tool.

Baseline Gene Expression Related to Pollution Level
In a WGCNA analysis of baseline gene expression patterns, module C9 showed opposing patterns in the low pollution vs. the high pollution sites (Figure 4). When compared to Gene Ontology (GO) categories upregulated under nutrient stress in Rosic et al. (2014), none of the categories matched Biological Processes GO terms in module C9. We also compared the GO terms of genes in module C9 to the generalized Acropora stress response GO terms from Dixon et al. (2020) and found five of the 190 genes shared GO terms. These five genes were primarily mini-collagen and calcium-binding proteins. While few of the GO terms in module C9 matched previous studies, many of the GO terms appear to be involved in production of methionine, which is an amino acid that has been shown to mitigate oxidative stress (Luo and Levine, 2009;Aguilar et al., 2017). Therefore, it is possible that pollution is inducing stress response genes that differ from those identified in the two previous studies. This may be due to the highly context-dependent nature of pollution, whereby differences in pollution may induce different stress response genes. This suggests that pollution may be inducing some level of cellular stress, even though pollution did not correlate with thermotolerance in 2019.

Gene Expression Under Heat Stress Related to Pollution Level, Symbiont Community, and Thermotolerance
During heat stress, gene expression correlated with pollution but showed more significant correlations with symbiont community (Supplementary Figure 6). Two gene modules showed opposite patterns in corals hosting entirely Durusdinium compared to those hosting Cladocopium and Durusdinium. This follows previous work showing that Symbiodiniaceae type can affect gene expression in the coral host (Yuyama et al., 2012;Barfield et al., 2018;Helmkampf et al., 2019). The modules that varied by Symbiodiniaceae type contained overrepresented Biological Processes (BP) Gene Ontology (GO) categories relating to RNA splicing, processing, and gene silencing. Barfield et al. (2018) also found that RNA processing and modification genes were upregulated in corals hosting Cladocopium compared to those hosting Durusdinium. This suggests that maintaining symbiosis with Cladocopium may require post-transcriptional modifications (Baumgarten et al., 2017;Barfield et al., 2018).
Gene expression under heat stress was also correlated to thermotolerance, with four modules related to the 10-20% most thermotolerant corals and two modules related to the 10-20% least thermotolerant corals (Supplementary Figure 6). Module H15, whose expression was upregulated in the least thermotolerant corals, was significantly enriched for MAPK signaling. Mitogen-activated protein kinases (MAPK) are signaling proteins involved in repairing oxidative damage that occurs during stress (Kültz, 2005). The least thermotolerant corals may express higher levels of stress response genes during heat stress because they are encountering greater macromolecular damage than more thermotolerant corals.
Interestingly, two of the modules that were upregulated in Cladocopium-containing corals under heat stress were also upregulated in the most thermotolerant corals. This is unexpected since previous work shows that colonies hosting entirely Durusdinium are typically more thermotolerant (Berkelmans and van Oppen, 2006;Stat and Gates, 2011;Howells et al., 2020). However, this pattern is driven by two individuals mentioned above hosting low levels (<1%) of Cladocopium which were among the top performing 10%.

CONCLUSION
This study explored the impact of pollution on coral thermotolerance, symbiont communities and gene expression in a field-based experiment. While our field-based experiment cannot disentangle the specific impacts of pollution on these variables, it is useful to study these effects in situ to determine how real-world environments vary. Symbiont communities showed trends with pollution level, with more polluted sites hosting higher proportions of heat-tolerant Durusdinium. Yet, by 2019, all sites overwhelmingly hosted Durusdinium, demonstrating a noticeable shift in symbiont communities from 2014, which contained much higher levels of heat-sensitive Cladocopium. Thermotolerance was not determined by symbiont communities nor pollution level by 2019, but did relate to baseline gene expression patterns prior to heat stress. This suggests that differences in baseline gene expression may allow some corals to better tolerate subsequent heat stress. We found that baseline expression of three gene modules was predictive of future thermotolerance. Future work should investigate what triggers these differences in baseline gene expression to better understand how management efforts can use them in conservation efforts. This study highlights how gene expression patterns could aid in the identification of heat-tolerant corals in a future where most corals are dominated by Durusdinium and symbiont-driven thermotolerance has reached an upper limit. Further, the work serves as a cautionary tale that we need to routinely re-examine our assumptions about patterns of coral thermotolerance as climate change effects become increasingly manifest.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: the raw RNA sequencing reads are available on the NCBI Sequence Read Archive (SRA) database (Bioproject accession: PRJNA762371; SRA accession: SRP339664). All data and code are available at https://github.com/melissanaugle/SCLERA_Tutuila_ Thermotolerance.

AUTHOR CONTRIBUTIONS
TO, DB, and CL collected the data in 2014. MN collected the data in 2019 and wrote the first draft of the manuscript. MN and TO performed the statistical analysis. MN, TO, and CL contributed to the interpretation of the results and wrote sections of the manuscript. All authors, excepting RG who sadly could not, contributed to manuscript revision, read, and approved the submitted version. All authors contributed to the conception and design of the study.

ACKNOWLEDGMENTS
We would like to thank Griffin Srednick, Hideyo Hattori, Johann Vollrath, Valentine Vaeoso, and Salvador Jorgensen for their roles in the 2014 field collections. We thank Jennifer Grossman and Casey Juliussen, who were supported by the Undergraduate Research Opportunities Center at California State University Monterey Bay (CSUMB), for their roles in 2019 field collections. We also thank J. Steve Ryan for his CBASS assistance. We additionally thank Melanie (Cuijuan) Yu, Jacoby Baker, Juliana Cornett, Arie Dash and the 2020 and 2021 Marine Experimental Physiology courses at CSUMB for their contributions to RNAseq sample preparation and preliminary data analysis. We would also like to thank two reviewers for their insightful comments and efforts toward improving our manuscript.