Ocean Warming Leads to Increases in Aerobic Demand and Changes to Gene Expression in the Pinfish (Lagodon rhomboides)

Anthropogenic climate change is causing increases in the frequency, intensity, and duration of marine heatwaves (MHWs). These short-term warming events can last for days to weeks and can produce severe disruptions in marine ecosystems, as many aquatic species are poikilotherms that depend on the conditions of the environment for physiological processes. It is crucial to investigate the effects of these thermal fluctuations on species that play a disproportionate ecological role in marine ecosystems, such as the pinfish (Lagodon rhomboides) in the Gulf of Mexico and western Atlantic. In this study, we exposed pinfish to a simulated MHW in aquaria and examined the impacts of acute warming on two life stages (juvenile and adult), measuring oxygen consumption and gene expression in two relevant tissue types (liver and muscle). We saw significant increases in routine metabolic rate with increasing temperature in both juveniles (24.58 mgO2/kg/h increase per 1°C of warming) and adults (10.01 mgO2/kg/h increase per 1°C of warming). These results indicate that exposure to increased temperatures was more metabolically costly for juveniles than for adults, on a mass-specific basis. This was also observed in the molecular analyses, where the largest number of differentially expressed genes were observed in the juvenile pinfish. The analyses of gene expression suggest warming produces changes to immune function, cell proliferation, muscle contraction, nervous system function, and oxygen transport. These results indicate that this ecologically relevant species will be significantly impacted by projected increases in frequency and magnitude of MHWs, particularly in the juvenile stage.


INTRODUCTION
Ocean warming is one of the principal impacts of anthropogenic climate change (IPCC, 2014;IPCC, 2021), representing a major threat to the persistence of marine ecosystems. From 1971 to 2010, ocean temperatures rose by an average of 0.11 • C per decade, and warming is projected to continue to increase in upcoming decades even under stringent emission mitigation scenarios (IPCC, 2014). In addition to the long-term increase in average ocean temperature, there have been increases in the frequency, intensity, and duration of acute thermal fluctuations, termed marine heatwaves (MHWs; Frölicher et al., 2018;Oliver et al., 2018;IPCC, 2021). These events are defined as a period of five or more consecutive days in which the local temperature exceeds the 90th percentile based on a 30-year historical period (Hobday et al., 2016). MHWs have approximately doubled in frequency since the 1980s (IPCC, 2021) and will likely continue to become more frequent over the next century (Frölicher et al., 2018;Oliver et al., 2018). The biological impacts of these thermal fluctuations are often extreme, as vital processes of marine poikilotherms, such as swimming, growth, development, and reproduction, are strongly dependent on ambient water temperatures Green and Fisher, 2004;Munday et al., 2008a,b;Hoegh-Guldberg and Bruno, 2010). Given the considerable effects that MHWs can have on marine species, it is crucial to understand and predict the specific impacts that extreme temperature fluctuations may have on keystone species in a given ecosystem.
At the physiological level, elevated temperatures are metabolically costly for ectotherms. To compensate for heightened energetic demands, organisms must increase their resting metabolic rate as temperatures increase (Ege and Krogh, 1914;Fry and Hart, 1948;Clarke and Johnston, 1999;Chabot et al., 2016). Resting (or standard) metabolic rate is an estimate of an organism's minimum energetic requirements under specific environmental conditions (Chabot et al., 2016). As resting metabolic rate increases with temperature, it becomes challenging for an organism to maintain homeostasis (Fry and Hart, 1948). Eventually, past a threshold temperature, aerobic scope (i.e., the capacity for an organism to increase its metabolic rate above maintenance levels, calculated as maximum metabolic rate -resting metabolic rate) declines to zero as oxygen demands surpass the maximum capacity for oxygen uptake, due to physical limitations of the cardiovascular system (Pörtner, 2001;Pörtner and Knust, 2007;Pörtner et al., 2017; but see Jutfelt et al., 2018). Such changes in metabolic demand and oxygen availability at high temperatures can have considerable effects at the cellular and molecular level, leading to activation of genes associated with metabolism, oxidoreductase activity, immunity, and cellular stress response (Veilleux et al., 2015;Qian and Xue, 2016;Bernal et al., 2020;Huang et al., 2020). These detrimental molecular and physiological consequences of warming can ultimately lead to loss of biodiversity and modifications to local fish communities (Wernberg et al., 2013;Cavole et al., 2016;Hughes et al., 2017). Thus, it is fundamental to understand how changes in oxygen consumption and gene expression under thermal stress will affect species that play a disproportionate ecological role in marine ecosystems.
Notably, coastal ecosystems of the Gulf of Mexico and western Atlantic Ocean support high levels of biodiversity but are increasingly threatened by climate change and MHWs (Costello et al., 2010). For example, the extreme oceanic warming event of 2012 affected the northwest Atlantic Ocean, leading to northward shifts in species distributions and alterations to local marine systems along the east coast of the United States (Mills et al., 2013). One key ecological player in coastal systems of this region is the pinfish (Lagodon rhomboides). This broadly distributed marine fish is found in nearshore areas from the Yucatan Peninsula in Mexico to Cape Cod, Massachusetts (Darcy, 1985). During winter and spring, juvenile pinfish inhabit estuaries and seagrass beds close to the shore, migrating offshore to reproduce once they develop into adults (Nelson, 1998). This species is an important link in marine food webs, as it is consumed by other fish species of commercial importance, as well as marine mammals and birds (Darcy, 1985). Furthermore, this migration from estuaries to the ocean makes pinfish a critical component of nutrient transfer among coastal and oceanic systems, as pinfish are widely consumed by predatory fishes as they move offshore (Nelson et al., 2012). Studies have estimated that pinfish alone can comprise nearly a quarter of the total nitrogen available to higher trophic levels in offshore food webs in the Gulf of Mexico (Nelson et al., 2013). Despite the clear ecological importance of this species, there is a gap in our understanding of how pinfish will respond to warming waters. Previous studies indicate that elevated temperatures result in heightened metabolic rate of pinfish (Wohlschlag et al., 1968;Cameron, 1969;Wohlschlag and Cech, 1970;Hoss, 1974). However, the bulk of this work was conducted decades ago, prior to advancements in oxygen monitoring technologies (Clark et al., 2013). Additionally, none of the aforementioned studies have made comparisons between juveniles and adults and the respective changes to their metabolism under heat stress. Since juvenile and adult pinfish are often found in distinct habitats (Darcy, 1985), and different life stages can often have varying thermal sensitivity, even within the same species (Dahlke et al., 2020), it is necessary to compare how energetic requirements and molecular processes may be differentially affected by warming.
In the present study, we aim to understand how pinfish from the northern Gulf of Mexico will be affected by MHWs, by coupling analyses of routine metabolic rate (RMR) with gene expression of the liver and muscle. Specifically, we ask whether physiological and molecular responses to heat stress are conserved across life stages (juvenile and adult) and tissue types (muscle and liver). We predict that exposure to acute heat stress will lead to significant increases in oxygen consumption in both juvenile and adult pinfish, but that thermal stress will be more metabolically costly for juveniles. We also predict that there will be different molecular mechanisms activated to compensate for the effects of warming among different life stages and tissue types, but core genes associated with aerobic metabolism and cellular stress response should be consistently upregulated. This study also presents a high-quality reference transcriptome for L. rhomboides for the first time. Given the relevance of L. rhomboides to estuarine and coastal ecosystems of the Gulf of Mexico and western Atlantic, this study advances our understanding of how marine communities may be affected by changing climatic conditions in coming decades.

