Comparative transcriptomics of spotted seatrout (Cynoscion nebulosus) populations to cold and heat stress

Abstract Resilience to climate change depends on a species' adaptive potential and phenotypic plasticity. The latter can enhance survival of individual organisms during short periods of extreme environmental perturbations, allowing genetic adaptation to take place over generations. Along the U.S. East Coast, estuarine‐dependent spotted seatrout (Cynoscion nebulosus) populations span a steep temperature gradient that provides an ideal opportunity to explore the molecular basis of phenotypic plasticity. Genetically distinct spotted seatrout sampled from a northern and a southern population were exposed to acute cold and heat stress (5 biological replicates in each treatment and control group), and their transcriptomic responses were compared using RNA‐sequencing (RNA‐seq). The southern population showed a larger transcriptomic response to acute cold stress, whereas the northern population showed a larger transcriptomic response to acute heat stress compared with their respective population controls. Shared transcripts showing significant differences in expression levels were predominantly enriched in pathways that included metabolism, transcriptional regulation, and immune response. In response to heat stress, only the northern population significantly upregulated genes in the apoptosis pathway, which could suggest greater vulnerability to future heat waves in this population as compared to the southern population. Genes showing population‐specific patterns of expression, including hpt, acot, hspa5, and hsc71, are candidates for future studies aiming to monitor intraspecific differences in temperature stress responses in spotted seatrout. Our findings contribute to the current understanding of phenotypic plasticity and provide a basis for predicting the response of a eurythermal fish species to future extreme temperatures.

(RNA-seq). The southern population showed a larger transcriptomic response to acute cold stress, whereas the northern population showed a larger transcriptomic response to acute heat stress compared with their respective population controls. Shared transcripts showing significant differences in expression levels were predominantly enriched in pathways that included metabolism, transcriptional regulation, and immune response. In response to heat stress, only the northern population significantly upregulated genes in the apoptosis pathway, which could suggest greater vulnerability to future heat waves in this population as compared to the southern population. Genes showing population-specific patterns of expression, including hpt, acot, hspa5, and hsc71, are candidates for future studies aiming to monitor intraspecific differences in temperature stress responses in spotted seatrout. Our findings contribute to the current understanding of phenotypic plasticity and provide a basis for predicting the response of a eurythermal fish species to future extreme temperatures.

K E Y W O R D S
climate change, Cynoscion nebulosus, phenotypic plasticity, RNA-seq, temperature stress, transcriptome

