Transcriptional flexibility during thermal challenge corresponds with expanded thermal tolerance in an invasive compared to native fish

Capacity to cope with warming temperatures is a key determining factor of species' persistence under global climate change. Many successful invasive species have heightened thermal tolerance relative to their native counterparts, which may provide competitive advantages for habitat utilization and resource acquisition under warming scenarios, ultimately contributing to radically altered community composition. Enhanced transcriptional plasticity may be an important factor conferring superior abilities to cope with environmental stress, but the molecular mechanisms driving key differences of organismal traits in invasive versus native species are not well known. Although it is predicted that established invaders will evolve canalized phenotypes well‐adapted to new environments, it is not clear whether the same expectations are true for invaders of variable thermal environments or under climate warming scenarios where abilities to cope with fluctuating and increasing temperatures may provide fitness advantages. Here, we compare a highly successful invasive fish and a sympatric endangered native fish living in a dynamic estuarine environment that is projected to warm under climate change. We linked organismal physiological limits with global transcriptional responses at multiple common relative and absolute temperature thresholds and determined that heightened thermal tolerance of invasive Inland Silversides (Menidia beryllina) is associated with transcriptional changes that are greater both in the number of genes and the magnitude of response relative to native Delta Smelt (Hypomesus transpacificus). Modulated genes contributed to the enrichment of biological processes that differed between species and generally increased with temperature. These results are in concordance with the hypothesis that transcriptional plasticity may play a key role in determining population persistence, species interactions, and shaping community assemblages under climate change. Future studies encompassing a wider range of species and taxa are needed to determine whether this is a general pattern found between native and invasive species more broadly.