Sample Collection
We collected juvenile pinfish (average mass: 2.65 g, ±0.72 g, SD) via seine nets at depths of approximately 1 m from Perdido Pass (30.285 N,87.544 W), Alabama, United States on March 12, 2021. Collections were made in mid-March as this is when juvenile pinfish are highly abundant in coastal seagrass beds. We transported individuals back to aquarium facilities at Auburn University (Auburn, AL, United States), and upon arrival, fish were housed in compartments of 120-gallon tanks (10 fish per compartment, 2 compartments per tank, 2 tanks per recirculating aquarium system, and a total of 6 independent recirculating aquarium systems). Fish were acclimated to the systems for at least 21 days prior to the start of the experiment. Separately, we collected adult pinfish (average mass: 48.9 g, ±20.7 g, SD) on hook and line at depths of approximately 2-5 m from Dauphin Island (30.251 N,88.079 W) and Perdido Pass (30.285 N,87.544 W), Alabama, United States on November 13, 2020. We chose to collect adult pinfish in November because they are abundant in nearshore areas during that time of year. Individuals were transported back to aquarium facilities at Auburn University (Auburn, AL, United States). Adult pinfish were housed separately in subdivided compartments of 120gallon tanks (one fish per compartment, maximum of four fish per tank, two tanks per recirculating aquarium system, and a total of six independent recirculating aquarium systems). Fish were acclimated to the aquarium system for at least 10 days prior to the start of the experiment. The main difference in experimental setup between the juvenile and adult pinfish is that adults were physically separated into individual compartments by barriers in the aquaria, whereas juveniles were maintained in groups of 10, as adult pinfish displayed aggressive behavior when housed together.
All collections were done with the permission of the Alabama Department of Natural Resources (Permit 2021-3-01), and all experiments were done following the guidelines and permits of the Institutional Animal Care and Use Committee (IACUC) of Auburn University (Permit 2020-3708).

Experimental Design
The experimental approaches were similar for both juvenile and adult pinfish but were carried out separately. All individuals were acclimated to the aquarium systems and held at a constant temperature of 21 • C, which is slightly above the seasonal average temperature in Perdido Pass, Alabama from 2011 to 2019 (approximately 19 • C; National Data Buoy Center [NDBC], 1971). The acclimation temperature was set at 21 • C (rather than 19 • C) to ensure that individuals were not placed under additional stress during the acclimation period, because temperature conditions in Perdido Pass at the time of collection were above the seasonal average (adults were collected at 24.1 • C and juveniles at 20.6 • C). Following acclimation to the aquarium systems, we assigned individuals to one of three temperature treatments (Figure 1). The control treatment was maintained at 21 • C, which is near the temperature at the time of collection for both juveniles and adults. We had two acute warming treatments: +3 • C (24 • C) and +6 • C (27 • C), representing two simulated MHWs of varying intensity. The experimental design encompasses the 90th percentile of recorded temperatures at Perdido Pass, Alabama in March and November of 2011-2020 which were 22.1 and 23.1 • C, respectively (National Data Buoy Center [NDBC], 1971). Additionally, the maximum recorded sea surface temperature at Perdido Pass for March is 26.8 • C, and for November is 26.6 • C (National Data Buoy Center [NDBC], 1971). By choosing 24 and 27 • C as the temperatures for our experimental groups, we can accurately and realistically simulate the conditions of a MHW in the wild, making the results from this assay relevant for natural populations. Additionally, using the same temperature scheme for both the juveniles and adults allows us to directly compare the results of pinfish at different life stages. Temperature data from Perdido Pass are not available from years prior to 2011, and as such, a longer seasonal average temperature could not be calculated. For each experiment, temperature treatments were replicated across two independent recirculating aquarium systems. For fish in the +3 and +6 • C treatments, we increased the temperature gradually, at a rate of +1.5 • C/day (similar to natural variations in temperature that might be experienced in the wild; National Data Buoy Center [NDBC], 1971), to minimize thermal shock (Figure 1). For juveniles a total of n = 22 fish were in the control treatment, n = 23 in the +3 • C treatment, and n = 26 in the +6 • C treatment. A similar protocol was followed for the adult pinfish, with n = 9 samples in the control treatment, n = 5 in the +3 • C treatment, and n = 6 in the +6 • C treatment. Sample size varied by treatment due to unexpected mortality during respirometry experiments as a consequence of handling stress.
Fish were maintained at their assigned temperature for 8-11 days (Figure 1). The variation in exposure time was due to the limited number of oxygen consumption chambers, where a maximum of eight fish could be measured over the course of 2 days. Following temperature exposure, 12-16 juvenile pinfish per treatment were taken for measurements of oxygen consumption, and 10 different individuals per treatment had liver and muscle tissue sampled and preserved in DNA/RNAShield (Zymo) for analyses of gene expression. For the adults, approximately half of the individuals in each treatment were used for measurements of oxygen consumption, and all individuals had liver and muscle tissue taken and preserved in DNA/RNAShield for gene expression analyses.