| INTRODUC TI ON
Temperature has direct and pervasive effects on fish physiology (Angilletta et al., 2010;Fry, 1947). Since most fish are ectothermic and their body temperature tracks ambient water temperature, temperature governs the rate of biochemical reactions (Allen et al., 2002), metabolic rates (Chabot et al., 2016), and ultimately the distribution of species (Pinsky et al., 2013). There is ample evidence that fish populations from contrasting thermal regimes show divergent physiological responses to water temperatures, likely due to local adaptation. For example, sockeye salmon (Oncorhynchus nerka) populations were found to differ in their cardiovascular physiology, which correlated with the historical river temperatures each population encountered during upriver migration (Eliason et al., 2011). Common killifish (Fundulus heteroclitus) collected from their northern and southern range limit along the western Atlantic coast show different thermal tolerance; the critical thermal maximum was significantly higher in the southern population (~1.5°C; Fangue, 2006). The underlying molecular mechanisms, however, are complex and only starting to be explored for fishes (Oomen & Hutchings, 2017). This line of research has historically been difficult as quantifying gene expression in nonmodel species was limited to a few candidate genes or required significant investment of time and expertise to generate species-specific resources (e.g., microarrays).
Next-generation sequencing and RNA-sequencing (RNA-seq) have become more accessible and allowed the transcriptome (whole set of the messenger RNA molecules in a cell or tissue) of any organisms to be studied (Alvarez et al., 2015;Todd et al., 2016). Previous studies looking at gene expression differences have generally focused on sublethal temperature stress lasting from weeks to months (Narum & Campbell, 2015;Newton et al., 2013;Veilleux et al., 2018).
How gene expression responds to acute temperature stress on the scale of hours to days is less well understood (Buckley, 2006;Healy et al., 2010). It is predicted that extreme episodes of water temperature change, such as marine heat waves, will increase in frequency due to climate change (Trenberth & Fasullo, 2012). Transcriptomic studies can improve our understanding of the molecular mechanisms underlying these short-term events and help predict which population might be more vulnerable to these events.
Spotted seatrout, Cynoscion nebulosus, is a teleostean fish distributed from coastal waters in New York to the Gulf of Mexico (Bortone, 2002). This species is uncommon north of Chesapeake Bay, most likely due to its low survival in water temperatures below 5°C . Spotted seatrout is comprised of several genetically distinct populations throughout its geographic range (Anderson & Karel, 2009;Seyoum et al., 2018;Weinstein & Yerger, 1976), which provides an ideal opportunity to study differences in thermal plasticity among populations. Along the U.S.
East Coast, spotted seatrout in Chesapeake Bay are genetically distinct from those in South Carolina and further south (McDowell et al., 2015;Wiley & Chapman, 2003). Based on genome-wide single nucleotide polymorphism markers (Song, 2020), the level of genetic differentiation between spotted seatrout sampled from Corrotoman River, Virginia, and those sampled from James Island,  (Smith et al., 2008), size at maturity (Brown-Peterson, 2003;Ihde, 2000), and metabolic rate (Song et al., 2019), but whether these populations also respond differently to thermal stress is unknown.
Fish populations that experience contrasting temperature regimes show distinct transcriptomic responses, likely as a result of local adaptation (Newton et al., 2013;Veilleux et al., 2018). Spotted seatrout living at its northern range limit in Chesapeake Bay encounter near-freezing water temperature yearly (McGrath & Hilton, 2017), while its southern counterparts rarely experience temperatures below 5°C. In summer, however, maximum water temperature in Chesapeake Bay can reach above 30°C, which is similar to temperatures experienced by fish in more southern regions (Lewisetta, VA, and Charleston, SC, National Centers for Environmental Information, National Oceanic and Atmospheric Administration). Episodic winterkills of spotted seatrout are common in Virginia, but rare in South Carolina and further south on the US East Coast .
The differences in winter severity suggest that the two populations are under different selective pressure by temperature and this may have led to distinct transcriptomic responses to acute cold and heat stresses.
Previous study has shown that the northern spotted seatrout population has greater capacity for metabolic plasticity than the southern population from 5°C to 30°C (Song et al., 2019). The purpose of this study was to compare the underlying transcriptomic response to acute temperature stress in these two genetically distinct and physiologically divergent spotted seatrout populations to better understand the observed metabolic differences. The objectives were threefold: (a) to construct a high-quality transcriptome for spotted seatrout, (b) to discover and quantify shared transcriptomic responses to cold and heat stress in both populations, and (c) to discover and quantify unique transcriptomic responses to cold and heat stress in each population.
The F1 generation from wild SC parents was also included to supplement the sample size of the SC wild-caught fish due to logistical difficulties with sample collection and transport (Table S1).
Fish were acclimated in flow-through 10,000-L circular tanks with brackish water (salinity 20-22 ppt) from the York River, Virginia.
Water temperature was maintained at 15°C prior to cold stress experiments and at 20°C prior to heat stress experiments. Different acclimation temperatures were used for heat and cold stress experiments to achieve the same level of temperature change (10°C decrease or increase; see next section). Cold-stressed fish from both populations were collected and experimented upon during the fall/winter months, whereas the heat-stressed fish were collected during the fall or spring and experimented upon during the summer. The acclimation temperatures were chosen because they represent common water temperatures that spotted seatrout commonly experience in the fall and spring, respectively. Differences in acclimation lengths (30-105 days; Table S1) were due to different sampling dates, as well as experimental dates (Song et al., 2019).
Nevertheless, all fish were acclimated at least 30 days prior to cold and heat stress. All fish were fed frozen and thawed bay anchovies (Anchoa mitchilli) every 2 days to satiation.