| INTRODUC TI ON
Rapid changes in environmental conditions due to climate change are having profound effects on biodiversity across biological scales and are forecasted to continue increasing in coming decades (Hoegh-Guldberg & Bruno, 2010;Hughes, 2000;Parmesan & Yohe, 2003). Temperature plays a central role in determining organismal performance, fitness, and geographic distributions, such that climate warming is modifying selective pressures for many populations (Hochachka & Somero, 2002;Pörtner, 2002;Pörtner & Farrell, 2008). Survival and successful reproduction to maintain population viability under climate warming will largely depend on existing phenotypic traits such as thermal tolerance, and organisms' capacity to respond to thermal regime shifts via phenotypic plasticity (Bradshaw, 1965), rapid evolutionary adaptation, or geographic range and phenological shifts (Brown, O'Connor, et al., 2016;Hoffmann & Sgro, 2011;Hofmann & Todgham, 2009).
Thermal tolerance is a key trait frequently subjected to strong natural selection and has been correlated with thermal regimes at species and population levels (Clusella-Trullas et al., 2011;Eliason et al., 2011;Fangue et al., 2006;Stillman, 2002;Sunday et al., 2011). Both temperature tolerance and temperature-dependent performance can differ across ontogenetic stages and thermal history conditions (Komoroske et al., 2014;Schulte et al., 2011), and thermal tolerance breadth is often proportional to the magnitude of in situ temperature variation (Addo-Bediako et al., 2000;Deutsch et al., 2008;Ghalambor et al., 2006;Janzen, 1967). Additionally, the magnitude and direction of plastic responses to temperature changes can evolve in response to selection (Reed et al., 2011;Schlichting & Smith, 2002). These observations suggest that capacity for phenotypic plasticity may confer fitness advantages in variable and rapidly changing environments and will play a central role in species' abilities to respond to environmental warming (Calosi et al., 2008;Helmuth, 2009;Stillman, 2003). Divergent population-specific patterns of thermal plasticity have been documented in species inhabiting different in situ thermal environments (Narum & Campbell, 2015), and the costs of maintaining plastic capacity may only be outweighed by benefits in certain environmental contexts (Hendry, 2016). Theoretical work has also suggested that phenotypic plasticity can facilitate evolutionary adaptation via genetic accommodation (Levis & Pfennig, 2016), serving as an important source of adaptive capacity of environmental stress tolerance (Ghalambor et al., 2015;Schneider & Meyer, 2017). Thus, understanding differences in phenotypic plasticity among species and populations is important to understanding both proximate and long-term biological responses and resilience to changes in thermal regimes.
Changes in environmental conditions also alter species interactions, community assemblages, and ecological function (Gilman et al., 2010;Hooper et al., 2005). The devastating impacts of introduced species on native biodiversity and ecological regime shifts are well-documented around the world (Molnar et al., 2008;Vitousek et al., 1997), and increased thermal tolerance is associated with invasion success in both freshwater and marine ecosystems (Bates et al., 2013;Kolar & Lodge, 2002). Thus, in addition to direct negative impacts on native populations, this thermal advantage can result in a competitive edge for eurythermal nonindigenous species in habitat utilization or resource acquisition (Carveth et al., 2006;Cheng et al., 2016;Sorte et al., 2013;Stachowicz et al., 2002) and, in turn, have indirect impacts at community and ecosystem levels (Hendry, 2016). Combined abiotic-biotic stressors under global climate change may further exacerbate sublethal impacts of thermal stress on less tolerant native species under warming scenarios (Dukes & Mooney, 1999). Although such trends have been described across a diversity of taxa and ecological contexts and it is thought that physiological capabilities are key factors in the successful transport, establishment and spread of non-native species (Bates et al., 2013;Kelley, 2014), the underlying mechanisms driving these patterns are not well understood. Invasion biology has largely focused on factors such as propagule pressure (i.e., the intensity and frequency of arriving recruits), while treating species' traits such as physiological tolerances as fixed. However, understanding the extent to which plasticity aids colonization into new environments, as well as dealing with in situ environmental change, remains a key question in evolution and ecology (Hendry, 2016) and there have been recent efforts to consider the roles of plasticity and evolutionary potential when assessing invasion risk (Briski et al., 2018;Wellband & Heath, 2017;Whitney & Gabler, 2008). There is also evidence that natural selection and rapid adaptation of non-native species can result in the establishment of "pre-adapted" individuals better poised to colonize and spread in novel environments, even despite expected reductions in genetic diversity in founding relative to source populations (Briski et al., 2018;Whitney & Gabler, 2008). This likely includes selection for high plasticity especially under harsh or variable thermal conditions frequently experienced during the transport phase and establishment in novel environments (Ghalambor et al., 2007).
Enhanced phenotypic plasticity has been reported as a critical factor in the post-establishment range expansion of a highly successful invader (Wellband & Heath, 2017). However, general support for the hypothesis that invaders are more plastic than native species has been inconsistent, perhaps as a result of contextual and methodological differences in how plasticity is assessed (Davidson et al., 2011;Godoy et al., 2011;Palacio-López & Gianoli, 2011).
It has been predicted that plasticity will be greatest during early  (Lande, 2015) where there is little benefit to maintaining a highly plastic phenotype that is energetically costly. This makes it challenging to control for temporal effects in empirical comparisons of plasticity in native versus invasive species. However, it is not clear whether the same expectations hold for invaders of variable thermal environments or under climate warming scenarios where capacity to cope with dynamic and increasing acute and chronic temperatures may continue to provide fitness advantages. There is empirical support that greater thermal plasticity can evolve under conditions of higher environmental variability relative to stable ancestral conditions (Morris et al., 2014) and that it can be maintained in spatial and temporal heterogeneous environments (Gianoli & Gonzalez-Teuber, 2005;Hendry, 2016). Additionally, tools to measure global changes at the molecular level that complement whole-organism metrics of phenotypic plasticity have only recently become available for most nonmodel organisms, without which it has been challenging to assess whether differences (or lack thereof) detected at the organismal level are driven by compensatory physiological responses (i.e., alterations at the molecular level). In order to forecast when invasive species are likely to establish and outcompete native species under climate warming, studies are needed that not only assess differences in thermal traits at the organismal level, but also link these with the underlying genomic mechanisms driving these patterns.
Transcriptional modulation of genes in key biological pathways is one way that organisms can respond to new or changing environmental conditions and may be at the heart of capacity for phenotypic plasticity (Aubin-Horth & Renn, 2009;Schlichting & Smith, 2002) and underpin responses to environmental stressors and adaptive resilience (Logan & Buckley, 2015;Sandoval-Castillo et al., 2020;Whitehead, 2012). Linking molecular changes with physiological phenotypes (e.g., thermal tolerance and performance) can provide insight into the mechanisms underlying an organism's capacity to cope with environmental stress (Feder & Hofmann, 1999). For example, Wellband and Heath (2017) found evidence that adaptive transcriptional responses associated with regaining homeostasis under thermal stress correlated with greater invasion success among two nonindigenous freshwater fishes, and Jeffries et al. (2016) linked transcriptional responses with metabolic performance and thermal tolerance among two native fishes to examine differential climate warming sensitivity. Differential transcriptional plasticity (i.e., changes in the transcriptional profile in response to different thermal environments) to thermal stress can be heritable (McCairns et al., 2016;Tedeschi et al., 2016) such that selection can act upon these attributes. For example, Sandoval-Castillo et al. (2020) recently found evidence that adaptive transcriptional plasticity has evolved in warmer climates in Australian rainbowfish, contributing to greater thermal resilience in the subtropics. Moreover, plasticity can evolve rapidly over contemporary timescales, particularly for labile traits such as gene and protein expression, to potentially drive adaptive evolutionary responses under climate warming (Hendry, 2016;Oomen & Hutchings, 2017). It is also plausible that eurythermal invasive species have evolved greater capacity for transcriptional plasticity, underlying superior abilities to cope with thermal stress relative to mesothermal natives under climate change (Komoroske et al., 2015). Comparing transcriptional thermal plasticity among these species and assessing differences in the magnitude and identity of genes involved in these responses may offer key insight into the cellular and biochemical pathways through which thermally tolerant invaders can outcompete native species, and how expectations may contrast under different environmental contexts (Hendry, 2016;McCairns & Bernatchez, 2010;Studivan & Voss, 2020).
Environmental temperature is critical for aquatic ectotherms because their body temperatures typically conform quickly to water temperatures (Beitinger et al., 2000). Collectively, teleost fishes are the most diverse group of vertebrates and inhabit nearly every aquatic system on the planet (Nelson, 2006). Introduced species pose a strong threat to native fish biodiversity and their ecological communities; almost 70% of fish extinctions in North America in the 20th century have been associated with invasive species (Miller et al., 1989;Vitousek et al., 1997). Thermal physiology plays a critical role in the establishment of invasive ectotherms, including broader thermal tolerance and higher transcriptional plasticity relative to native species (Kelley, 2014). Complex and differential evolutionary histories among fishes have altered the "genomic tools" with which fishes can respond to contemporary thermal stress, resulting in some species more successfully responding to environmental change relative to others (Oomen & Hutchings, 2017;Schulte, 2004).
Historically, the majority of our knowledge of the mechanisms underlying thermal plasticity and adaptation in fishes has come from a few well-studied eurythermal or stenothermal species, but recent work capitalizing on advancements in genomic technologies in nonmodel organisms has begun to expand our understanding in a wider diversity of species with a broad range of physiological adaptations, life histories, and ecological contexts (Connon et al., 2018). Research to date has supported conserved global coordination of key pathways and biological functions, combined with species-specific adaptations (Oomen & Hutchings, 2017). There are relatively few empirical data that directly assess the mechanisms of thermal tolerance and plasticity on sympatric native and invasive fishes. Such studies inherently deviate from a traditional comparative framework because native and invasive species are often phylogenetically distant from one another, and native species of high conservation concern may be endemic or have other constraints (e.g., lack of availability of multiple populations or sister species for comparisons). However, despite these limitations, they provide important opportunities to understand whether the same affected biological pathways and functions are broadly conserved and, in conjunction with organismal trait information, can inform assessment of the species' abilities to persist under climate warming.
The Delta Smelt (Hypomesus transpacificus) is an endemic fish in the San Francisco Estuary, California USA, and has declined since the 1980s and precipitously since the early 2000s (Brown et al., 2009;Feyrer et al., 2007). The underlying causes of Delta Smelt decline along with several other forage fishes have been potentially attributed to a multitude of environmental stressors including habitat degradation, temperature stress, and interspecific interactions with invasive species (Brooks et al., 2012;Feyrer et al., 2007;Sommer et al., 2007;Winder & Jassby, 2011). Delta Smelt are mesothermal fish with moderate thermal tolerance and limited acclimation capacity (Komoroske et al., 2014(Komoroske et al., , 2015, and temperature increases due to climate warming are predicted to further reduce the availability of adequate habitat in the absence of rapid thermal adaptation (Brown, Komoroske, et al., 2016). Delta Smelt have reduced transcriptional capacity to respond to heat stress relative to well-studied eurythermal species (Komoroske et al., 2015). Conversely, the Inland Silverside (Menidia beryllina) is an invasive fish in the San Francisco Estuary whose abundance has steadily increased since introduction to the system in 1967 (Cook & Moore, 1970). Characteristic of successful invasive species, Inland Silversides have high fecundity and are thought to be relatively stress tolerant (Lutterschmidt & Hutchison, 1997;Moyle, 2002;Swanson et al., 2000). Native to the Atlantic coastline of North America and the Gulf of Mexico that is the strongest temperature gradient in the world (Baumann & Doherty, 2013), it is likely that this evolutionary history contributes to their current eurythermal capabilities (Hendry, 2016). Both species are small-bodied, primarily annual fishes with sympatric distributions in the San Francisco Estuary. In addition to competing with Delta Smelt for habitat and resources, Inland Silversides also prey upon Delta Smelt especially in offshore habitats (Baerwald et al., 2012;Schreier et al., 2016).
This offers an ideal system to explore the underlying molecular drivers of enhanced plasticity and thermal tolerance of an invasive versus a native fish competing in a variable environment.
RNA sequencing coupled with organismal phenotype metrics offers a powerful method to examine environmental responses in natural populations (Alvarez et al., 2015;Connon et al., 2018;Oomen & Hutchings, 2017). Integrating these data can be particularly useful in sensitive species where traditional methods are challenging, contributing to effective species management by determining critical thresholds for optimal growth or stress responses (Connon et al., 2018). The overarching goal of our work was to understand the molecular mechanisms underlying thermal tolerance, and how these may contribute to the physiological resilience of a prolific invader to cope with thermal stress and outcompete a critically endangered native fish under future climate warming scenarios. Thermal plasticity can manifest differently across sublethal thresholds and timescales; therefore, we first quantified acute upper thermal tolerance (i.e., the physiological maximum temperature limit) in both species and then examined the effects of temperature on the transcriptional responses at several common temperatures below their physiological limits ( Figure 1). Specifically, we used this coupled thermal tolerance-global transcriptional response data (a) to examine whether thermal stress invokes similar cellular mechanisms that are highly conserved between these phylogenetically distant species and/or whether Inland Silversides' enhanced capacity for tolerating environmental stress is facilitated by transcriptional responses of distinct genes and (b) to assess whether the two species initiate stress responses at similar relative thresholds (e.g., CT Max -2°C) or whether the endangered native shows signatures of a sublethal stress response at lower thresholds indicative of greater sensitivity to climate warming. We hypothesized that if Inland Silversides displayed enhanced thermal tolerance relative to Delta Smelt, they would also display superior abilities to elicit transcriptional responses to cope with thermal stress through the magnitude of their responses or differential biological functions associated with affected genes. Identifying the underlying drivers responsible for differences in thermal tolerance offers insight into the evolution of this critical trait in a declining native versus an abundant, successful invader. Contrasts of these two species sheds light on the potential molecular mechanisms and roles F I G U R E 1 (a) Experimental design of temperature treatments selected based on thresholds relative to (b) critical thermal maximum for each species (mean + SEM of all individuals assayed; n = 20 for Delta Smelt and n = 15 for Inland Silverside). Asterisk denotes the significant difference between the two species (p < 0.0001; linear model, see Table S1) of adaptive thermal plasticity in influencing population persistence, interspecific competition, and, ultimately, fish community responses to climate warming.

| Fish collection, culture, and holding conditions
Juvenile Inland Silversides were collected by beach seine from the Napa River in July 2014 in Napa, California, and then transported to the University of California, Davis Bodega Marine Laboratory, Bodega Bay, California, and held at 16°C and 2.5 psu saltwater in 68-L tanks. Fish were fed live artemia, weaned onto partial pellets and given at least 3 weeks to acclimatize to laboratory conditions prior to experimentation. Species identity was visually confirmed based on morphological diagnostics, and we targeted fish 50-60 mm fork length and classified them as juveniles based on collection timing and previous work correlating age (via otoliths) and body size (Barkman & Bengtson, 1987;Brander et al., 2013). All handling, care, and experimental procedures used were reviewed and approved by the UC Davis Institutional Animal Care and Use Committee (IACUC Protocol # 16845).
Due to their threatened status, Delta Smelt cannot ethically be collected from the wild; therefore, fish were obtained from a refuge population   Spectrum Brands, Inc.) and fed an ad libitum pellet and plankton mixture. Environmental acclimation conditions (particularly temperature and salinity) were therefore very similar for the two species. We note that due to the specific needs of Delta Smelt culturing , there were some differences between other holding conditions of the two species, specifically in the exact food composition, holding tank size, and experimental timing. However, given the strong physiological and transcriptomic effects of temperature (reviewed in Logan & Buckley, 2015) and the consistency of thermal tolerance measures in these species (Davis et al., 2019;Jeffries et al., 2016), these differences likely exert small effects on our results relative to temperature. Additionally, instead of direct comparisons of gene expression levels between the species, we focus our objectives and analyses on contrasting the species' responses to thermal challenges that should be largely robust to the differences between the two experiments (i.e., contrasts of treatments vs. handling controls within each species, quantifying the slopes of the reaction norms rather than the intercepts).

| Critical thermal maximum experiments and analyses
The Delta Smelt juvenile thermal tolerance data presented here are a subset of the total dataset across life stages and acclimation temperatures in Komoroske et al. (2014), while all other data presented in this study have not yet been published. We determined Inland Silverside acute temperature tolerance by replicating previous methodology employed to establish thermal tolerance for Delta Smelt (Komoroske et al., 2014). In brief, we used critical thermal methodology (Beitinger et al., 2000), which employs a thermal increase of 0.3°C/min for critical thermal maximum (CT max ) trials to have fish core temperatures closely track changes in water temperature without allowing time for fish to thermally acclimate during the experiments (Becker & Genoway, 1979). We employed loss of equilibrium (LOE; the inability to maintain an upright position in the water column) as the endpoint determining CT max , signifying "ecological death" (the inability to survive in an ecological context; Becker & Genoway, 1979;Beitinger et al., 2000;Cox, 1974). Once LOE was reached, we recorded temperature and immediately returned fish to adjacent chambers containing water at the fish's original acclimation temperature and allowed them to recover. For each CT max trial, we placed a randomly selected fish from the larger holding tanks in a 2 L black chamber filled with water at the acclimation temperature and covered with black mesh. Recovered fish were weighed (wet mass ± 0.1 g) and measured (fork length ± 0.5 mm), and returned to separate holding tanks to ensure they would not be selected for subsequent CT max trials (Delta Smelt n = 20 and Inland Silverside n = 15). We conducted statistical analyses using R (version 2.15.2; R-CoreTeam, 2012) and associated packages (e.g., ggplot2; Wickham, 2009), calculating CT max as the arithmetic mean of the LOE temperatures for each species (Beitinger et al., 2000;Cox, 1974) and assessing the effects of species and fish size on thermal tolerance by constructing a linear model (LOE temperature ~ species*fork length; Table S1). We evaluated data assumptions and model fit graphically (residual versus fitted values, residual versus predictor values, and residual histograms; Zuur et al., 2009; see Data Availability Statement for details and access to all code used for analyses).

| Temperature challenge exposures
Acute thermal exposures followed methodology of Komoroske et al. (2015), with modified temperature treatments selected (a) relative to the acute thermal tolerance for each species and (b) to establish common absolute temperatures between species (see experimental design details in Figure 1). This design allowed us to assess both how the two species might cope with certain environmental conditions in situ (e.g., if summer water temperatures warm to 26°C), but also to understand if similar or different physiological functions occur at comparable thresholds (e.g., do both species exhibit molecular signatures of cellular stress response at CT Max -2°C?). In brief, for each acute thermal exposure, we placed randomly selected fish into a 7.6-L round black chamber filled with water at the acclimation temperature and covered with black mesh. Fish were given a 30-to 45-min habituation period prior to the start of the temperature increase, followed by increasing temperature (0.3°C/min) to the target. When each treatment temperature was reached, fish were removed and either immediately sacrificed and sampled (0 min recovery) or transferred into 35.2-L black round recovery tanks at the original acclimation temperature, and allowed to recover for 60 min prior to sampling. Handling control trials were conducted by subjecting fish to the same conditions, except they were held at the acclimation temperature for the average time of the treatment trials, followed by subsequent transfer into appropriate recovery tanks. Delta Smelt used for these experiments were of the same generation and multifamily group as those used for the experiments in Komoroske et al. (2014), and Inland Silversides were all from the same field collection in 2014 described above.
A total of 5-6 replicates (i.e., individual fish) were used at each acute challenge temperature × recovery time treatment for each species to account for increased biological variation expected in nonmodel organisms (Oomen & Hutchings, 2017;Robles et al., 2012) for a total of 96 fish used for the gene expression work.
Fish were sacrificed with an overdose of tricaine methanesulfonate (MS-222; Finquel), weighed (wet mass ± 0.1 g), and measured (fork length ± 0.5 mm), and gill tissue was immediately dissected and frozen in liquid nitrogen. We chose gill tissue because in teleost fishes it serves critical functions underlying thermal physiology, including oxygen uptake, osmotic and ionic regulation, nitrogenous waste excretion, and other critical organismal functions (Evans et al., 2005;Somero et al., 2017). It has also been extensively demonstrated that gills undergo large changes in gene and protein expression related to phenotypic plasticity in other fishes (Buckley et al., 2006;Logan & Somero, 2011). Samples were stored at −80°C until RNA extraction.

| RNA extraction, library preparation, and sequencing
Total RNA was extracted from gill tissue using Qiagen RNeasy Kits (Qiagen, Inc.) according to manufacturer's instructions. RNA concentrations (ng/µl) and purity (A260/A280 and A260/A230 ratios) were determined using a NanoDrop ND1000 Spectrophotometer (NanoDrop Technologies, Inc.), and integrity was verified through electrophoresis on a Bioanalyzer (Agilent; RIN = 9.5 ± 0.34 SD). Total RNA was then used to prepare cDNA libraries for 96 total samples (n = 5-6 replicates for each temperature × recovery time treatment)

| Sequence data processing and quality control
We demultiplexed and combined sequences across lanes for each sample and used FASTQC to assess initial data quality (see Data Availability Statement for details and access to all code used for analyses). Raw reads per sample were then filtered via adaptor trimming with Scythe (Buffalo, 2014) and adaptive trimming with Sickle (Joshi & Fass, 2011). Divergence between Delta Smelt and Inland Silversides precluded aligning reads for both species to a common reference, so data analyses for each species were conducted separately but following the same methodology until merging at the ortholog level as described below. Trimmed sequences were aligned to previously generated reference transcriptomes for each species, respectively (Jeffries et al., , 2016, using BWA (Burrows-Wheeler Aligner, version 0.7.15-r1140; ). We conducted additional filtering to remove secondary and supplementary alignments to avoid issues of paralog and chimeric sequences. Resulting filtered alignments were processed to counts per transcript contig using SAMtools version 1.3 idxstats . We used Transdecoder (v. 5.5.0; Haas et al., 2013)  Finally, we also generated separate count matrices of filtered transcriptomes for each species to explore responses to thermal stress in genes lacking an orthologues.

| Differential gene expression analyses
To achieve our first specific objective examining whether thermal stress invokes highly conserved cellular mechanisms and/or whether Inland Silversides' higher thermal tolerance is facilitated by differential transcriptional responses, we performed differential expression analyses using the limma R package (Ritchie et al., 2015;v3.30.13).
Prior to running analyses, we removed genes with less than one count per million in five samples and visually inspected pre-and post-filtering count distributions and voom plots of mean-variance relationships to ensure filtered datasets met model assumptions. We also visually confirmed library size consistency across treatments and species and conducted library size normalization within limma prior to implementing models. Because our central research questions were focused on identifying how each species responds to sublethal thermal stress at multiple levels relative to control conditions, we focused analyses on the main effects. We therefore employed the limma-voom transformation with a general linear model as recommended by Ritchie et al. (2015), where gene expression is modeled as a function of the species-temperature treatment (see Data Availability Statement for details). A priori contrasts of interest were designed to assess transcriptional responses at each temperature threshold and time point as a combined treatment relative to the handling control within the species. We also evaluated contrasts of handling controls between the two species to explore differences in baseline constitutive expression; however, some differences in experimental conditions as detailed above limit robust interpretation of these comparisons (Conesa et al., 2016). Genes were considered differentially expressed at a Benjamini-Hochberg corrected false discovery rate (FDR) <0.05 (Benjamini & Hochberg, 1995). Finally, we repeated the analyses as described above for the species-specific count matrices and quantified the number of differentially expressed genes that were private to either Delta Smelt or Inland Silversides (i.e., those that did not have an ortholog identified in the other species). Graphical visualization of differentially expressed genes was conducted within limma, ggplot2 (Wickham, 2009), gplots (Warnes et al., 2019), and UpSetR (Conway & Gehlenborg, 2017).

| Functional enrichment analyses
To achieve our second specific objective assessing if the two species initiate sublethal stress and other biological responses at similar or differential thermal thresholds, we performed functional enrichment analyses. Blast2GO (Conesa et al., 2005) was previously used to identify Gene Ontology (GO) categories associated with annotated genes in the Delta Smelt and Inland Silverside reference transcriptomes (Jeffries et al., , 2016, providing a common metric to compare affected biological processes between the two species. We combined and removed duplicate GO terms for annotated transcripts within orthogroups and then used GO information associated with each orthogroup to assess functional enrichment for a priori contrasts specified above using Fisher's exact tests as implemented in topGO (Alexa & Rahnenfuhrer, 2019) using the weight01 algorithm.
In brief, this approach uses the GO annotation of genes in a set of interest and compares them to those found in the entire background set, while also taking GO term hierarchy into account. The authors of topGO recommend against using adjusted p-values due to the high amount of nesting among GO terms (Alexa et al., 2006), so we employ a conservative significance threshold of p ≤ 0.01 for enrichment and visualization. However, when comparing affected biological processes between treatments, we also consider GO terms at p ≤ 0.05 to avoid possible downward biases in reporting only shared terms under a strict p ≤ 0.01 threshold. We used the orthogroup transcriptome as the background gene set and designated orthogroups with FDR < 0.05 as our gene set of interest for each comparison. We focused our analysis on biological processes GO terms because this ontology was the most relevant to our research questions. GO terms significantly enriched in at least one treatment for either species were included in a similarity matrix using GoSemSim (Yu et al., 2010) and used to perform hierarchical clustering for visualization of results (sensu Wellband & Heath, 2017). Informative grouping categories were determined by examining a range of clustering cutoffs to determine the threshold that included all of the major branches within which GO terms were related to a clear parent GO category; parent terms and descriptions for GO terms within each group were determined using AmiGO (Carbon et al., 2009).

| Thermal tolerance and shared orthology
Acute thermal tolerance was 33.9°C ± 0.15 for juvenile Inland Silverside and 28.2°C ± 0.09 for Delta Smelt (mean ± SEM) when acclimated to 16°C. Species was the only significant predictor of thermal tolerance (p < 0.0001), whereas fish morphometrics had no effect (Figure 1; see Table S1 for model estimates). This confirms our hypothesis that Inland Silversides have a higher thermal tolerance relative to Delta Smelt. A total of 11,829 orthogroups were identified with 11,782 having at least one transcript from both the Delta Smelt and Inland Silverside transcriptomes. An overall total of 41.8% of the peptide sequences translated from the transcriptomes were assigned to an orthogroup (50.2% and 37.2% for Delta Smelt and Inland Silversides, respectively), and the mean orthogroup size was 3.3 transcripts.

| Transcriptional responses
Differential expression analyses revealed support for higher thermal tolerance of Inland Silversides being likely facilitated by greater transcriptional responses relative to Delta Smelt, although thermal stress did invoke transcriptional changes in some common genes across species. In both species, the number of differentially expressed orthogroups increased with temperature ( Figure 2; Table S2); however, at analogous sublethal thermal stress treatments the number of orthogroups and overall magnitude of response was stronger in Inland Silversides relative to Delta Smelt (Figures 3 and 4). For Delta Smelt, very few orthogroups were differentially expressed at the lowest temperature (CT max -6°C) in the initial time treatment, and none were differentially expressed after 1 h of recovery. Inland Silversides also had relatively few orthogroups differentially expressed at the two lower temperatures (CT max -8°C and CT max -12°C) at the initial time point, but a considerable number after 1 h of recovery, particularly at CT max -8°C. We observed similar patterns across treatments in the number of thermally responsive genes private to each species (i.e., those without an identified ortholog; Table S3). Thus, although Delta Smelt clearly have the capacity for some modulation of gene expression in the face of sublethal thermal stress, at many of the same relative temperatures Inland Silversides displayed stronger up-and down-regulation of suites of orthogroups ( Figure 5).
The two highest sublethal stress thresholds (CT max -2°C and CT max -4°C) had a substantial overlap of orthogroups differentially expressed between the two species, with Inland Silversides generally displaying stronger log-fold changes (Figures 5 and 6). For example, at CT max -2°C both species differentially expressed orthologs of genes commonly involved in cellular stress response such as HSP70, SERPINH1, ULK2, DDIT4, and UBE2 (see Table S2 for full annotated results for each contrast). A number of these genes were also differentially expressed at CT max -4°C for both species, but not at lower relative temperature treatments. Interestingly, many differentially expressed orthogroups were also exclusive to CT max -2°C and CT max -4°C at both time points in Inland Silversides and CT max -2°C at the initial time point for Delta Smelt (Figure 6).
Though both species displayed transcriptional responses at common absolute temperatures (22°C and 26°C), these responses strongly differed between the species with regard to which orthogroups were differentially expressed (Figures 5 and 6) and may reflect different physiological processes (e.g., acclimation vs. cellular stress response). At 26°C, Delta Smelt had substantially more orthogroups differentially expressed at the initial time point, but after the 1 h of recovery this pattern was reversed. However, at 22°C as well as 26°C, few orthogroups were shared between the species at either time point (Figure 6; Table S2), despite the fact that both species were acclimated to and began thermal ramps at the same temperature (16°C). Exploratory contrasts between the handling controls of Delta Smelt and Inland Silversides also revealed a large number of differentially expressed genes between the species. F I G U R E 2 Reaction norms of differentially expressed orthogroups in response to temperature challenges in Delta Smelt (red) and Inland Silversides (blue) at two time points (triangle = initial, circle = recovery). Points are given a 0.5 position dodge for visualization of data

| Affected biological processes
Our comparison of sublethal responses between the two species using functional enrichment analyses provided evidence that greater detrimental physiological impacts likely occur at lower temperatures in Delta Smelt relative to Inland Silversides. However, importantly it also revealed more nuanced patterns of the specific biological processes affected by thermal stress within and across the species.
Functional enrichment analyses identified a total of 469 GO terms for biological processes that were enriched at p ≤ 0.01 in at least one treatment (Figure 7; Table S4). We determined a clustering threshold of 0.92 to be a representative cutoff for the dendrogram visualization, consisting of 21 groups of broader biological processes.
Concordant with transcriptional response results, the number of enriched GO terms increased with temperature for both species, and Inland Silversides generally had more enriched GO terms relative to Delta Smelt at the higher relative temperature treatments (Table S5).
Within each species, more GO terms were generally enriched at the two highest temperature treatments (CT max -4°C and CT max -2°C).
For Delta Smelt, at the lowest temperature treatment (CT max -6°C), many of the enriched terms were within groups reflective of cytokine production, amino acid and ion transport, and immune cell differentiation (Table S4). At CT max -4°C, enriched terms included metabolic and biosynthetic processes, protein folding/chaperone-mediated folding, and mitotic cell cycle regulation. All three of these groups were also represented in the CT max -2°C treatment, as well as DNA damage and cellular stress response. In Inland Silversides, at CT max -12°C enriched processes related to signaling, cell differentiation and development, fatty acid elongation, biosynthetic processes, and DNA repair. At CT max -8°C, particularly after recovery, cell differentiation and development and biosynthetic processes were enriched as in CT max -12°C, as well as a number of toll-like receptor pathways, histone acetylation, cell cycle regulation, and some cellular response to DNA damage and ubiquitination. At the two higher temperatures (CT max -4°C and CT max -2°C), an increasing proportion of enriched processes were related to stress response, DNA damage repair, cell cycle regulation and cell division, ubiquitination, mRNA transcription and processes, apoptosis, and oxygen homeostasis. As in CT max -8°C, toll receptor signaling and histone acetylation were consistently enriched at both temperatures, predominantly after recovery.
Although Delta Smelt and Inland Silversides displayed several to many enriched biological processes for treatments independently, F I G U R E 3 Log-fold change relative to the average log expression for each treatment and species; red and blue indicate orthogroups that were higher and lower (FDR < 0.05), respectively, relative to handling control treatments at each time point relatively few GO terms were shared between the species at common absolute or relative thermal treatments across either time points (Table S6). At 22°C, no terms were shared at the p ≤ 0.01 threshold and only cellular heat acclimation and copper ion transport were shared at the p ≤ 0.05 threshold. At 26°C, only one was shared at the p ≤ 0.01 threshold (circadian rhythm), but an additional 19 were shared at the p ≤ 0.05 threshold (Tables S4 and S6), including histone acetylation, responses to corticosterone and cytokines, and toll-like receptor signaling pathways. Similarly, in the relative treatments, no GO terms were shared at the p ≤ 0.01 threshold for the CT max -6/8°C and only four were shared at the p ≤ 0.05 threshold. Both CT max -4°C and CT max -2°C also had few shared terms between species at p ≤ 0.01, though this increased substantially at p ≤ 0.05 (22 and 51 shared GO terms, respectively, Table S6). Although TopGO takes into account GO hierarchy in Fisher's exact tests for enrichment, it is possible GO term specificity could limit overlap between treatments even when similar parent biological processes are affected.
Comparisons of hierarchical clustering groupings (Figure 7) represented by the GO terms enriched at p ≤ 0.01 does depict overlap between species for both absolute and relative temperature treatments (Tables S7 and S8). However, the degree of overlap displays a pattern largely similar to GO term comparisons, where greater number of these broader biological categories are shared at higher temperatures (e.g., at CT max -2°C 15 of the 21 groups were present in both species while at CT max -6/8°C only seven were shared).

| D ISCUSS I ON
Understanding the influence of thermal plasticity on organismal tolerance and interspecific competition is one key element of forecasting biological responses to climate change. There is evidence that heat tolerance directly facilitates introduced species in a warming climate (Bates et al., 2013) and in particular that acute thermal tolerance may play a critical role in extreme temperature events that result in altered community structure (Smale & Wernberg, 2013).
Under current and predicted future temperature regimes in the San Francisco Estuary, higher thermal tolerance of Inland Silversides likely contributes to the success of this abundant, widespread invasive species over the declining endangered native Delta Smelt. In combination with other traits that promote invasion success such as greater diet breadth and dispersal potential, these capabilities F I G U R E 4 Distributions of log-fold change in expression including all orthogroups that were differentially expressed (FDR < 0.05) in at least one treatment for Delta Smelt (left column) and Inland Silverside (right column). Initial and recovery data are illustrated as pairs of panels for each temperature, where top member is initial and bottom member is recovery may enable broader habitat and resource exploitation, ultimately providing competitive advantages that will further increase under projected warming scenarios in the San Francisco Estuary (Brown, Komoroske, et al., 2016). These findings may also be relevant for species interactions with other mesothermal fish in this system such as Longfin Smelt (Spirinchus thaleichthys), another native species of conservation concern that has been predicted to be at higher risk of climate warming impacts than Delta Smelt (Jeffries et al., 2016). Our study also provides comparative data for juvenile fishes, a critical life stage for which there is little information on thermal tolerance and stress response relative to adult stages (Bates et al., 2013). Due to life cycle timing, the critical growth phase for juveniles of many species occurs during summer months with exposure to the warmest conditions and extreme temperature events (Brown, Komoroske, et al., 2016;Nobriga et al., 2008), and may be an important link to species persistence under global climate change.
Capacity for transcriptional plasticity may be a key determinant of physiological resilience to environmental stress (López-Maury et al., 2008), and transcriptional profiling linked with phenotypic data offers important insight into evolutionary and ecological processes as well as conservation applications (Alvarez et al., 2015;Connon et al., 2018;Whitehead, 2012). Broadly, Inland Silversides responded to sublethal temperature stress with greater transcriptional changes relative to Delta Smelt, both in terms of the number F I G U R E 5 Log-fold change in expression including all orthogroups that were differentially expressed (FDR < 0.05) in at least one treatment for Delta Smelt (left) and Inland Silverside (right). Initial and recovery data are illustrated as pairs of columns for each temperature, where left member is initial and right member is recovery of differentially expressed genes and in the magnitude of the response (as measured by log-fold changes). Additionally, although there was some overlap in the identity of genes that responded to thermal stress between the two species, many genes were responsive only in Inland Silversides. These differential thermal responses support our hypothesis that if Inland Silversides displayed enhanced thermal tolerance relative to Delta Smelt, they would also show superior abilities to elicit transcriptional responses to cope with thermal stress. These results are also in concordance with previous work showing that enhanced transcriptional plasticity is associated with increased thermal tolerance (Garvin et al., 2015;Narum & Campbell, 2015;Sandoval-Castillo et al., 2020;Wellband & Heath, 2017), as well as limited overlap among the identities of genes involved in adaptive transcriptional plasticity across ecotypes (Sandoval-Castillo et al., 2020). However, other studies have found evidence that physiological resilience is correlated with large transcriptional responses that are transient in time (i.e., less resilient populations have large transcriptional responses that persist for longer, perhaps due to failure to attain a new level of homeostasis; Brennan et al., 2015;Franssen et al., 2011). Though we found that Inland Silversides had similar or greater differential expression at the second time point relative to the first for multiple treatments, we only examined transcriptional responses to acute thermal stress over a very short time period. Future work using comparative time-course experiments explicitly designed to test these ideas is needed to understand the temporal dynamics F I G U R E 6 Shared and exclusive differentially expressed (FDR < 0.05) orthogroups among treatments with five or greater orthogroups within the set. Initial and recovery data are illustrated as pairs of rows for each temperature, where the first member is initial and second member is recovery F I G U R E 7 Gene Ontology (GO) enrichment for biological functions including GO terms that were p ≤ 0.01 in the Fisher's exact tests for at least one treatment. Each row is a GO term, clustered in the dendrogram by similarity. Initial and recovery data are illustrated as pairs of columns for each temperature, where left member is initial and right member is recovery. Bar on the right indicates grouping for each section with a clustering threshold of 0.92 and labels based on descriptive categories encompassing GO terms within each group of transcriptional plasticity and how it translates to physiological resilience. It is also important to recognize that a strong transcriptional response alone may not indicate greater tolerance, and conversely, smaller magnitude responses can have large phenotypic effects and fitness consequences (Evans, 2015;Oomen & Hutchings, 2017). Additionally, non-or maladaptive plastic responses can occur in new environments at multiple biological levels (Ghalambor et al., 2007(Ghalambor et al., , 2015, so interpreting how transcriptional responses relate to physiological resilience without phenotypic trait data is challenging. Coupling complementary molecular and organismal metrics as we have done in this study is critical for distinguishing compensatory responses from those that indicate irreversibly detrimental impacts on performance and fitness (Connon et al., 2018). While only a small proportion of studies to date have made such connections across biological levels (Alvarez et al., 2015), this approach is becoming more common and is essential for advancing our understanding of how alterations at the molecular level affect individual performance, survival, and reproductive success. The results of our study linking thermal tolerance and transcriptional plasticity support the hypothesis that successful invaders may have a particular set of mechanisms conferring heat tolerance that differ from their native counterparts (Bates et al., 2013). Since first proposed, this idea has largely remained untested, and further studies are needed to determine whether this pattern can be broadly generalized, as we only compared one invasive and native species. However, our results also align with previous work demonstrating the role of enhanced thermal tolerance in favoring invasive species under climate warming (Bates et al., 2013;Dukes & Mooney, 1999).
The majority of research to date investigating the role of transcriptional plasticity in physiological resilience has focused on a single temperature and time point in a small number of individuals, but it is well-established that physiological responses are dependent on both the magnitude and duration of the environmental stressor and can vary across individuals, particularly in wild populations that may have mixed individual genotypes (reviewed in Connon et al., 2018;Oomen & Hutchings, 2017). In line with our second objective to assess whether the endangered native displayed signatures of a stress response at lower sublethal thresholds, we quantified global transcriptional responses across a gradient of thermal conditions, enabling comparisons at multiple thresholds that provide a more nuanced perspective on how thermal responses are conserved or diverge between these two species with different evolutionary histories. Detailed comparisons at mutual relative thresholds revealed that while Delta Smelt initiated most transcriptional responses at 2-4°C below their acute thermal limits, Inland Silversides also initiated responses at 8°C below their acute thermal limit, supporting a prior study suggesting that mesothermal Delta Smelt differ in their abilities to mount transcriptional responses and maintain homeostasis under sublethal thermal challenges relative to eurythermal fishes (Komoroske et al., 2015). At the highest relative thermal challenge 2°C below their respective upper tolerance limits, both species displayed strong transcriptional changes in line with cellular stress responses that are well known to be widely conserved (Kassahn et al., 2009). However, the heightened magnitude of change in Inland Silversides for many of the differentially expressed genes common to both species suggests that they may have greater abilities to mount compensatory physiological responses under thermal stress.
In addition to the observed conserved response, the fact that both species, but particularly Inland Silversides, also exhibited transcriptional changes in a large number of genes that had orthologs in Delta Smelt but were not thermally responsive, suggest that the species may also employ different molecular pathways to cope with thermal stress. Genes that show transcriptional plasticity in response to different environments have been found to be often linked to phenotypic variation and involved in adaptive divergence between populations (reviewed in Gibert, 2017), so these genes responsive to thermal challenges in Inland Silversides may be of high interest for future studies of rapid thermal evolution. The hypothesis that Inland Silversides may employ different molecular pathways to cope with thermal stress is further supported by the higher number of differentially expressed private genes (i.e., those that did not have an ortholog identified in Delta Smelt) in Inland Silversides relative to Delta Smelt. However, further study of these potential species-specific genes with high-quality reference genomes is needed to fully understand their origins, molecular functions, and roles in influencing thermal tolerance. As genomic resources are rapidly improving and becoming available for larger-scale comparative analyses, further examination of these questions should deepen our understanding of the diverse genomic and physiological "tools" that different species possess (or not) to cope with global change. Further, as plasticity is increasingly incorporated into population viability models and compared between greater numbers of species showing population declines versus resilience (Hendry, 2016), the inclusion of transcriptional plasticity and molecular mechanisms underlying organismal phenotypes will be useful to assess the likely prevalence of plastic rescue under future global change (Kovach-Orr & Fussmann, 2013).
The comparisons discussed above focus on similar relative sublethal thresholds, which occur inherently at different temperatures when two species differ in their thermal tolerance. Additionally, it is important to consider what drives differences in the absolute temperatures for when homeostasis is disrupted, and physiological responses are required. This is particularly relevant under different climate change scenarios, as recent work has shown that plastic capacity can be a target of climatic selection that may confer differential adaptive resilience to warming (Sandoval-Castillo et al., 2020). For example, why is it that acute exposure to 26°C, a 10 degree increase from acclimation conditions, is highly stressful for Delta Smelt but not for Inland Silversides? As described above, it is not until 30 or 32°C that Inland Silversides differentially express many of the cellular stress response genes that Delta Smelt up-or down-regulate at 24 and 26°C, so what drives this expanded range of thermal resilience of Inland Silversides? In the 26°C treatments, as well as at 22°C, we observed roughly similar numbers of genes being differentially expressed in both species, but the identity and logfold change greatly differed. These differences may reflect different physiological processes (e.g., capacity for rapid acclimation in Inland Silversides versus cellular stress response in Delta Smelt) and/or the superior ability of Inland Silversides to sustain homeostasis via rapid and/or large magnitude transcriptional changes. It is also plausible that Inland Silversides are able to cope with thermal challenges at this level through means aside from transcriptional plasticity.
Alternate thermal adaptations, such as more thermally resilient proteins through amino acid substitution, post-translational modifications, or higher constitutive expression of key proteins like molecular chaperones, can also be important in determining tolerance limits (DeBiasse & Kelly, 2016;Somero et al., 2017). Altered constitutive gene expression relative to physiological tolerance has been documented among locally adapted populations (Gleason & Burton, 2015;Maynard et al., 2018) and closely related species (Wellband & Heath, 2017). However, continued transcription and protein production can be energetically costly (Feder & Hofmann, 1999)  despite the challenges of working with sensitive species that often present limitations for data interpretation that are important to acknowledge, much insight can still be gained from these studies and their exclusion from study may bias our understanding of the molecular processes underlying the vast phenotypic diversity observed in nature (Komoroske et al., 2015).
Both species displayed greater enrichment of affected biological processes with increased temperature. Although the patterns in our data were more nuanced than the clear physiological shifts that have been previously described in other species (Kassahn et al., 2009;Logan & Somero, 2011), biological functions within categories classically associated with thermal stress such as chaperone-mediated protein folding and DNA and cellular damage detection (Logan & Buckley, 2015) were present in both species at the highest temperatures. These patterns reinforce the findings discussed above that greater detrimental physiological impacts occur at lower absolute temperatures in Delta Smelt. However, some processes related to cellular response to DNA damage and ubiquitination did occur at 8°C below the acute thermal limit for Inland Silversides, so it is possible that the eurythermal species does still experience some level of stress at this lower relative threshold that needs to be compensated for. Nonetheless, the prevalence of enrichment of toll-like receptor pathways and histone acetylation particularly in almost all treatments for Inland Silversides (while only enriched at higher relative thresholds and to a lesser degree in Delta Smelt) suggests these processes might play key roles in their expanded thermal resilience.
The low proportions of enriched GO terms shared between species at common absolute temperatures was expected and in agreement with previous work comparing molecular mechanisms at common temperatures in species with different thermal tolerances (Wellband & Heath, 2017). However, it was surprising that there were also relatively few shared enriched biological processes at the lower relative thermal thresholds. It is possible this could be a consequence of the smaller transcriptional responses in Delta Smelt in these treatments. Despite this, overall, the increase in shared GO terms and broader functional categories between species at the highest relative thermal challenges is concordant with the idea that transcriptional responses are reflective of biological responses at thresholds relative to physiological limits rather than absolute temperature, as supported by data in this study and others (Komoroske et al., 2015;Logan & Somero, 2011). Thus, when species differing in physiological tolerance are exposed to one common temperature, responses are inherently at different sublethal relative thresholds and may provide only a partial picture of the biological processes involved in the organisms' response to thermal stress (Jeffries et al., 2016).
Finally, functional annotation still remains one of the largest challenges in applying RNA-sequencing approaches to nonmodel organisms and ecological interpretations (Oomen & Hutchings, 2017). As a consequence, many studies, including ours, principally focus on the largest or most well-defined categories thought to be biologically conserved across taxa. However, this may omit important components of the multifaceted molecular networks involved in environmental stress responses or incorrectly attribute biological processes to genes where the functions may be different in the species of study. As functional annotation databases improve, more in-depth, detailed analyses and interpretations may be necessary to refine our understanding of how complex, coordinated transcriptional changes manifest to differential physiological performance, survival, and species interactions.
Identifying the mechanisms that underpin physiological resilience is a key goal to fundamentally understand how organisms cope with environmental challenges and determine the demographic and ecological consequences of climate change. Transcriptional plasticity may be one factor that can play an important function in facilitating heightened thermal tolerance and, in turn, population persistence, species interactions, and shaping community assemblages under warming, particularly in variable environments. In addition to providing important information on the sensitivity of a critically endangered native fish to climate warming, the findings of our study in these two species can inform future studies encompassing a wider range of species with varied phylogenomic divergence to more broadly quantify the roles of plasticity at multiple biological levels in shaping future ecosystems.

ACK N OWLED G EM ENTS
We would like to thank staff at the University of California, Davis K. Wellband for sharing R code for GO enrichment analyses. We are grateful to two anonymous reviewers and the editor for constructive feedback on a previous draft that improved the final version of this manuscript.

CO N FLI C T O F I NTE R E S T S
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Illumina raw reads and sample metadata are deposited in NCBI Sequence Read Archive (Bioproject PRJNA512860). Code for all analyses can be accessed at https://github.com/lkomo ro/Nativ e-Invas ive_trans cript ional_plast icity.