Oxygen Consumption
We measured oxygen consumption using the Loligo fiber optic respirometry system with AutoResp software ver 2.3.0. Prior to respiration measurements, pinfish were fasted for 24 h. Individuals were placed in acrylic chambers (0.6 L for juveniles, 2.5 L for adults), which were then submersed in a trough containing 30 ppt saltwater (300 L for juveniles, 150 L for adults). The temperature of the trough was maintained at the same experimental temperature to which fish had been exposed (i.e., 21 • C for control, 24 • C for +3 • C, and 27 • C for +6 • C), using an AquaLogic Nema 4X temperature controller and heating element. Each acrylic chamber was equipped with two circulating pumps for intermittent respirometry: a flush pump that allowed for influx of freshly oxygenated water to the chamber, and a recirculating pump that allowed water movement within the chamber when the flush pump was turned off. We selected the intermittent respirometry parameters to ensure that the measure period was short enough such that oxygen levels in the chamber FIGURE 1 | Simplified schematic of the experimental design. Wild pinfish were acclimated to the aquaria for 10-21 days. Temperatures for the experimental groups were then raised from the control (21 • C) at a rate of 1.5 • C per day until the appropriate experimental temperature was reached (24 or 27 • C). After 8-11 days of exposure, tissues were taken for analyses of gene expression, and routine metabolic rate was measured.
did not fall lower than 80% of the maximum oxygen saturation, and the flush period was long enough to ensure return to complete oxygen saturation. For juveniles, conditions were set to 800 s flush, followed by a 60 s waiting period (in which the flush pump shuts off and dissolved oxygen content begins to decline non-linearly), and a 500 s measurement period. RMR was estimated for n = 12 fish in the control group, n = 13 fish in the +3 • C group, and n = 16 fish in the +6 • C group. For adults, intermittent respirometry was conducted in a similar fashion, but the cycling parameters were as follows: 800 or 1500 s flush (dependent on whether 800 s was sufficient to ensure complete oxygen saturation of the chamber), 60 s wait, and 300 s measure, and we had sample sizes of n = 4 control fish, n = 4 fish in the +3 • C treatment, and n = 2 fish in the +6 • C treatment.
For both juveniles and adults, fish were allowed to acclimate to the intermittent respirometry conditions overnight, and measurements were taken the following morning for approximately 2.5 h. To account for bacterial oxygen consumption, we carried out 1-h closed respirometry control runs prior to moving fish into the chambers and after the fish were removed (i.e., once measurements were completed). We estimated bacterial respiration over the course of the experiment using a linear regression, and calculated values of bacterial respiration for the duration of the measurement period, which were then used to correct raw measurements of oxygen consumption. We estimated RMR for each fish as the average oxygen consumption rate, or MO 2 (mgO 2 /kg/h) over the measurement period. MO 2 values were estimated by AutoResp software as is the oxygen concentration at time point zero (mgO 2 /L), (O 2 ) t 1 is the oxygen concentration at time point one (mgO 2 /L), V is the respirometer volume (in L) corrected for the volume of the experimental animal, t is the time (in hours) between time point zero and time point one, and B is the weight of the fish (kg). We analyzed the association between temperature and RMR by building general linear models in R ver. 4.0.2 (R Core Team, 2020). Temperature, fish weight, and exposure length were used as independent variables, and RMR was used as the response variable. To assess the possibility of a significant deviation from linearity in the relationship between temperature and RMR, we also built ANOVA models in R ver. 4.0.2 (R Core Team, 2020) using the same independent and response variables as described above, as ANOVAs have no assumption of linearity. The Akaike's Information Criterion adjusted for sample size (AICc) was used to rank models, and the best-fit model was chosen based on the lowest AICc value. We also calculated Q 10 values, which describe the proportional change in respiration rate over a 10 • C temperature increase, based on the average RMR at 21 • C (control group) and 27 • C (+6 • C group) for both juveniles and adults. Q 10 's were estimated using the R package "respirometry, " which uses the following formula: Q 10 = (RMR 2 /RMR 1 ) × [10/(T 2 −T 1 )] (Birk, 2021).

RNA Extraction and Sequencing
We sampled both liver and muscle tissue from juvenile and adult pinfish at the conclusion of the experiment. For the juveniles, we took liver and muscle tissue after 8 days of temperature exposure from n = 10 fish per treatment, and tissue samples were only taken from individuals that had not undergone respirometry measurements. For the adults, however, liver and muscle tissue were taken from the majority of fish (n = 6 in the control treatment, n = 5 in the +3 • C treatment, and n = 6 in the +6 • C treatment) after 8-11 days of thermal stress. Due to the small number of adult pinfish, tissues were taken from both individuals that had undergone respirometry analyses, and those that had not. All tissue samples were preserved in 500 µL of DNA/RNA Shield (Zymo) and stored at −80 • C.
We extracted total RNA from each preserved sample, using the Zymo Quick-RNA MiniPrep Plus Kit, following the manufacturer's protocol with minor modifications. In brief, approximately 50 mg of tissue was mechanically homogenized in 250 µL DNA/RNA Shield, and 30 µL of Proteinase K Digestion Buffer and 15 µL Proteinase K were added to each sample. Lysis was carried out overnight at 32 • C for muscle tissue, and at room temperature (approximately 21 • C) for liver tissue. Following overnight digestion, extractions proceeded following the manufacturer's protocol. Total RNA was eluted in 60 µL nuclease-free water. We quantified and quality-checked the extracted RNA using a Nanodrop 2000 spectrophotometer and visualized it on a 1% agarose gel. Samples with low A260/A230 ratios (<1.50) were purified using the Zymo RNA Clean and Concentrator Kit, following the manufacturer's protocol. Purified total RNA was sent for library preparation and Tag-Seq sequencing at the University of Texas at Austin Genomic Sequencing and Analysis Facility. Tag-Seq is an efficient alternative to traditional RNA-Seq, in which the 3 end of the input mRNA is preferentially sequenced, rather than the entire molecule (Meyer et al., 2011). Libraries were sequenced on an Illumina NovaSeq6000, generating 100-bp single-end reads.