| Experimental design
Experimental fish were held in cylindrical respirometers immersed in water as part of a previous study measuring metabolic rates of spotted seatrout at various temperatures (Song et al., 2019). Briefly, fish from both VA and SC in the cold stress group were consecutively exposed to decreasing temperatures starting from the acclimation temperature at 15°C (20 hr), 10°C (5 hr), and 5°C (18 hr) at a rate of −2.5°C/hr. The rate of temperature decrease was chosen because winterkills of spotted seatrout typically involve rapid temperature decreases due to cold fronts or snow melt . Fish in the heat stress group were exposed to increasing water temperatures: 20°C (20 hr), 25°C (5 hr), and 30°C (18 hr) at a rate of 2.5°C/hr. This rate of temperature increase was chosen to mirror the cold shock experiment and to simulate potential heat shocks in the shallow estuarine habitats for spotted seatrout. At the end of each respirometry experiment, fish were euthanized by applying cranial concussion followed by pithing. This 2-step protocol follows the American Veterinary Medical Association Guideline for the Euthanasia of Animals (Leary & Johnson, 2020). This protocol is faster than other euthanasia methods, such as immersing fish in icy water mix or using tricaine methanesulfonate (MS-222), and therefore minimized the impact of euthanasia on gene expression. A small piece of liver tissue from each fish was collected after euthanasia using sterilized scissors. Liver tissue was used in this study because it is a key regulator for metabolic processes and produces many stress-responsive proteins (Currie & Schulte, 2013). The tissues were stored in cryovials and flash-frozen in liquid nitrogen until RNA extraction. Spotted seatrout in the control groups were directly netted out of the holding tanks (15°C), and liver tissues were preserved as

| RNA extraction and RNA-seq
Liver tissue (~20 mg) from each sample was weighed and pulverized in liquid nitrogen using a mortar and pestle. The mortar and pestle were thoroughly cleaned by rinsing under deionized water and wiped with RNase AWAY solution (Thermo Scientific) between samples to ensure no RNA contamination occurred. Total RNA was extracted from the pulverized tissues using the RNeasy Mini Kit (Qiagen) according to the manufacturer's protocol including the on-column DNase digestion step to eliminate genomic DNA. RNA concentration was determined using a Qubit 2.0 fluorometer (BR RNA Assay, Invitrogen). RNA quality was assessed in two ways. First, the ratio of absorbance at 260 and 280 nm (260/280) was evaluated using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies).
Samples with a 260/280 around 2 suggest high purity of RNA.   (Andrews, 2015). To reduce the quantity of the input reads for de novo assembly and to improve assembly efficiency (Brown et al., 2012;Fletcher et al., 2013), forward and reverse reads from each sample were first concatenated and in silico normalization was performed using the default setting in Trinity v2.8.4 (Haas et al., 2013). Different assembly algorithms can complement each other in discovering genes that might be missed using a To assess the quality of the newly generated spotted seatrout transcriptome, QUAST v4.6.3 (Gurevich et al., 2013) was used to calculate common genomic metrics, including transcript length summary, N50 (length of at least half of all the transcripts) and