Reference Transcriptome Generation
To generate a reference transcriptome for L. rhomboides, we prepared RNA-Seq libraries from liver, muscle, gill, fin, brain, heart, and eye tissue from a single pinfish collected from Perdido Pass, AL, in February 2020. We extracted total RNA using the Zymo Quick RNA Mini Prep Plus Kit, following the manufacturer's protocol, with the exception that all samples were incubated for 30 min at room temperature during the lysis step. Extracted total RNA was sent to Novogene Corporation (Beijing, China) for library preparation and sequencing (150-bp paired-end, Illumina NovaSeq6000). Raw reads were quality-checked using FastQC ver. 0.10.1 before and after adapter trimming (Andrews, 2010). Low quality sequences were removed, and adapters were trimmed using Trimmomatic ver. 0.38 (Bolger et al., 2014), with the following parameters: ILLUMINACLIP:Trimmomatic_adapters.fa:2:35:10 LEADING:30 TRAILING:30 MINLEN:25. We assembled these trimmed paired-end reads into a preliminary transcriptome using Trinity ver. 2.9.1 (Grabherr et al., 2011), setting the minimum assembled contig length to 300 bp to eliminate short contigs. Contigs were clustered based on sequence similarity using the program CD-HIT-EST ver. 4.8.1 (Li and Godzik, 2006) to eliminate redundancy (parameters: -c 0.95 -n 10 -d 200). We assessed the completeness of the final transcriptome using BUSCO ver. 4.1.2 (Simão et al., 2015) comparing to the actinopterygii_odb10 database with a minimum BLAST e-value of 1e−6.
We annotated the assembled transcriptome using a series of iterative BLAST searches against both nucleotide and protein databases. Using the program BLAST+ ver. 2.6.0 (Altschul et al., 1990;Camacho et al., 2008), the assembled transcriptome was searched for sequence similarity against the unspliced gene set of the closely related gilthead seabream (Sparus aurata) on Ensembl (Howe et al., 2021 Ensembl release 104 1 ), using the command BLASTn, with a minimum e-value of 1e−10. We filtered matches based on quality, removing any matches that had a query coverage ≤50% or a percent identity <75%. The best match for 1 www.ensembl.org each transcript was identified based on the highest bit-score, using a custom bash script. Transcripts with no high-quality match were isolated from the full transcriptome using SeqKit ver. 0.10.1 (Shen et al., 2016). This process was repeated iteratively for remaining transcripts without annotation information, using (in the following order) the large yellow croaker (Larimichthys crocea), Nile tilapia (Oreochromis niloticus), and zebrafish (Danio rerio) Ensembl unspliced gene sets as references. Any remaining unannotated contigs were then searched using BLASTx with a minimum e-value of 1e−10 against the UniProt Swiss-Prot and TrEMBL protein databases (The UniProt Consortium, 2021). 2 We removed low-quality matches that had query coverage ≤48%, percent identity <75%, or matches to non-chordates, and the best match for each transcript was again identified using a custom bash script. The full pipeline for assembly and annotation of the transcriptome is available online at: https://github.com/kmeaton/ Pinfish_Transcriptome.

Tag-Seq Read Processing and Differential Gene Expression Analysis
We analyzed the raw Illumina Tag-Seq reads for both juveniles and adults separately, using the same pipeline. Reads were quality-checked using FastQC ver. 0.10.1 (Andrews, 2010). The 5 leader sequence was clipped, and duplicated reads with the same header and the same first 20 bases of sequence were removed using a custom Perl script. Illumina adapters were trimmed using the program cutadapt v2.10 (Martin, 2011; parameters:a AAAAAAAA -a AGATCGG -q 15 -m 25 -j 8). Trimmed reads were again quality-checked using FastQC, and then mapped to the L. rhomboides reference transcriptome with Bowtie 2 ver. 2.2.9 (Langmead and Salzberg, 2012;parameters: -local -p 8 -nohd -no-sq -no-unal -k 5). We used custom scripts to generate counts of aligned reads mapping to each transcript for each sample, and to compile read counts for each sample into a single table. An overview of the read-processing pipeline, as well as all the scripts used are available at: https://github.com/z0on/tag-based_RNAseq.
We separated read count data by life stage (juvenile/adult) and tissue type (liver/muscle). Patterns of differential gene expression for each of the four datasets were analyzed separately, using the package DESeq2 ver. 1.28.1 (Love et al., 2014) in R ver. 4.0.2 (R Core Team, 2020). For each of these comparisons, the data was analyzed using two separate tests: the likelihood ratio test (LRT) and the Wald test. Differentially expressed genes (DEGs) were identified as those with an adjusted p-value < 0.05 based on a Benjamini-Hochberg correction. The LRT was used to identify genes that were significantly differentially expressed based on the effect of temperature, controlling for whether a fish had undergone respirometry (adults). To visualize the results of this test, we transformed the count data matrix using a variance stabilizing transformation. We extracted the transformed count data for genes that were significantly differentially expressed based on the effect of temperature (p-adj < 0.05), and these data were used to generate a principal component analysis (PCA) plot of expression patterns. The Wald test was used to identify significantly DEGs for each of the following pairwise comparisons of temperature treatments (controlling for the effect of respirometry in adults): +3 • C vs. control treatment, +6 • C vs. control treatment, and +6 • C vs. +3 • C treatment. The full R script used to identify DEGs is available at: https://github.com/ kmeaton/Pinfish_Transcriptome. Significantly enriched Gene Ontology (GO) categories for each of these pairwise comparisons were identified with a Mann-Whitney U test, following the pipeline provided by Wright et al. (2015) 3 . Significantly upregulated or downregulated GO categories were identified as those with a false discovery rate (FDR) < 0.05.
For both Tag-Seq and the transcriptome reads, sequencing data was deposited in the NCBI Sequence Read Archive (SRA) and Transcriptome Shotgun Assembly (TSA) under BioProject number PRJNA776622.

Oxygen Consumption
For both juvenile and adult pinfish, we observed a significant increase in routine oxygen consumption with increasing temperature (Figure 2). Interestingly, per-unit mass increases in RMR were far higher for juveniles than for adults across the same gradient of temperature change (Figure 2), although proportional increases in oxygen consumption (i.e., Q 10 ) with increasing temperature were similar between the two life stages.
For juvenile pinfish, AICc showed that the best-fit model was a linear model which included temperature and fish mass as covariates. Though the length of exposure varied from 8-11 days, this factor was not significantly associated with changes in routine oxygen consumption and including it as a covariate did not significantly improve the fit of the model. We observed a 24.58 mgO 2 /kg/h increase in RMR for each 1 • C increase in temperature (13.72-35.45 mgO 2 /kg/h, 95% CI, p = 4.89 × 10 −5 , R 2 = 0.36) (Figure 2). On average, pinfish in the +3 • C treatment had 73.74 mgO 2 /kg/h (41.15-106.34 mgO 2 /kg/h, 95% CI) higher routine oxygen consumption than pinfish in the control group, and pinfish in the +6 • C treatment had an average of 147.49 mgO 2 /kg/h (82.29-212.68 mgO 2 /kg/h, 95% CI) higher routine oxygen consumption than fish in the control group.
For the adult pinfish, AICc showed that the best-fit model was a linear model with temperature as the independent variable. Even though mass-specific respiration rate is known to decrease with increased fish mass, including mass as an additional parameter in our model did not significantly improve the fit of the model, and therefore, we excluded this variable from our analysis of adults. Likewise, including duration of exposure as a covariate also did not improve the fit of the model, nor was it significantly associated with changes in routine oxygen consumption. We observed a 10.01 mgO 2 /kg/h increase in RMR for each 1 • C increase in temperature (3.55-16.47 mgO 2 /kg/h, 95% CI; p = 0.00725, R 2 = 0.61; Figure 2). In other words, pinfish in the +3 • C treatment had on average 30.03 mgO 2 /kg/h (10.65-49.40 mgO 2 /kg/h, 95% CI) higher 3 https://github.com/z0on/GO_MWU FIGURE 2 | Routine metabolic rate (mgO 2 /kg/h) of juvenile and adult pinfish. Data for the adult pinfish are shown in hashed boxplots, and data for the juveniles are shown in solid boxplots. Boxes are colored by experimental treatment, where blue corresponds to the control treatment, yellow corresponds to the +3 • C treatment, and red corresponds to the +6 • C treatment.
routine oxygen consumption than pinfish in the control group, while fish at +6 • C had an average of 60.06 mgO 2 /kg/h (21.30-98.79 mgO 2 /kg/h, 95% CI) higher routine oxygen consumption than pinfish in the control group.
Proportional increases in metabolic costs with elevated temperatures were extremely similar between adults and juveniles. Q 10 calculations showed that adult pinfish increase their RMR 1.816-fold over a 10 • C increase in temperature, whereas juveniles increase their RMR 1.820-fold over a 10 • C temperature increase.

Reference Transcriptome
Illumina sequencing generated 213,136,134 paired end reads, for a total of 31,970,420,100 bp of sequence data. The final assembled transcriptome contained a total of 169,541 contigs, with an average length of 1299 bp and a contig N50 of 2,368 bp (Supplementary Table 1).
Estimates of completeness using BUSCO suggest that the de novo reference transcriptome is of high quality, as it contained the following: 3720 (89.9%) complete orthologs (60.2% complete and single copy, 29.7% complete and duplicated), 116 (3.2%) fragmented orthologs, and just 254 (6.9%) missing orthologs. The annotation process resulted in 86,567 annotated transcripts (51.06% of all assembled contigs). The complete annotation table is available as Supplementary File 1.