| Bioinformatic analyses
GC content. BUSCO v3.0.2 (Benchmarking Universal Single-Copy Orthologs) was used to assess the completeness of the assembled transcriptome in terms of expected gene content (Simão et al., 2015).
Specifically, the spotted seatrout liver transcriptome was searched against the database actinopterygii_odb9, which contains 4,584 evolutionarily conserved genes expected to be found as single-copy orthologs in at least 90% of the Actinopterygii (ray-finned fishes).
To assess gene expression levels for each transcript in the transcriptome for all 30 fish, original sequences were mapped back to the new transcriptome using kallisto v0.43.1 (Bray et al., 2016).
Kallisto uses a novel "pseudoalignment" approach to eliminate the need for perfectly aligning individual bases, thereby reducing the computing time by two orders of magnitude compared with alternative programs while achieving similar mapping accuracy (Bray et al., 2016). Transcript abundances were normalized and reported in transcripts per million (TPM).
To assess differential expression, two programs with different statistical frameworks were used: DESeq2 (Love et al., 2014) and limma + voom (Law et al., 2014;Ritchie et al., 2015). Currently, there is no clear consensus on which differential expression algorithm achieves the best balance between type I and type II errors (Costa-Silva et al., 2017;Soneson & Delorenzi, 2013). Similar to our use of different de novo assemblers, we used more than one approach to detect consensus transcripts because a false positive is less likely to be identified twice. DESeq2 uses a negative binomial distribution and a shrinkage estimator for dispersion and fold change.
Limma + voom uses normal distribution and fits a linear model with all genes and samples combined. Both methods have been shown to achieve a good balance between accuracy and sensitivity when the number of biological replicates is at least three ( (Benjamini & Hochberg, 1995). DEGs identified by both methods were retained and used in downstream analyses.
The new transcriptome was annotated using Trinotate v3.1.1 (Bryant et al., 2017), which conducts BLASTx and BLASTp searches against the Swiss-Prot database to identify matches to known proteins using default E-value cutoffs (Altschul et al., 1990 et al., 2006) was used with the classic Fisher's exact test, p < .01. All GO terms found using Trinotate were used as the gene universe ("reference") and DEGs were used as target lists ("interesting genes").

| Validation using RT-qPCR
To validate the accuracy of RNA-seq results, real-time quantitative polymerase chain reaction (RT-qPCR) was used to calculate log 2 fold change of a subset of the transcripts (n = 6). Primers were designed based on transcripts assembled in the de novo transcriptome using Primer-BLAST (Ye et al., 2012), with primer length set to 20 bp and melting temperature set to 60°C (Table S2). 18S ribosomal RNA was used as an internal control using primers from Brewton et al., 2013, and VA and SC control group samples were used as reference samples. and 95°C (15 s). Each reaction was performed in triplicate, and the mean threshold cycle (C T ) was used for subsequent analyses. The comparative C T method was used to present log 2 fold change in order to make results directly comparable to RNA-seq results. To assess the validity of using the comparative C T method, PCR efficiencies of all primers were calculated by running a standard curve with five 1:10 serial dilution points (Schmittgen & Livak, 2008).

| Transcriptome assembly
All 30 samples had an RNA Integration Number (RIN) value greater than 8, indicating high-quality RNA (Figure 1). Sequencing reads per sample ranged from 13.6 million to 18.6 million (mean ± SD = 15.7 ± 1.2), and visual examination of fastQC reports confirmed good average quality scores across all bases (average quality score > 28). After in silico normalization, 6.6 million pairedend reads were retained for de novo transcriptome assembly. Four different assemblers produced the following numbers of transcripts with the results from different k-mers combined: Trinity,185,556;SOAPdenovo,403,848;TransAbyss,209,799;and Velvet,416,310.
The final spotted seatrout liver transcriptome consisted of 37,398 transcripts.
The quality of the de novo-assembled spotted seatrout liver transcriptome was assessed using QUAST and BUSCO. QUAST indicated that this transcriptome contained 21,316 transcripts longer than 500 bp, an N50 of 3,121 bp, and a GC content of 48.98% (Table   S3). BUSCO analysis found that this transcriptome assembly contained 81.3% complete gene sequence information and 6.1% fragmented BUSCOs of the essential genes in bony fish genomes (lineage dataset: actinopterygii_odb9, 20 species, 4,584 total BUSCO groups searched).

| Gene expression plasticity among and within populations
Gene expression levels showed good correlation among biological replicates within treatment and control groups (Pearson's r = .89, Figure 2). While the two populations clearly separated in the control groups, they did not form distinct clusters in heat and cold treatment groups (Figure 2, Figure S1). A total of 1,653 DEGs were found in VA and SC samples combined in response to cold and heat stress treatments. Of these, 1,281 DEGs were found in the cold treatment groups, and 570 DEGs were found in the heat treatment groups as compared to controls. A set of 20 DEGs were responsive to both heat and cold stress in both populations (File S3, Figure S2).

| Functional analysis
There were 53 common pathways with 29 unique pathways in VA and 18 unique pathways in SC. A list of all the molecular pathways can be found in the supplemental information (File S6).
GO enrichment analyses revealed common and unique biological pathways between populations in response to thermal stress (Table 3) Figure S5).

| D ISCUSS I ON
Transcriptomic studies of fish populations originating from contrasting temperature regimes have shown distinct responses to chronic thermal stress, but little is known about the effects of acute thermal stress. We conducted RNA-seq on two genetically distinct and physiologically divergent populations of spotted seatrout exposed to acute temperature stress and de novo-assembled the first liver transcriptome for this species. Based on this transcriptome, we compared the molecular mechanisms underlying the response to temperature stress between the two populations. Log2 fold change of a subset of genes between RT-qPCR and RNAseq were highly consistent, validating the accuracy of the RNA-seq results.
A high-quality transcriptome is essential for differential expression analyses and functional annotation (Grabherr et al., 2011 indicates that this population experienced greater physiological stress than the southern population when temperatures were elevated, likely due to a mismatch between oxygen demand (high SMR) and oxygen supply (low dissolved oxygen) (Pörtner & Knust, 2007).   (Long et al., 2013). In common killifish, cold acclimation at 5°C induced 5,460 upregulated genes compared with 1,746 downregulated genes in muscle transcriptomes (Healy et al., 2017). This trend is predicted because low temperature depresses the rate of biochemical processes, and increased expression can compensate for this kinetic restraint (Currie & Schulte, 2013 Protein processing in endoplasmic reticulum (9) lco00100 Steroid biosynthesis (8) lco04010 MAPK signaling pathway (8) lco04920 Adipocytokine signaling pathway (7) lco00561 Glycerolipid metabolism (7) lco04210 Apoptosis (7) lco04371 Apelin signaling pathway (7) lco03320 PPAR signaling pathway (7) lco04530 Tight junction (6) lco04910 Insulin signaling pathway (6) lco04110 Cell cycle (6) lco04530 Tight junction (6) Heat stress We found 20 shared DEGs in both cold and heat stress in both populations ( Figures S2 and S3). These genes most likely play important roles in the cellular stress response (CSR). Seven out of the 20 genes were annotated, and associations included the heat-shock protein 90-alpha (hs90A), apolipoprotein (apom), and haptoglobin (hpt). Heat-shock proteins are a group of well-studied gene families, which act as generic molecular chaperones to help maintain protein integrity under a range of stressors such as temperature, salinity, and disease (Iwama et al., 2001). Apolipoproteins bind to lipids and play a major role in lipid transport and have been shown to be a part of the innate immunity in fishes (Causey et al., 2018;Concha et al., 2004;Pereiro et al., 2012). Hpt is one of the acute-phase stress proteins synthesized by the liver (Cordero et al., 2017;Windisch et al., 2014) and can bind to hemoglobin in the plasma and reduce oxidative stress (Alayash, 2011). While most of the 20 genes showed similar levels of differential expression between VA and SC ( Figures S2 and S3), hpt showed the largest differences in log2 fold change ( Figure S4, cold stress: 10 versus. 3, heat stress: 10 versus. 5.5). This substantial difference was in part due to the low average expression level of hpt in SC controls. We suggest that hpt warrants further investigation to confirm this finding, and if the results are repeatable, it would be an excellent genetic marker for estimating population-level temperature stress in spotted seatrout. Unfortunately, 13 out of the 20 generic stress genes did not return a BLAST hit, highlighting the problem with poor functional annotation in nonmodel organisms (Pavey et al., 2012).

| Shared KEGG pathways
A few shared molecular pathways between the northern and southern populations accounted for a disproportionately large number of DEGs (Table 2). Metabolic pathways accounted for the most DEGs.
The MAPK signaling pathway was one of the top pathways observed in response to heat stress in spotted seatrout and is one of the hallmarks of CSR in other organisms (Kültz, 2004). MAPKs are kinases that are involved in protein phosphorylation cascades in order to regulate the expression of many downstream targets (Cowan & Storey, 2003;Huang et al., 2017). This pathway has been shown to be particularly important in modulating gene expression in gill epithelium cells in killifish during osmoregulation (Kültz & Avila, 2001) and in the swim bladders of channel catfish (Ictalurus punctatus) in response to low dissolved oxygen (Yang et al., 2018

TA B L E 3 (Continued)
stress, the FoxO signaling pathway was a top pathway in both populations. Similar to the MAPK pathway, the FoxO signaling pathway can effect global transcriptomic change via a group of transcription factors that target genes involved in apoptosis, oxidative stress resistance, and cell cycle control (Puig & Mattila, 2011;Ronnebaum & Patterson, 2010).

| Population-specific genes and KEGG pathways
In response to cold stress, the northern population downregulated acyl-coenzyme A thioesterase (acot) in the fatty acid elongation (lco00062) and biosynthesis of unsaturated fatty acid (lco01040) pathways. Acot plays a critical role in fatty acid metabolism and ATP generation (Tillander et al., 2017). This suggests that fat metabolism is suppressed, at least in the short term, in response to cold stress.
In addition, genes involved in protein synthesis and transport were upregulated in pathways such as ribosome biogenesis in eukaryotes (lco03008) and protein export (lco03060). These genes included ribonuclease P protein subunit (pop4), U4/U6 small nuclear ribonucleoprotein (snu13), and transport protein subunit alpha (sec61a).
26S proteasome regulatory subunit T3 (psmc4) in the proteasome (lco03050) pathway is also upregulated. The proteasome is a large molecular complex where protein degradation occurs (Glickman & Ciechanover, 2002). Protein synthesis, export, and degradation are all energetically expensive processes requiring ATP. The upregulation of genes in these pathways suggests the northern population is capable of producing sufficient ATP to meet energy requirements at low temperature.
Under acute cold stress, the southern population was found to uniquely upregulate long-chain acyl-CoA synthetase (acsl) and long-chain fatty acid-CoA ligase (acsbg) in the fatty acid biosynthesis (lco00061) and fatty acid metabolism (lco01212) pathways.
Both genes play critical roles in the fatty acids oxidation (Cheng et al., 2017). The large subunit ribosomal protein L23Ae (rp-l23ae) in the ribosome (lco03010) pathway is upregulated, suggesting enhanced protein production. However, genes in the protein degradation pathway, lysosome (lco04142), are downregulated: deoxyribonuclease II (dnase2), lysosomal alpha-glucosidase (gaa), and acid ceramidase (asah1). Lysosomes are organelles in eukaryotic cells containing hydrolytic enzymes that degrade biomolecules (Levine & Klionsky, 2004). A mismatch between protein synthesis and protein degradation may suggest the increased accumulation of misfolded proteins in the cell and can lead to apoptosis (Fribley et al., 2009).
The northern and southern populations of spotted seatrout also showed distinct responses to heat stress. Under acute heat stress, the northern population uniquely upregulated heat shock 70 kDa protein (hspa5), a member of the heat-shock protein 70 family (Roberts et al., 2010). Acetyl-CoA synthetase (acss1) was upregulated in the pyruvate metabolism (lco00620) and carbon metabolism (lco01200) pathways. Acss1 catalyzes the reaction that produces the raw material, acetyl-CoA, for the citric acid cycle (aka. tricarboxylic acid cycle) in order to produce ATP. However, this pathway is only activated when cells have depleted their normal carbon source (pyruvate) (Wolfe, 2005). Isocitrate dehydrogenase (idh1), which is involved in the citric acid cycle, was also uniquely upregulated. The higher metabolic demands of the northern fish observed in previous physiology studies (Song et al., 2019) may have quickly exhausted pyruvate, and the upregulation of acss1 serves as a short-term solution to supply acetyl-CoA to generate ATP. In addition, genes in the N-glycan biosynthesis (lco00510) pathway were upregulated: alpha-1,2-mannosyltransferase (alg9) and mannosyl-oligosaccharide glucosidase (mogs). These enzymes are involved in N-linked glycosylation (addition of oligosaccharides to proteins) and play important roles in protein stability and cell signaling (Sinclair & Elliott, 2005). GO analysis also identified glycosylation-related processes to be enriched in VA fish (Table 3. GO:0006487. GO:0009101, GO:0009100).
Similarly, under heat stress, the southern population uniquely upregulated heat-shock cognate 71 kDa protein (hsc71). Hsc71 belongs to the same protein family as hspa5, which was uniquely upregulated in the northern population, and performs similar functions; however, hsc71 is known to be constitutively expressed in unstressed cells and hspa5 is only upregulated during stress (Goldfarb et al., 2006;Roberts et al., 2010). Hsc71 was also found to be inducible under heat stress in gill tissues of Gillichthys mirabili (Buckley, 2006). RAC serine/threonine-protein kinase (akt) is downregulated in the following pathways: mTOR signaling pathway (lco04150), insulin signaling pathway (lco04910), and apelin signaling pathway (lco04371).
G1/S-specific cyclin-D2 (ccnd2) was downregulated in the hedgehog signaling pathway (lco04340). Both akt and ccnd2 play a role in cell proliferation and the cell cycle (Katoh & Katoh, 2009;Manning & Cantley, 2007); thus, their downregulation may indicate heat stress also induces cell growth arrest in the southern population.
This study, which examined transcriptomic responses to acute temperature stress in spotted seatrout, provided novel insights into stress response. For example, cytoskeletal re-organization is among the significantly enriched biological processes in studies, which applied a chronic temperature stress (Healy et al., 2017;Metzger & Schulte, 2017;Narum & Campbell, 2015;Newton et al., 2013;Sun et al., 2019); however, we did not find this process to play a significant role under acute temperature stress. This finding is consistent with a recent study that applied an acute heat stress protocol (2°C/ hr, 20-32°C) to small yellow croaker, Larimichthys polyactis (Chu et al., 2020). This discrepancy suggests that gene expression underlying transcriptomic plasticity may act on different timescales (Metzger & Schulte, 2017) and requires more studies to understand the commonalities of response to acute temperature stress in fishes.
The fish used in this study were part of an ongoing respirometry study, and there are potential caveats that should be considered when interpreting the results. One limitation concerns the origin of the SC control samples, which were not wild-caught, but instead were the F1 generation of wild parents. Although the controls were from the F1 generation, the rearing tanks received water directly from Charleston Harbor, and thus experienced the same natural temperature fluctuations as their wild counterparts. Admittedly, other aspects of being in captivity may have altered gene expression, but all fish in the treatment group experienced acclimation in this same environment. Therefore, we contend that the fish in the control group, although not optimal, were appropriate for this study.
Secondly, gene expression patterns are known to be affected by seasonality (Tian et al., 2019) and sex (Selmoni et al., 2020). Due to logistic difficulties, all fish were captured in fall/winter (October/ November) with the exception of the SC fish used in the heat stress experiments, which were captured in spring (March) ( Table S1).
Despite the evidence that eurythermal fishes require much less time (days to a week) to acclimate than more stenothermal fishes (weeks to months) (Healy & Schulte, 2012), it is not well understood how differences in acclimation times affect gene expression in eurythermal fishes. Finally, this study used a mix of sexes based on visual examination of gonads. Thus, it is possible that there were sexspecific transcriptomic responses to temperature stress that were unaccounted for in this study.
In summary, we present evidence that differences in thermal tolerance, as reflected by differential gene expression, contribute to the population-level divergence observed in spotted seatrout from the U.S. South Atlantic and provide mechanistic insights into physiological responses to acute temperature stress. Our results indicate that the northern population shows transcriptional signatures consistent with cold adaptation, yet they may be more vulnerable to acute heat waves than the southern population. Likewise, the transcriptomic profiles of the southern fish indicate that they have a greater stress response than northern fish when exposed to cold temperatures. These differences may play a role in maintaining the observed genetic differences between these two populations. The liver transcriptome represents a valuable resource for future genetic monitoring studies. Candidate genes (e.g., hpt, acot, hsc71, hspa5) identified in this study should be functionally validated or screened more broadly across the species' range to verify their ecological importance. Furthermore, genomic-level and cross-generational investigations would also complement this study by discovering genes under selection and improve our understanding of adaptive evolution in general.

ACK N OWLED G M ENT
We would like to thank Steve Arnott, Tim O'Donnell, and John We are grateful to the three anonymous reviewers for their comments and suggestions, which significantly improved this manuscript.

CO N FLI C T O F I NTE R E S T
None declared.