Differential Gene Expression
On average, Tag-Seq sequencing resulted in 6,366,537 reads per sample (±1,141,544; SD), for a total of 643,020,218 bp of sequencing data per sample (±115,295,951; SD). After adapter trimming and removal of PCR duplicates, there remained an average of 3,021,545 reads per sample (±582,585; SD), for a total of 263,856,222 bp per sample (±50,413,880; SD). The average read mapping rate to the L. rhomboides de novo transcriptome was 95%.
The number of significantly DEGs for each comparison (i.e., +3 • C as compared to the control treatment, +6 • C vs. control, and +6 • C vs. +3 • C) in each tissue type is shown in Table 1 for the juvenile pinfish, and Table 2 for the adult pinfish. Overall, the number of DEGs in juvenile pinfish was higher than the number of DEGs in adults, even after removal of outliers in the juvenile liver dataset (see below, Tables 1, 2 and Figure 3).
The juvenile liver gene expression results suggest an unusually high number of DEGs between the +3 • C and control groups (1735 DEGs) and +6 • C and control groups (2446 DEGs; Table 1 and Figure 3A). This large number of DEGs was likely due to unexpectedly high within-group variance among the control samples ( Figure 4A). Considering this variance, we reanalyzed the juvenile liver gene expression data removing the outliers, which decreased the number of DEGs across both the +3 • C vs. control (42 DEGs) and +6 • C vs. control (161 DEGs) comparisons and increased the number of DEGs in the +6 • C vs. +3 • C comparison (45 DEGs; Supplementary Table 2). There were few major functional differences in the  processes that were significantly up-or downregulated under elevated temperatures when the outliers were removed. Since the major molecular processes did not change extensively with the exclusion of outliers, the results reported here are based on the full dataset (i.e., with outliers included, see Supplementary Files 2, 3 for the raw results). Results with outliers excluded are also available in Supplementary Files 4, 5. Juvenile muscle gene expression results showed similar high-within group variability in expression patterns, particularly among the samples in the +3 • C group (Figure 4B). However, in this case the distribution of samples was continuous, and it was not possible to identify outliers ( Figure 4B). Interestingly, the four individuals which were clear outliers in the juvenile liver dataset ( Figure 4A) were not identified as outliers in the muscle dataset ( Figure 4B).
In the adult pinfish, we saw far fewer DEGs than in the juvenile pinfish (Tables 1, 2 and Figure 3), and there was far less withingroup variance in expression patterns (Figures 4C,D). In the muscle tissue, the vast majority of DEGs were only significantly differentially expressed in the +6 • C vs. control comparison (85 DEGs; Table 2 and Figure 3D), with relatively few genes being up-or downregulated in the +3 • C treatment. In the liver tissue of adults, the number of DEGs was generally similar across all comparisons (Table 2 and Figure 3C).
While there were specific functional categories that were activated in different tissues or life stages, there were several DEGs that overlapped in their expression pattern under heat stress across several comparisons. Specifically, a gene which encodes the MHC class I alpha domain (trov-UBA) was significantly downregulated in juvenile liver (at both +3 and +6 • C) and adult liver (at +6 • C). Several heat shock proteins (HSPs) were significantly upregulated at the +6 • C treatment across either both life stages or tissue types, including hspa1b (member of the HSP70 family), hsp90ab1 (member of the HSP90 family), and dnaja3a (member of the HSP40 family; Figure 5).
In the +3 • C treatment, GO terms associated with immune function, stress response, cell replication, and DNA damage repair had significant changes in expression (FDR < 0.05) when compared to the control fish, particularly for juveniles. Specifically, immunity-associated GO terms such as "humoral immune response, " "immune effector process, " and "positive regulation of immune response" were significantly upregulated in both juvenile liver and muscle tissue at +3 • C. GO categories associated with cellular stress response, including "regulation of cell death" and "response to stress" were significantly deactivated in juveniles in the +3 • C treatment. A similar trend was observed for GO terms related to DNA replication and cell proliferation (e.g., "pre-replicative complex assembly, " "DNA unwinding involved in DNA replication, " "cell cycle, " and "cell division"). DNA damage repair was also significantly deactivated at +3 • C, as we saw downregulation of GO terms such as "doublestrand break repair via break-induced replication." Interestingly, we saw very few significantly enriched GO term categories for adult pinfish at +3 • C. We did observe significant downregulation of ubiquitin B (ubb), a gene coding for a protein involved in degradation of misfolded proteins, in adult muscle tissue at +3 • C. However, this same gene was significantly upregulated in juvenile muscle tissue at +3 • C, which may suggest different mechanisms of response to heat stress across different life stages.
Specific molecular categories that were enriched in the +6 • C treatment (as compared to the control group, FDR < 0.05) were development, transmembrane ion transport, immune system function, muscle contraction, cell replication, and DNA damage repair. We saw significant upregulation of GO terms associated with development for fish at +6 • C when compared to the control group. Specifically, "forebrain neuron development" was activated in the juvenile liver under heat stress, which may suggest changes to nervous system function at elevated temperatures. GO categories related to transmembrane ion transport, such as "extracellular ligand-gated ion channel activity, " "voltage-gated calcium channel activity, " and "potassium ion transmembrane transporter activity" were also significantly upregulated at +6 • C across juvenile and adult liver. In contrast to our findings in the +3 • C treatment, GO terms associated with immune system function (such as "antigen processing and presentation, " "humoral immune response, " "innate immune response, " and "immune effector process") were heavily downregulated at +6 • C in juvenile liver and muscle, as well as in adult liver. In the adult muscle tissue, however, we saw significant upregulation of immune system GO categories (e.g., "positive regulation of immune response, " "humoral immune response, " and "immune effector process"). Categories related to muscle structure and function were also significantly downregulated at +6 • C, particularly in juveniles, as we saw deactivation of GO terms such as "muscle filament sliding, " "actin-mediated cell contraction, " and "muscle myosin complex." We also saw significant deactivation of GO categories associated with DNA replication and cell division (e.g., "prereplicative complex assembly, " "DNA recombination, " and "cell division"), particularly in juveniles, which parallels our findings at +3 • C. GO term categories involved in DNA damage repair, such as "double-strand break repair via break-induced replication, " "cellular response to DNA damage stimulus, " and "site of DNA damage" were also significantly downregulated.
Comparisons among the +3 and +6 • C treatments showed significant differences in processes such as immune system function, muscle activity, oxygen transport, nervous system function, and transmembrane transport (FDR < 0.05). Several of these functional categories were also significantly differentially expressed between the +6 • C and control groups, providing additional evidence that changes in these functions are specific to the highest temperature treatment. When compared to the +3 • C treatment, immune system function was significantly downregulated in the +6 • C treatment for juvenile muscle and adult liver (however, "humoral immune response" was upregulated at +6 • C in the adult muscle). Muscle contraction was also significantly downregulated at +6 • C in juveniles (e.g., "muscle filament sliding, " "actin-mediated cell contraction"). The GO terms "oxygen binding" and "gas transport" were also strongly deactivated in adult liver tissue at +6 • C in comparison to the +3 • C group, which may suggest an inability to cope with elevated oxygen demands in fish exposed to the +6 • C thermal stressor. Nervous system function and transmembrane ion transport were both upregulated at +6 • C as compared to +3 • C, which we also found when comparing the +6 • C group to the control group. Specifically, nervous system-associated GO terms such as "neuron differentiation, " "GABA receptor activity, " and "synaptic signaling" were significantly activated in juvenile and adult liver tissue, suggesting that neural processes are influenced by warming in both life stages.

DISCUSSION
Pinfish are key members of coastal and offshore ecosystems of the western Atlantic, including the Gulf of Mexico, yet there is a gap in our understanding of how this species responds to ocean warming. Increases in RMR were observed in both juvenile and adult pinfish during a simulated MHW, and we FIGURE 5 | Heatmap of variance stabilized expression data for three heat shock proteins. The values for each gene represent the average of the variance stabilized data for all samples in a specific treatment. Note that juvenile liver, adult liver, juvenile muscle, and adult muscle datasets were transformed separately, so comparisons of expression levels for a particular gene should be made within rows, and not columns, of the matrix. Each gene has an individual color scale, where darker reds indicate higher levels of expression, and lighter yellows indicate lower levels of expression. also saw changes in expression of genes associated with gas transport, musculoskeletal system function, cell division, nervous system function, and immunity. The results of this study indicate that L. rhomboides is susceptible to both moderate and extreme warming, and that juvenile pinfish may be more vulnerable to these stressors than adults. Given the ecological relevance of this species, the observed susceptibility of the pinfish to acute warming may have detrimental consequences for other species of the region.

Elevated Routine Metabolic Rate Under Thermal Stress
Results of the respirometry analyses indicate that both juvenile and adult pinfish had significant increases in RMR under elevated temperatures, consistent with the results of previous studies (Wohlschlag et al., 1968;Cameron, 1969;Wohlschlag and Cech, 1970). While the observed fold-change in oxygen consumption (i.e., Q 10 ) was the same for juveniles and adults over a six-degree temperature increase from 21 to 27 • C, juveniles had far more extreme increases in per-kilogram metabolic rate than adults (Figure 2). Considering the elevated per-kilogram oxygen demand that we saw in juveniles, exposure to acute temperature stressors such as MHWs is more metabolically costly for juveniles than adults. These heightened metabolic costs also had consequences at the molecular level, as energy needs of the organism shift away from growth and move toward maintenance (see section "Molecular Responses to Acute Warming in the Pinfish").
One of the limitations of our study is that maximum oxygen consumption was not measured, and aerobic scope could not be estimated. However, the observed increases in metabolic demands and the physiological sensitivity of juveniles to thermal stress aligns well with the natural history of the pinfish. Juvenile abundance declines in hot, shallow, nearshore areas during the warmest months of the summer (Hellier, 1962;Adams, 1976;Darcy, 1985), which may be due to the high metabolic costs associated with exposure to warm waters.

Molecular Responses to Acute Warming in the Pinfish
While molecular responses to heat stress were often tissueand/or life stage-specific, certain functional changes overlapped across multiple comparisons. Notably, HSPs (e.g., hsp90ab1 and hsp90b1) were activated at +3 • C in juvenile pinfish, and at +6 • C in both juvenile and adult pinfish (e.g., dnaja3a, hsp90aa1, and hspa1b; Figure 5). These highly conserved proteins are hallmarks of cellular stress across eukaryotes, and have a wide array of functions, including protein folding and stabilization under conditions of stress (Welch, 1993;Iwama et al., 1998;Feder and Hofmann, 1999). Activation of HSPs in juveniles at a moderately elevated temperature (24 • C, the +3 • C treatment) further confirms that juveniles are more susceptible to warming than adults. Exposure to 27 • C (the +6 • C treatment temperature) elicited HSP upregulation in both life stages, indicating that conditions only slightly above the recorded seasonal maximum for the collection site will have detrimental consequences for the species as a whole.
The heightened metabolic demands that both juvenile and adult pinfish experience during acute heat stress lead to changes in molecular function and re-allocation of energy toward homeostatic maintenance and away from growth. Specifically, cell division and replication were consistently downregulated at both +3 and +6 • C, when compared to the control group.
We saw downregulation of cyclin-A2 (ccna2), which promotes transition through G1/S and G2/M phases of the cell cycle (Pagano et al., 1992), as well as polo like kinase 1 (plk1), which is involved in mitotic spindle formation and function (Golsteyn et al., 1995). Similar deactivation of genes involved in growth and proliferation under warming has been observed in the muscle tissue of the goby Gillichthys mirabilis and the cyprinid Squalius torgalensis (Buckley et al., 2006;Jesus et al., 2016). This deactivation of growth and proliferation may be associated with changes in energy availability and increased metabolic demands during thermal stress (Figure 2), resulting in the re-allocation of resources toward maintaining homeostasis and repairing damaged molecules (Buckley et al., 2006).
Similarly, expression of genes associated with muscle contraction and function was significantly downregulated at +6 • C in both juveniles and adults. We saw significant downregulation of genes encoding the cardiac myosin heavy chain 6 (myh6) and troponin T (tnni2a.3) at +6 • C, which are key proteins involved in cardiac and skeletal muscle function, respectively (Perry, 1998;Abu-Daya et al., 2009). These changes to muscle contraction under elevated temperatures may have implications for swimming speed and performance of the pinfish. For example, studies on rainbow smelt (Osmerus mordax) and rainbow trout (Oncorhynchus mykiss) acclimated to warm temperatures show declines in maximum swimming speed and slower muscle contraction (Woytanowski and Coughlin, 2013;Coughlin et al., 2016;Shuman and Coughlin, 2018). Decreases in swimming capacity under elevated temperatures can have cascading effects, as they could be involved with changes in predator/prey interactions and dispersal potential of this ecologically relevant species.
The molecular response to warming also differed based on the intensity of the heat stressor, as we observed differential expression of genes associated with nervous system function and oxygen transport between the +3 and +6 • C treatments. Downregulation of GO categories associated with oxygen transport was only observed at +6 • C, suggesting that the ability of pinfish to supply oxygen to tissues declines as warming increases. Molecular challenges to oxygen transport under heat stress have been previously reported in catfish (Liu et al., 2013), and it has been proposed that cardiovascular function becomes impaired at elevated temperatures (Pörtner, 2001;Pörtner and Knust, 2007;Pörtner et al., 2017; but see Jutfelt et al., 2018). Changes to nervous system function were also unique to the +6 • C group, indicating that there may be systemic changes to transmission of neuronal signals at extremely elevated temperatures. Considering that both of our experimental temperatures represent realistic MHW scenarios that may be experienced by pinfish in the wild today, these results highlight the importance of studying the differential responses based on the intensity of warming.

Immune Response Under Elevated Temperatures
Immune function of the pinfish was significantly modified at both +3 and +6 • C, though the intensity of warming caused different effects. Generally, immune response was activated in the +3 • C treatment and downregulated in the +6 • C treatment. The complex relationship between temperature and immune function has also been reported in previous studies. Elevated expression of immune-associated genes under heat stress has been reported in freshwater (O. mykiss; Lewis et al., 2010;Rebl et al., 2013;Verleih et al., 2015), and marine fishes (Pomacentrus moluccensis; Kassahn et al., 2007; Acanthochromis polyacanthus, Ostorhinchus cyanosoma, Ostorhinchus doederlini, Cheilodipterus quinquelineatus; Bernal et al., 2020). However, other studies have found declines in immune function at elevated temperatures in turbot (Scophthalmus maximus; Zhang and Sun, 2017) or conflicting results in sticklebacks (Gasterosteus aculeatus; Dittmar et al., 2014) and carp (Ctenopharyngodon idellus; Yang et al., 2016). Ultimately, the interaction between water temperature and immunity appears to be complex, and there is no single clear effect of temperature on immunity that is common to all species (Bowden, 2008;Tort, 2011). However, given the extreme deactivation of immune-associated genes in pinfish in the +6 • C treatment, these results suggest warming results in higher vulnerability to disease in this species.

Consequences of Warming Across Multiple Tissues and Life Stages
Gene expression analyses revealed major differences in functional responses to heat stress by tissue type. Specifically, changes in oxygen binding and gas transport were only observed in the liver, and downregulation of muscle contraction was only observed in the muscle. Changes in these relevant functions would have been missed if only one tissue type had been examined. Since studies of gene expression in the context of warming usually examine only a single tissue type, typically one that is closely linked to the phenotype being evaluated (Kassahn et al., 2007;Lewis et al., 2010;Rebl et al., 2013;Verleih et al., 2015;Qian and Xue, 2016;Yang et al., 2016;Zhang and Sun, 2017;Bernal et al., 2020;Huang et al., 2020) we recommend the inclusion of more than one tissue type whenever possible. This will enhance our understanding of the potential effects of warming at the organismal level.
The differences in both oxygen consumption and gene expression between juvenile and adult pinfish highlight the importance of analyzing the effects of heat stress on different life stages. Increases in routine metabolic demand were more severe over the same thermal gradient in juvenile pinfish, which coincided with the activation of stress-response associated genes (including HSPs) at a lower temperature threshold than the adults. Furthermore, the total number of DEGs was far higher for juvenile pinfish than for adults, across nearly all comparisons (Tables 1, 2), indicating that juveniles must change their patterns of gene expression more than adults to deal with the same temperature change. The combination of these results shows that juvenile pinfish are more susceptible to warming than adults, once again highlighting that different life stages have different vulnerabilities to the same stressor (Dahlke et al., 2020).
Overall, the observed thermal sensitivity of juvenile pinfish aligns well with the natural history of the species, as youngof-the-year pinfish are highly abundant in shallow estuaries and eelgrass beds during late winter, spring, and early summer (Darcy, 1985). Pinfish abundance in nearshore areas generally declines during mid-late summer (i.e., July-August; Hellier, 1962;Adams, 1976), during which time these habitats generally attain their maximum yearly temperature. While it is difficult to ascertain the exact driver of the movement of pinfish during the warmest months of the year (e.g., developmental cues, predation, environmental factors, or a combination thereof), it is possible that juvenile pinfish may be emigrating from shallow areas to reduce the costly metabolic demands imposed by seasonal temperature extremes. The results of our study suggest that the susceptibility of juvenile pinfish to rising temperatures presents a major concern, as sea surface temperatures in the northern Gulf of Mexico are projected to continue to rise by 2-2.5 • C by the end of the century (Biasutti et al., 2012). These rising temperatures can cause changes in molecular functions and energetic demands at the organismal level, which may lead to behavioral modifications to avoid exposure to costly environmental conditions (Muñoz and Bodensteiner, 2019). Rising sea surface temperatures may cause juvenile pinfish to egress from nearshore ecosystems earlier than normal, which could produce changes in seasonal abundance, potentially influencing coastal food webs.
Ultimately, the results of this study provide evidence that the ecologically important L. rhomboides may be more susceptible to warming in the Gulf of Mexico and western Atlantic than previously assumed. Physiological and molecular analyses show that juvenile pinfish have significant increases in metabolic demand under elevated temperatures, and they show differences in molecular pathways related to cell growth and proliferation, immune response, nervous system function, and muscle contraction. Adults showed similar functional changes in gene expression but had milder responses with fewer total DEGs. Despite the clear costs associated with exposure to elevated temperatures in this species, juvenile and adult pinfish were able to survive for 8-11 days at the +6 • C treatment. The observed changes to gene expression suggest that plasticity allows pinfish to survive under heat stress, despite the clear energetic costs associated with such responses. Still, the full potential for developmental plasticity and/or adaptation to warming for this species remains to be evaluated. Considering that the northern Gulf of Mexico is projected to warm by several degrees by the end of the century (Biasutti et al., 2012), and MHWs are increasing in frequency and intensity (Frölicher et al., 2018), both juveniles and adults of L. rhomboides could see detrimental consequences with warming in coming decades. The molecular and physiological effects of exposure to elevated temperatures could result in changes in behavior, such as seeking out cooler waters earlier in the spring and changes in seasonal abundance of this ecologically important species, which could have broad implications for the future of coastal ecosystems of the Gulf of Mexico.

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 in the article/ Supplementary Material.

ETHICS STATEMENT
The animal study was reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) of Auburn University.

AUTHOR CONTRIBUTIONS
KE and MB designed the experiment. KE, AH, and MB collected the samples. KE and AH conducted the experiment and collected the data, with advice from JS and MB. KE and MB analyzed the data, with advice from JS. All authors contributed to the drafting and editing of the manuscript.

FUNDING
This study was made possible with funds from the Department of Biological Sciences at Auburn University to MB.