Skip to main content
Download PDF

Transcriptomic analysis provides insights into molecular mechanisms of thermal physiology

BMC Genomics volume 23, Article number: 421 (2022) Cite this article

Abstract

Physiological trait variation underlies health, responses to global climate change, and ecological performance. Yet, most physiological traits are complex, and we have little understanding of the genes and genomic architectures that define their variation. To provide insight into the genetic architecture of physiological processes, we related physiological traits to heart and brain mRNA expression using a weighted gene co-expression network analysis. mRNA expression was used to explain variation in six physiological traits (whole animal metabolism (WAM), critical thermal maximum (CTmax), and four substrate specific cardiac metabolic rates (CaM)) under 12 °C and 28 °C acclimation conditions. Notably, the physiological trait variations among the three geographically close (within 15 km) and genetically similar F. heteroclitus populations are similar to those found among 77 aquatic species spanning 15–20° of latitude (~ 2,000 km). These large physiological trait variations among genetically similar individuals provide a powerful approach to determine the relationship between mRNA expression and heritable fitness related traits unconfounded by interspecific differences. Expression patterns explained up to 82% of metabolic trait variation and were enriched for multiple signaling pathways known to impact metabolic and thermal tolerance (e.g., AMPK, PPAR, mTOR, FoxO, and MAPK) but also contained several unexpected pathways (e.g., apoptosis, cellular senescence), suggesting that physiological trait variation is affected by many diverse genes.

Peer Review reports

Introduction

Physiological traits define how species live, the habitats they can exploit, health, and energy allocation between ecological interactions or reproduction [1,2,3,4,5,6,7,8,9,10]. To understand the molecular mechanisms and evolutionary forces altering physiological traits, it is critical to understand the genes involved. Quantifying mRNA expression provides information for both heritable and plastic responses [11,12,13,14,15,16,17], offering insight where genome wide association studies (GWAS) or similar approaches may be less informative due to context dependent genetic architecture [19, 20]. Studying mRNA expression is also likely to provide a greater mechanistic understanding of complex traits in comparison to genome wide nucleotide variation because mRNAs are more often defined genes associated with biochemical or physiological pathways (e.g., through Kyoto Encyclopedia of Genes and Genomes [KEGG] or Gene Ontology [GO] terms or molecular investigation). Overall, an improved mechanistic understanding of traits enables us to parse redundancy across biological organization levels and reveal the evolutionary processes and genetic architectures contributing to phenotypic variation.

Here, we used mRNA expression patterns to identify molecular mechanisms underlying physiological traits in the small saltmarsh teleost fish, Fundulus heteroclitus. Thermal physiology is often studied in F. heteroclitus because they are highly plastic and have adapted to live in the highly variable temperate coastal salt marshes where temperature fluctuates 15 °C on daily and seasonal timescales [18, 21,22,23,24,25,26,27,28,29]. Yet, few studies have examined the molecular and genetic basis of physiological trait variation related to thermal responses in this species beyond specific gene expression (e.g., heat shock protein expression, [26, 30]), limiting our understanding of physiological response to temperature, which is likely to include 10 s or 100 s of expressed genes and may differ in response to variable acclimation conditions. To better understand thermal response molecular mechanisms in F. heteroclitus, we used individuals captured from three wild populations less than 15 km apart with little overall genetic distance among them [31,32,33] to examine physiological trait variation and then related physiological trait variation to heart and brain mRNA expression under 12 °C and 28 °C acclimation conditions (Fig. 1, [18, 34]). Traits included whole animal metabolic rate (WAM), critical thermal maximum (CTmax), and four substrate specific cardiac metabolic rates (CaM substrates: glucose [Glu], fatty acids [FA], lactate + ketones + ethanol [LKA], and endogenous [END]) [18]. Notably, there were few differences among the three populations in the six physiological traits measured. This allowed us to address two hypotheses: 1) genetically similar populations with little divergence in physiological traits will have a shared mRNA expression response to thermal acclimation, and 2) physiological traits will be related to temperature and tissue specific mRNA expression, with variance among traits in the mRNAs that explain trait variation.

Fig. 1
figure 1

Metabolic and thermal tolerance trait variation among individual Fundulus heteroclitus. Relative physiological trait performance at 12 °C (blue) and 28 °C (red) acclimation and assay temperatures for A)Whole animal metabolic rate versus thermal plasticity between acclimation temperatures, B critical thermal maximum versus thermal plasticity between acclimation temperatures, and C-F cardiac metabolism under 12 °C and 28 °C acclimation and assay temperatures with substrates C glucose, D fatty acids, E lactate ketones and ethanol (LKA), and F endogenous (no substrate). Previously, difference between individuals when acclimated and assayed at 12 °C (blue) and 28 °C (red) revealing up to 14-fold variance in mass corrected whole animal metabolism (for the most plastic individual), 0.25-fold variance in critical thermal maximum, and 13-fold variance in fatty acid cardiac metabolism (the most variable substrate). Adapted from Drown et. al. 2021. [18]

Full size image

Using mRNA expression in combination with physiological trait variation we found co-variation between physiological traits and heart and brain mRNA expression among 86 individuals across the three populations. Among individuals, we previously found mass corrected physiological traits to be highly variable with up to 14-fold difference (range spread) in metabolic rate (mass corrected), 0.25-fold difference in CTmax (mass corrected), and 13-fold difference in fatty-acid CaM (heart mass corrected, the most variable substrate) [18] (Fig. 1, [18, 34]). This degree of trait variance is comparable to that found in metabolic rates [35] among 77 fish species spanning thousands of kilometers (15–20° latitude, ~ 2000 km (Fig. S1, data from [35]), and for thermal tolerance among polar and temperate crabs [36]. This large physiological trait variation among individuals within a species provides a powerful approach to determine the relationship between mRNA expression and these heritable fitness related traits [16, 37,38,39,40,41]. We show that among these individuals, acclimation induced differential expression of 362 heart mRNAs and 528 brain mRNAs across all three populations, yet few differentially expressed mRNAs (one or less) were shared across all three populations. Within each acclimation temperature, co-expressed mRNA modules were significantly associated with WAM, CTmax, and CaM. Using KEGG and GO enrichment, we identify biologically relevant networks among co-expressed mRNA modules that explain these traits. These data link a simpler molecular phenotype (mRNA expression) to complex trait variation to enhance our understanding of biological pathways that affect these traits and may be important for evolutionary adaptation.

Results

Differential expression analysis

Sequencing statistics and sample sizes are summarized in Table 1. Differential expression patterns among populations and acclimation temperatures were identified using DESeq2 [42]. First, to examine population and temperature specific expression we used model design: (~ Population + Acclimation-Temperature + Population*Acclimation-Temperature). This analysis revealed significant Population*Acclimation-Temperature interactions, suggesting acclimation temperature specific mRNA expression patterns among populations. Because of the significant interactions, we analyzed individuals acclimated to 12 °C or 28 °C separately with model design: ~ Population to identify differentially expressed mRNAs among populations within an acclimation temperature. Similar to other species, there were significant differentially expressed mRNAs between acclimation temperatures, reflecting changes in response to environmental temperature ([14, 43, 44], Table S1). Across all 3 populations, hearts had 362 mRNAs (3.5% of total) that were significantly different between the two acclimation temperatures (FDR < 0.05) with an equal number of up and down regulated mRNAs for 12 °C versus 28 °C (180 up and 182 down). For brains, 528 mRNAs (4.8% of total) were significantly differently expressed between the two acclimation temperatures (FDR < 0.05) with ~ 2.5-fold more down regulated at 28 °C relative to 12 °C (148 up and 380 down). While all three populations showed acclimation effects for heart and brain mRNAs, the affected mRNAs were not shared among all populations (Fig. 2, Table S1). Additionally, acclimation significant mRNAs were unique to either heart or brain with no shared (0 mRNAs) acclimation response between tissues. This reflected different expression patterns between tissues, previously identified with PCA analysis (Fig. S2).

Table 1 Sequencing statistics. Sequencing statistics and sample size distribution among tissues, acclimation temperatures, and populations
Full size table
Fig. 2
figure 2

Population and tissue specific transcriptomic response to acclimation temperature. Number of differentially expressed mRNAs (DEGs) within each population between 12 °C and 28 °C acclimated hearts (A and C) and 12 °C and 28 °C acclimated brains (B and D). For both heart and brain, populations had many unique DEGs (C and D, upregulated = black, downregulated = grey) that were differentially expressed between acclimation temperatures and few shared DEGs (C and D, shared = white), with only 1 DEG shared among all three populations for brain (LOC118561484 in brains). Population and number of DEGs are not independent, Chi-Squared test, heart p = 1.64 × 10–10, brain p = 2 × 10–16

Full size image

At each acclimation temperature, populations also had significant expression differences (Fig. 3, Table S2). Hearts at 12 °C and 28 °C had 158 or 153 differentially expressed mRNAs among populations, respectively (Table S2). These represent 1.50% or 1.45% of all expressed heart mRNAs at 12 °C and 28 °C, respectively; brains had 242 or 330 differentially expressed mRNAs among populations at 12 °C and 28 °C, respectively. These represent 2.21% or 3.02% of all expressed brain mRNAs at 12 °C and 28 °C, respectively. None of the population effects were significant across all three populations (Fig. S4) for any acclimation temperature or tissue.

Fig. 3
figure 3

Temperature specific differential mRNA expression among populations. A Differential expression for 158 heart mRNAs and B 242 brain mRNAs between any two populations at 12 °C. In hearts, 43.0% of mRNAs (68/158) are shared among any two population comparisons. In brains, 32.2% of mRNAs (78/242) are shared among any two population comparisons

Full size image

Importantly, there is an adaptive hypothesis: that the mRNAs in the anthropogenically warmed population, TE, are uniquely different, where TE is significantly different from both northern and southern reference populations with no significant differences between the references [14, 45]. For heart mRNAs at 12 °C there are 10 mRNAs (6.33% of significant mRNAs), and at 28 °C there are 3 mRNAs (1.96%) where the TE population is uniquely different from both references (Fig S4, Table S2). For brain mRNAs at 12 °C there are 11 mRNAs (4.55% of all significant mRNAs), and at 28 °C there are 27 mRNAs (8.18%) where the TE population is uniquely different from both reference populations. While the overall frequency of differentially expressed genes is small (1.45% to 3.02% vs. total 10 k mRNAs), the pattern where the TE population is different from both northern and southern reference populations but the two reference populations are not different aligns with an adaptative hypothesis.

Variation in mRNA expression

In addition to differential expression analysis, we were interested in the degree of mRNA expression variance. Previously, we found that variation in WAM, CTmax, and substrate specific CaM was greater at 12 °C than at 28 °C. Here we found that both heart and brain tissue had greater average coefficient of variation (CV, standard deviation/mean*100) across all mRNAs at 12 °C than at 28 °C (T-test, heart p = 9.587e-05, brain p = 0.02014), similar to our findings of greater variation in physiological traits measured at 12 °C.

Weighted gene co-expression network analysis

We used weighted gene co-expression network analysis (WGCNA, [46]) to detect co-expressed mRNA clusters. WGCNA approaches group mRNAs with similar expression patterns into independent modules. Expression patterns for all mRNAs within a module were then summarized into principal components called module eigengenes (MEs, Table 2 for heart mRNAs and Table 4 for brain mRNAs), and these MEs were correlated to each of the six physiological traits (Table 3 shows significant heart ME trait correlations, and Table 5 shows significant brain ME trait correlations). Each ME has a “hub-MM”, the mRNA with the highest correlation to the ME, with MM being the correlation coefficient (Tables 2 and 4).

Table 2 Heart Significant Modules
Full size table
Table 3 Heart Significant Module Trait Correlations
Full size table
Table 4 Brain Significant Modules
Full size table

MEs were correlated to the body mass residuals for the six traits (WAM, CTmax, and the four substrate specific CaM). These analyses were done across all three populations because populations did not have any significant differences among traits. Each of the ME-trait correlations had a “hub-GS” – the mRNA in the module with the highest correlation to the trait, with GS (gene specific) being the correlation coefficient for this single mRNA. Both heart and brain WGCNA analysis used a minimum module size of 30 mRNAs per module and combined modules with correlation > 75% (see methods). To verify that trait versus ME correlations were not driven by spurious outliers, we used a jack-knife approach to subsample 90% of individuals and repeat ME-trait correlations 100 times. Correlations that were significant in at least 70 out of 100 repeated correlations in the same direction (positive or negative correlation coefficient) were retained for further analysis (see methods).

Heart WGCNA

For heart mRNAs, we found 39 co-expression modules with 90 to 554 mRNAs in each module (Table 2), and these heart MEs (first principal component of co-expressed mRNAs) were correlated to six physiological traits at each acclimation temperature. There were 12 significant ME-trait correlations: 9 heart MEs with 5 temperature specific traits (FDR < 0.05, Table 3, Figs. 3 and 4). Traits correlated with at least one of these 9 heart MEs included: at 12 °C WAM, CTmax, FA CaM, heart mass, LKA CaM, and at 28 °C, CTmax (Table 3, Figs. 4, 5). Two of these modules (ME4_heart, ME5_heart) were correlated with both WAM at 12 °C and FA CaM at 12 °C, and one module (ME6_heart) was correlated with both LKA CaM at 12 °C and heart mass at 12 °C. WAM at 12 °C had the most significant ME correlations (4 total), followed by FA CaM at 12 °C (3 total) and LKA CaM at 12 °C (2 total). The other three traits were each correlated with a single module (Table 3). On average, a single heart module explained 55% of variance for one trait with a minimum of 48.5% (ME1_heart with CTmax 12 °C) and a maximum of 65% (ME6_heart with LKA 12 °C).

Fig. 4
figure 4

Significant cardiac metabolism-heart module correlations from weighted gene co-expression network analysis. Significant correlation of fatty acid cardiac metabolic rate at 12 °C (N = 16) with ME3_heart (A), ME4_heart (B), and ME5_heart (C). Significant correlation of lactate, ketone, and alcohol (LKA) cardiac metabolic rate at 12 °C (N = 19) with ME6_heart (D) and ME7_heart (E). Pearson correlation coefficients (Cor) and FDR p-values are displayed for each significant correlation

Full size image
Fig. 5
figure 5

Significant whole animal trait-heart module correlations from weighted gene co-expression network analysis. Significant correlation of critical thermal maximum at 12 °C (N = 17) with ME1_heart (A), critical thermal maximum at 28 °C (N = 19) with ME2_heart (B), whole animal metabolic rate at 12 °C (N = 16) with ME8_heart (C), ME9_heart (D) ME4_heart (E), and ME5_heart (F). Pearson correlation coefficients (Cor) and FDR p-values are displayed for each significant correlation

Full size image

For traits that were significantly correlated with more than one ME, a multiple correlation coefficient was calculated. For WAM at 12 °C, the four significant MEs together had a multiple correlation coefficient of 82%, the three significant MEs for FA CaM at 12 °C had a multiple correlation coefficient of 79.5%, and the two significant MEs for LKA CaM at 12 °C had a multiple correlation coefficient of 75.5%. All modules correlated with FA CaM at 12 °C and CTmax at 12 °C had positive correlation coefficients while all other significant trait versus ME correlations in hearts had negative correlation coefficients.

Brain WGCNA

For brain mRNAs, we found 42 co-expressed modules with 142 to 393 mRNAs per module (Table 4). There were 6 significant ME-trait correlations (FDR < 0.05) that included 4 unique modules and 4 temperature specific traits (Table 5, Fig. 6): at 12 °C, body mass and at 28 °C, CTmax, WAM, and body mass. The trait with the most significant correlations was CTmax at 28 °C (3 significant ME’s), and two of these were also significant with body mass at 28 °C (ME3_brain) or WAM at 28 °C (ME4_brain). On average, the correlation coefficient for a brain ME was 62% with a minimum of 56.8% (ME3_brain with CTmax 28 °C) and a maximum of 70.2% (ME4_brain with WAM 28 °C). For CTmax at 28 °C, which was correlated with three MEs, the multiple correlation coefficient was 71.7%. All correlations between traits and brain MEs were negative except for body mass at 28 °C.

Table 5 Brain Significant Module Trait Correlations
Full size table
Fig. 6
figure 6

Significant whole animal trait-brain module correlations from weighted gene co-expression network analysis. Significant correlation of critical thermal maximum at 28 °C (N = 17) with ME3_brain (A), ME4_brain (B), and ME2_brain (C). Significant correlation of whole animal metabolic rate at 28 °C (N = 16) with ME4_brain (D). Pearson correlation coefficients (Cor) and FDR p-values are displayed for each significant correlation

Full size image

KEGG and GO enrichment

Critical thermal maximum enriched terms

MEs were tested for KEGG and GO term enrichment using the complete set of tissue specific mRNAs as the gene universal or reference set (10,535 for heart, 10,932 for brain; Tables S3 and S4). For CTmax, all 5 MEs were significantly enriched for KEGG pathways (ME1_heart, ME2_heart, ME2_brain, ME3_brain, and ME4_brain) and included MAPK signaling, mTOR signaling, glyoxylate and dicarboxylate metabolism, insulin signaling pathway, glutathione metabolism, metabolic pathways, carbon metabolism, and tryptophan metabolism. Enriched GO terms included regulation of organ growth and cellular stress response in heart and AMP-activated protein kinase (AMPK) activity in brain. ME for CTmax at 28 °C contained substantial overlap in enriched KEGG pathways related to metabolism with 5 out of 8 terms enriched in both heart and brain modules. So, although there were different mRNAs in heart and brain modules correlated with CTmax at 28 °C, the KEGG terms related to metabolism were shared (5 out of 8), suggesting that different mRNAs in heart and brain belonged to similar pathways impacting CTmax at 28 °C.

Metabolic rate enriched terms

Modules significantly correlated with WAM were enriched for KEGG pathways including oxidative phosphorylation (heart only), glutathione metabolism (brain only), and metabolic pathways (both heart and brain). In addition, ME3_heart was correlated with FA CaM at 12 °C and enriched for KEGG terms including metabolic pathways, forkhead protein (FoxO) signaling, and metabolism of NADH derivatives (nicotinate and nicotinamide). One module, ME5_heart, was significantly correlated with WAM at 12 °C and FA CaM at 12 °C and contained several KEGG pathways directly related to fatty acid metabolism as well as known transcription factors like PPAR that impact metabolic homeostasis by controlling expression of many metabolism related genes [47]. Previously, partial correlation coefficients between FA CaM at 12 °C and WAM at 12 °C were negatively correlated [18], and similarly ME_5 had opposite correlation coefficients for these two traits (Table 3, Figs. 4, 5). This correlation between traits and their correlations to MEs also occurs for CTmax and WAM at 28 °C [18] and ME4_brain (enriched for glutathione metabolism and metabolic pathways). The fact that the same MEs are associated with traits that have significant partial correlations emphasizes the biological relevance of the MEs.

Discussion

Understanding the genes that affect physiological performance is one of the more important research goals for human health [48,49,50] and evolutionary ecology, especially with global climate change [2, 6, 10, 51,52,53]. Here, to provide insight into the genes that affect physiological performance, we relate the mRNA variation with the variation in six physiological traits (Fig. 1, [18, 34]). This approach is strengthened by the large variation among F. heteroclitus individuals used in this study: the variation (range spread) is 14-fold for WAM, 13-fold for FA CaM, and 0.25-fold for CTmax [18] (Fig. 1). These physiological trait variations are larger than found among different species spanning temperate to tropical waters (Fig. S1). Additionally, five of six traits are not significantly different among the three populations (F-values < 2), with most trait variation found among individuals and not among populations [18]. This substantial interindividual variation in traits allowed us to examine interesting physiological patterns related to acclimation response and temperature performance. For example, individuals with low WAM at 12 °C tended to have high WAM at 28 °C and visa-versa; the substrates that supported CaM differed among individuals; and, acclimation response eliminated temperature effects such that Q10 (change in rate for every 10 °C increase in temperature) was approximately 1.0 for cardiac metabolic rates (nearly equal rate when measured at 12 °C and 28 °C, [18]). The large physiological variation among individuals, the difference in substrate use, and the strong acclimation responses provide a unique resource to begin to understand which mRNAs are associated with physiological performance and thus the genes responsible for health and response to environmental change.

The role of mRNA expression in acclimation and evolution

For mRNA expression in both heart and brain tissues, we found significant interactions between acclimation and population effects: the expression of several hundred mRNAs differed between acclimation temperatures, but these were not shared among all three populations (Fig. 2). Previously (Fig. 1., [18]), in these same individuals, acclimation response to 12 °C and 28 °C affected all six physiological traits (WAM, CTmax, and the four substrate specific CaM). For CTmax, there was an expected enhancement: higher CTmax in individuals experiencing warmer environments. For metabolic rates (WAM and CaM), acclimation to 12 °C and 28 °C mitigated the effect of temperature [18]. Specifically, without physiological acclimation there is an expected ~ threefold increase in metabolic rates with the 16 °C increase in acclimation and assay temperature (i.e., with a doubling for every 10 °C) [54]. Yet, WAM had only ~ 1.2-fold increase [43] from 12 °C to 28 °C, and CaM had no significant increase between temperatures [18]. Presented here, across all three populations, acclimation produced significant differential mRNA expression (FDR < 0.05) in hearts (362 mRNAs) and brains (528 mRNAs, (Table S1)). These mRNA expression changes associated with acclimation responses are similar to prior studies among ectotherms where transcriptomic response to temperature acclimation enhances thermal performance [43, 44, 55, 56]. For example, in three-spine stickleback (Gasterosteus aculeatus) and other fishes, metabolic enzyme expression and mitochondrial volume density increase in response to cold acclimation can compensate for reduced enzyme catalytic rate with decreased temperature [43, 57, 58]. Similarly, in eastern oysters (Crassostrea virginica) among a suite of environmental factors (temperature, pH, salinity, dissolved oxygen, etc.), temperature was the most important transcriptomic variation predictor with thermal stress increasing oxidative phosphorylation transcript expression [59]. Even Antarctic fish, which are adapted to extreme and invariable cold, show plasticity in metabolic transcripts with temperature acclimation that impacts whole animal performance [58, 60].

While quantitative gene expression changes are common with acclimation (affecting both mRNA and proteins [55, 61]), surprisingly, mRNA acclimation responses were different among populations—for heart mRNAs, no significant acclimation responses were shared among all three populations, and for brains only one mRNA was shared among the three populations. Further, 88–98% of significant acclimation responsive mRNAs are unique to each population. In contrast, the six physiological traits’ acclimation responses were not different among populations. All populations were subjected to a common environment for a long time (~ 1 year or nearly 30–50% of a fish’s expected life span) with acclimation to the 12 °C and 28C. Thus, the difference in acclimation mRNA response among populations was not due to short-term physiological effects and may be due to genetic polymorphisms driving acclimation responses but could also be due to irreversible developmental effects or trans-generational effects. Regardless of the genetic mechanisms responsible for the divergent mRNA acclimation responses among populations, these data suggest that multiple different mRNA expression patterns drive acclimation responses. This conclusion is similar to CaM measurements in Maine and Georgia populations, where the mRNAs that explain substrate specific metabolism varied among groups of individuals [16]. The observations, that plasticity in the six physiological traits between temperatures is similar among populations yet mRNA acclimation responses differ among populations, suggest that multiple redundant molecular mechanisms drive temperature compensation.

A single difference occurs among populations for the six physiological traits: endogenous CaM at 28 °C. Yet, populations had significant mRNA expression differences specific for each temperature, and none of the population significant mRNAs were shared at 12 °C or 28 °C (Fig. S4, Table S2). One pattern indicative of adaptation occurs where the anthropogenically warmed TE population is significantly different from both northern and southern reference populations (not heated by thermal effluent from nuclear power plant) but not different between these two references [14, 31, 45]. While we did not find this adaptive pattern for any of the physiological traits, a few mRNAs (2–8%) did show this pattern. Yet, this pattern is not unlikely given the number of differences among all populations. That is, it is probable that 2%-8% of mRNAs for either tissue or at either acclimation temperature would have this potentially adaptive pattern of mRNA expression. Thus, while some mRNA expression patterns are exclusive for the anthropogenically heated TE population relative to both references, an adaptive hypothesis is possible but not strongly supported by our mRNA data. Regardless, the differences here were for mRNA expression that can be affected by DNA polymorphisms or influenced by environment (i.e., GxE). Thus, the difference among populations in mRNA expression are dependent on the thermal environment, and, if adaptive, suggest that the different genetic polymorphisms are responsible for adaptive divergence at different temperatures. This conclusion is similar to a comparison within and among species across a wider geographic range where adaptive divergence in mRNA expression was dependent on the thermal environment [14]. For the TE population, Dayan et al. [31], concluded that there was adaptive divergence based on evolutionarily significant DNA polymorphisms. We would extend this to suggest that populations have evolved different mRNA expression patterns that are dependent on the thermal environment but that, nevertheless, produce similar physiological phenotypes. This is consistent with a redundant polygenic trait architecture, which may be expected for some complex traits that are affected by many physiologically processes [48]. While mRNA expression variation is largely attributed to DNA polymorphisms [62, 63], other possible mechanisms including epigenetic modification (e.g., DNA methylation) are also known to impact mRNA expression and may impact the physiological trait variation we have measured [64,65,66].

Biological relevance of co-expressed mRNAs

This study’s focus is to identify differential mRNA expression responsible for the large physiological trait variation previously characterized [18]. We examine co-variation among mRNAs and relate this to the six physiological traits using a weighted gene co-expression network analysis (WGCNA) [46]. WGCNA identified co-expressed mRNA modules, MEs, highly correlated with WAM, CTmax, FA CaM, LKA CaM, and body and heart mass, depending on the acclimation temperature (Figs, 4, 5, and 6). The average ME-trait correlation was 0.55 for heart and 0.62 for brains (Tables 3 and 5). These MEs, containing 90–554 mRNAs each, contained few (0–10) mRNAs with significant expression differences among populations, suggesting that the variation among individuals, and not differences among populations, is most important. This is in agreement with five of the six physiological traits, where most the variation is within and not between populations (exception of CaM END at 12 °C). WGCNA has been previously used to identify mRNA expression networks important for various pathologies including cardiovascular disease [67, 68], cancers [69,70,71,72,73], and diabetes [68, 74], among others. In non-human organisms, WGCNA has been used to characterize response to the environment, including heat stress in turbot [75], carotenoid metabolism in apricot fruit [76], and disease response in corals [77]. Although few studies, to our knowledge, have validated correlations using jack-knife subsampling to ensure that the correlations were consistent among most individuals and not driven by a few outliers as we did, these studies similarly identified potentially meaningful correlations between traits and co-expressed mRNAs.

Our jack-knife subsampling indicates that the high correlations between traits and co-expressed mRNAs are found among most individuals and not driven by a few outliers. More importantly, in this study, the correlation patterns between MEs and physiological traits are similar to the correlations among physiological traits. For example, at 12 °C, FA CaM and WAM were negatively correlated [18], and similarly, ME5_heart was significantly correlated with opposite signs with these two traits (i.e., positively correlated with FA CaM but negatively correlated with WAM, Table 3, Figs. 4, 5). Additionally, MEs correlated to WAM and CaM were enriched in KEGG metabolic pathways and GO terms related to metabolism. These data indicate that modules represent independent, biologically important mRNA networks.

The biological importance of co-expressed mRNA networks is also supported by their relation to metabolic processes. Eleven of the 13 significant heart or brain MEs were significantly enriched for at least one KEGG term, and 6 were significantly enriched for at least one GO term. KEGG terms mapped to biologically relevant KEGG pathways including metabolic pathways, mechanistic target of rapamycin (mTOR) signaling, mitogen activated protein kinase (MAPK) signaling, insulin signaling, and metabolism and biosynthesis of various macromolecules including glycogen, NADH precursors, amino sugars, and fatty acids (Table S3, S4). Importantly, 8 out of 11 modules with significantly enriched KEGG terms mapped to at least one metabolism related KEGG pathway. Among these KEGG terms are an abundance of signaling pathways that are known to affect physiological systems (see below). Yet, we also found several enriched KEGG pathways and GO terms that were uniquely enriched in only one or few modules and seemingly unrelated to the correlated trait(s) (e.g., cellular senescence). This could indicate a limited understanding of the complexity and interconnectedness among biological pathways and how different pathways affect a diversity of traits, mRNAs that are minimally annotated and missing relevant pathway involvement, or mRNA expression that impacts biological processes that indirectly impact the traits we have measured. The concept that there is a limited understanding of the interactions among pathways is supported by earlier mitochondrial respiration studies, specifically concerning the oxidative phosphorylation pathway (OxPhos) [78]. Among selectively important nuclear genes affecting OxPhos, none of the genes were among the 97 proteins in the OxPhos pathway; instead, they were in diverse pathways, some of which made sense (e.g., ADP transport- where ADP is a substrate for OxPhos) [78]. Thus, the few MEs associated with unexpected pathways may indicate a complexity in physiological traits where many pathways and the genes in these pathways affect trait variation.

Previously, data from our laboratory demonstrated that natural variation in substrate specific cardiac metabolism in F. heteroclitus could be explained by cardiac mRNA expression using microarray data [16]. Similar to the WGCNA approach presented here, the first principal component of mRNA expression from different metabolic pathways (oxidative phosphorylation, TCA cycle, glycolysis) explained substrate specific CaM among individuals with different mRNA pathways explaining substrate specific metabolisms in different individuals. Here, we found that mRNA expression explained a similar proportion of substrate specific CaM as previously reported (~ 80%) using three MEs to explain a single trait (FA CaM at 12 °C).

Few, if any, studies have examined the correlation between co-expressed mRNA and CTmax (although see [79]). Our analyses found that 341 heart mRNAs in two co-expressed modules and 682 brain mRNAs in three co-expressed modules were associated with CTmax at 12 °C or 28 °C, with different MEs at each temperature. Furthermore, heart and brain significant MEs for CTmax at 28 °C share enriched KEGG pathways, yet do not share any mRNAs, suggesting that different mRNAs affect a common set of pathways that impact CTmax. These data suggest that CTmax is polygenic and relies on different mRNAs in different tissues at different temperatures. There is prior evidence suggesting that CTmax is polygenic: a GBS study covering ~ 0.1% of the genome found up to 47 single nucleotide polymorphisms (SNPs) that explained 43.4% of variation in CTmax among F. heteroclitus individuals collected from Georgia, New Jersey, and New Hampshire, USA [28]. Here, a greater proportion of CTmax was explained with mRNA expression, up to 71.7% with 3 brain MEs. This increase in explained CTmax variance by mRNA expression is likely due to the combined heritable and physiologically inducible nature of mRNA expression. Few (0–10) of the mRNAs in MEs were differentially expressed between 12 °C and 28 °C, and thus MEs that explained CTmax variation within each of acclimation temperatures are not due to acclimation effects on mRNA expression. Instead, the CTmax variation within each acclimation temperature appears to be due to individual variation in mRNA expression, which may be explained by nucleotide variation driving differential expression.

Whole animal metabolism, WAM, is a fundamental physiological process that defines how animals live, niche space, evolutionary transition, and the human condition [2, 4, 6, 7, 80]. There is significant literature investigating metabolic rate variation (e.g. [41, 81, 82],); however, the relationship between metabolic rate and mRNA expression remains poorly understood. This may be due to the complex nature of whole animal metabolism, which is a sum of tissue specific metabolic demands and a balance between growth, maintenance, and energy storage. Yet, we find 82% of 12 °C WAM variation related to four heart MEs with 736 mRNAs and 50% of 28 °C WAM variation related to one brain ME with 142 mRNAs. These data indicate that a large proportion of WAM can be explained by mRNAs within a common pathway impacting cardiac metabolic processes and thus provides insight into the physiological relationship between cardiorespiratory performance and overall metabolism [83,84,85,86].

These WGCNA analyses suggest that many mRNAs in several biochemical pathways define the physiological state among individuals. Yet, the careful reader will note two substantial complexities: 1) heart MEs explain the variation in many physiological traits at 12 °C but few at 28 °C, and brain MEs explain the variation in many physiological traits at 28 °C but few at 12 °C and 2) mRNAs within MEs are enriched for many diverse and unexpected pathways as discussed above.

The difference between tissue specific MEs and their association with physiological traits is related to CaM, WAM, and CTmax having higher inter-individual variation at 12 °C than 28 °C, and similarly there is greater mRNA expression variation at 12 °C than at 28 °C. Thus, the more frequent explanation of physiological traits by 12 °C mRNA expression may simply result from greater statistical power due to the greater variance in both physiological traits and mRNA. Yet, in brains, mRNAs explain WAM and CTmax at 28 °C. While we can only speculate, these data suggest that at the higher temperature brain mRNA expression is more important than cardiac mRNA expression. There is evidence that acclimatory response to temperature in brain is greater than in hearts (more acclimatory mRNAs and greater decreased mitochondrial function in brain when compared to heart tissue [24]). This is similar to our data: brains at 28 °C have more acclimation responsive mRNAs than hearts and more population divergence than brains at 12 °C or hearts at either temperature. Together these data suggest that the variation in WAM and CTmax at 28 °C are more dependent on brain specific expression.

Many MEs are associated with KEGG pathways that impinge on metabolic processes. Physiological processes are within 9 heart MEs and 4 brain MEs (Tables 2 and 4), and these MEs each contain 90–554 mRNAs. Each of these MEs is significantly enriched for multiple KEGG and GO pathways (Tables S3 and S4). Surprisingly MEs that explain trait variance do not mainly include genes involved in primary metabolic pathways (e.g., glycolysis, TCA cycle, or oxidative phosphorylation) but instead are enriched for several signaling pathways including MAPK, mTOR, AMPK, PPAR, and FoxO. These signaling pathways are known to impact metabolic and thermal tolerance among ectotherms. For example, MAPK has been linked to adaptive cold tolerance [55, 87] and lipid metabolism [88]. Additionally, mTOR is involved with energy homeostasis, has been linked to growth and longevity, and may be sensitive to temperature variation [89,90,91]. AMPK induces cellular ATP production in mammals and is important for thermal stress response in ectotherms [92,93,94]. Furthermore, AMPK phosphorylation in Coho salmon and rainbow trout hearts has been correlated with exposure above optimum temperatures [95], suggesting a role of AMPK in fish thermal response. FoxO proteins, especially FoxO1, are involved in energy homeostasis and may aid in the switch from carbohydrate (glycolytic) to fatty acid metabolites [96]. Therefore, although these pathways may not be the “usual suspects”, their role in ectotherm metabolism and thermal tolerance is well documented, and we suggest that signaling pathways play an important role in explaining the trait variation examined here. In addition to these important signaling pathways, several modules are enriched for pathways not typically thought to be directly involved in metabolism or thermal tolerance. Similar to nuclear genes that impact Fundulus mitochondrial respiration [78], these data suggest that many metabolically distant genes affect physiological variation. This is important because too often publications have “just so stories” (e.g., [97, 98]) that only focus on a few preconceived genes to explain functional physiological variation [99]. While it is understandable to highlight a prior expectation, doing so limits our understanding of how genotypes effect phenotypes.

Conclusion

In summary, these data address two hypothesis, first, that genetically similar populations with little divergence in physiological traits would have a shared mRNA expression response to thermal acclimation. Here, we found that mRNAs important for acclimation are population specific and that divergence among geographically close populations did not include acclimation responsive mRNAs. This suggests that even genetically similar populations have distinct thermal response molecular mechanisms. Our second hypothesis was that physiological traits would be related to temperature and tissue specific mRNA expression, with variance among traits in the mRNAs that explain trait variation. We found this to be true, with tissue specific mRNA expression associated with physiological traits dependent on the thermal environment. We highlight that biologically important mRNA networks are related to 48–82% of variation in whole animal metabolism, thermal tolerance, or substrate specific cardiac metabolism and are different at different thermal environments. This suggests that mRNA variation among individuals within and among populations is important for explaining complex trait variation and, surprisingly, that while similar pathways can be important at different temperatures, the tissues where they are expressed differ: heart mRNA expression explains variation in more traits at 12 °C, and brain mRNA expression explains variation in more traits at 28 °C.

Methods

Animal care and use

Adult fish were collected in live traps in September 2018 at three sites in central New Jersey, USA near the Oyster Creek Nuclear generating station (OCNGS). Sites included one north reference (N.Ref; 39°52′28.000 N, 74°08′19.000 W), one south reference (S.Ref; 39°47′04.000 N, 74°11′07.000 W), and a central site located within the thermal effluent of the OCNGS (TE; 39°48′33.000 N, 74°10′51.000 W). All fish were transferred live to the University of Miami, FL where they were kept in accordance with the University of Miami Institutional Animal Care and Use Committee (IACUC) guidelines.

Individuals from all three populations were common gardened to 20 °C for three months (12:12 light dark cycle) and kept at 20 °C in a common recirculating seawater system (15ppt) at 12 h light:12 h dark, then subjected to pseudo-winter for 6 weeks at 8 °C (8:16 light dark cycle). Following the pseudo-winter, half of the fish from each population were acclimated to 12 °C and the other half to 28 °C (16:8 light dark cycle) for four weeks prior to determination of WAM and CTmax. Following this acclimation, fish originally acclimated to 12 °C were acclimated to 28 °C and vice versa for at least four weeks, and WAM and CTmax were measured at the new acclimation temperature. After a minimum one-week recovery period post- CTmax, fish were sacrificed, substrate specific CaM was measured at the second acclimation temperature, and mRNA was isolated. Thus, WAM and CTmax were measured in all individuals at 12 °C and 28 °C, but CaM and mRNA were sampled from half the individuals acclimated at 12 °C and the other half at 28 °C.

Quantifying metabolic and thermal tolerance traits

Individuals acclimated to 12 °C and 28 °C for at least 4 weeks were measured for whole animal metabolic rate (overnight intermittent flow respirometry) with a minimum of 20 replicate metabolic rate measures per individual used to determine the standard metabolic rate (SMR) in mgO2 hr−1 [34]. Critical thermal maximum (CTmax) was measured in a 10-gallon tank that was slowly heated at a rate of 0.3 °C min−1 as in [100] and was defined as the point when fish lost equilibrium in the water column for 5 consecutive seconds. Finally, substrate specific cardiac metabolic rate (CaM, substrates: 5 mM glucose, fatty acids – 1 mM Palmitic acid conjugated to fatty-acid-free bovine serum albumin, lactate + ketones + ethanol – 5 mM lactate, 5 mM hydroxybutyrate, 5 mM ethyl acetoacetate, 0.1% ethanol, endogenous – substrate free Ringers media) was measured using micro-respirometry. Heart ventricles were dissected out and splaying in Ringers media before being transferred to a custom 1 mL chamber system [101]. Each heart was measured for glucose (GLU), then fatty acids (FA), followed by lactate + ketones + ethanol (LKA), and endogenous (END) cardiac metabolic rate. All substrates except GLU used non-reversible glycolytic enzyme inhibitors, so the order of substrates did not differ among hearts. Each heart was measured for 6 min per substrate with the last 3 min used to calculated oxygen consumption in pmolO2 sec−1. All respirometry measurements used PreSens fiber optic oxygen sensors (POF) with flow through cells (FTC-Pst7-10,WAM) or sensor spots (SP-PSt7-10, CaM) (PreSens Precision Sensing, Regensburg, Germany). Both WAM and CaM were measured at both acclimation temperatures (12 °C and 28 °C) for the same individuals, and CaM was measured at a single acclimation temperature. For additional methods and analysis of physiological data see [18].

mRNA library preparation and sequencing

Tissues for mRNA expression were stored in chaotropic buffer at the time of CaM measurement and captured gene expression variation due to long term temperature acclimation rather than heat shock or acute temperature response (e.g., resulting from higher temperatures during CTmax measurements) because CTmax measurements were performed at least a week prior to tissue isolation. We extracted total RNA from homogenized heart and brain tissues using a phenol–chloroform isoamyl alcohol isolation and treated RNA samples with DNase to remove genomic DNA. For each sample we started with 50 ng of RNA and captured the 3’ mRNA ends using an NVdT primer with a poly-A tail for first strand cDNA synthesis (Table S5). This primer contained a unique barcode for each sample (1–96), which allowed all samples in a single plate to be pooled for the remaining library preparation steps. Nick translation was used to make double stranded cDNA that was digested with an in-house purified Tn5 transposase (as in [102]) loaded with partial adapter sequences to generate fragments of double stranded cDNA ranging from ~ 300-800 bp (Table S6). Libraries were amplified for 17 PCR cycles using primers complimentary to the inserted partial adapter sequence and a plate level barcode to fully multiplex samples.

mRNA data processing and analysis

A total of 219 libraries (110 individuals, 2 tissues per individual, 1 individual only heart was collected) were pooled and sequenced on 2 lanes of Illumina HiSeq4000 (dual end 150 bp reads) at the Genewiz LCC facility, South Plainfield NJ, USA. Raw reads were trimmed with BBDuk (from BBMap v38.87) to remove adapter sequences, aligned with STAR (v2.7.5) to the Fundulus heteroclitus genome, and counted with Featurecounts (v1.4.6-p5, parameters: -T 4 -s 2 -t gene -g gene_id).

The raw counts table was imported into R Studio (v1.4.1106) and all counts were normalized for library size using the median ratio method [103] with the “estimateSizeFactors” function in DESeq2 v1.6.3 [42]. Samples with a minimum of 1.5 million reads for hearts or 1 million reads for brains were retained and filtered to keep only mRNAs with at least 30 counts in 10% of individuals. After filtering, 53 heart samples and 58 brain samples remained. Principal component analysis (PCA) using “plotPCA” function from the DESeq2 package was used to examine variation among all samples. Among all samples, the 500 most variable mRNAs were used for PCA. In this analysis, 12 hearts and 13 brains were removed as outliers because they differed in expression from other same-tissue samples (i.e., some hearts had “brain-like” expression patterns and vice versa, Fig. S2). This reduced variation within a tissue and reduced sample size to 41 hearts and 45 brains (86 total individuals); however, the individuals removed were not from a single acclimation temperature or population so likely had little overall impact on further analyses. In a separate tissue specific principal component analysis, clustering of samples by biological effects including sex, habitat temperature, population, acclimation temperature, and date of tissue collection (possible batch effect) was examined. Batch effects did not split individuals along any of the principal components examined (PC1-PC4) for heart or brain. In this separate analysis, heart PC1 accounted for 18%, heart PC2 for 7%, brain PC1 for 11%, and brain PC2 for 7% of the variation among individuals. No biologically relevant clustering was detected among the first four principal components for either tissue (chi-squared test p > 0.05, Fig. S3A, B). To determine the degree of variation in mRNA expression among groups within a tissue, the coefficient of variation (CV, standard deviation/mean) for each expressed mRNA was calculated and the average CV compared among groups.

Differential expression analysis

DESeq2 v1.6.3 [42] package in R was used for differential expression analysis separately for heart and brain. To identify differentially expressed mRNAs between acclimation temperatures within populations, the DESeq model used was: ~ Population + Acclimation_Temperature + Population*Acclimation_Temperature. Additionally, due to significant interaction between population and acclimation temperature, a separate analysis was used to find differentially expressed mRNAs among populations within an acclimation temperature; individuals measured for CaM only at 12 °C or 28 °C were used with DESeq model: ~ Population. Multiple test correction across all comparisons made within a model used the Benjamin Hochberg false discovery rate with a significance threshold of 0.05.

Weighted gene co-expression network analysis

To identify sets of co-expressed mRNAs, weighted gene co-expression network analysis (WGCNA v1.70–3, [46]) was completed for heart and brain separately. Network calculation, used to group mRNAs into co-expressed modules for heart and brain, used soft thresholding to generate a scale free network with high similarity (soft thresholding power set to 5 in heart, 4 in brain) before calculating the topological overlap measure (TOM) and using dynamicTree with minimum module size of 30 and threshold for module merging of 0.75. The first principal component of each independent module (below the threshold for module merging), known as the module eigengene (ME), was then correlated with temperature specific quantitative traits using signed Pearson’s correlation. Multiple test correction across all correlations made for a single trait used the Benjamin Hochberg false discovery rate with a significance threshold of 0.05. To remove significant correlations potentially driven by outliers, a jack-knife approach was used to subsample 90% of individuals and repeat the signed Pearson correlation analysis 100 times. Correlations that were significant in > 70 out of 100 subsamples in the same direction were robust to outliers and reported as significant. A multivariate correlation coefficient was calculated for traits significantly correlated with more than one module by correlating the fitted values from a linear model with formula: trait ~ ME1 + ME2..ME# with trait data. This multivariate correlation coefficient represents the ability of the MEs together to accurately predict the trait. For significant modules, the mRNAs with the highest module membership (MM, correlation between mRNA expression and module eigengene) and the highest gene significance for a given trait within a module (GS, correlation between mRNA expression and a given quantitative trait) were also identified.

Gene ontology and Kyoto encyclopedia of genes and genomes enrichment

To identify biologically important networks within WGCNA modules that were significantly correlated with at least 1 trait, we used KEGG pathway and GO enrichment analyses. First, the genome was mapped to KEGG and GO terms using eggNOG mapper [104] with default parameters. The KEGG and GO terms were then matched to the set of expressed mRNAs for heart and brain. The list of mRNAs in each module was compared to the set of expressed mRNAs (set as the reference or gene universe) in each tissue for enrichment analysis in R using the clusterProfiler v3.16.1 package “enricher” function for KEGG terms [105]. To map enriched KEGG terms to KEGG pathways, the KEGG Mapper online tool was used with annotations from the closest relative, zebrafish (Danio rerio) [106, 107]. Cytoscape v3.8.2 with BiNGO v3.0.3 was used for GO enrichment using the set of expressed mRNAs as the reference to examined enrichment of biological process, molecular function, and cellular component GO terms [108]. Significant KEGG and GO terms are reported with FDR p-value threshold of 0.05.

Availability of data and materials

All sequence data is available in NCBI SRA: https://dataview.ncbi.nlm.nih.gov/object/PRJNA796010?reviewer=f1u5dbo46558lklv8m2cknap87 (public DOI available upon publication). All physiological data is available in Dryad: https://doi.org/10.5061/dryad.0gb5mkm0w. Code for processing of raw sequence files and all data analysis and visualization conducted in R is available in Github: https://github.com/mxd1288/OCNJ_F18_RNA.git.

References

  1. Bailly D, Cassemiro FA, Agostinho CS, Marques EE, Agostinho AA. The metabolic theory of ecology convincingly explains the latitudinal diversity gradient of neotropical freshwater fish. Ecology. 2014;95(2):553–62.

    Article  PubMed  Google Scholar 

  2. Brown JH, Gillooly JF, Allen AP, Savage VM, West GB. Toward a metabolic theory of ecology. Ecology. 2004;85(7):1771–89.

    Article  Google Scholar 

  3. Clarke A, Portner HO. Temperature, metabolic power and the evolution of endothermy. Biol Rev Camb Philos Soc. 2010;85(4):703–27.

    PubMed  Google Scholar 

  4. Killen SS, Glazier DS, Rezende EL, Clark TD, Atkinson D, Willener AS, Halsey LG. Ecological influences and morphological correlates of resting and maximal metabolic rates across teleost fish species. Am Nat. 2016;187(5):592–606.

    Article  PubMed  Google Scholar 

  5. Mathot KJ, Dingemanse NJ, Nakagawa S. The covariance between metabolic rate and behaviour varies across behaviours and thermal types: meta-analytic insights. Biol Rev Camb Philos Soc. 2019;94(3):1056–74.

    Article  PubMed  Google Scholar 

  6. Portner HO, Schulte PM, Wood CM, Schiemer F. Niche dimensions in fishes: an integrative view. Physiol Biochem Zool. 2010;83(5):808–26.

    Article  CAS  PubMed  Google Scholar 

  7. Wallace DC. A mitochondrial paradigm of metabolic and degenerative diseases, aging, and cancer: a dawn for evolutionary medicine. Annu Rev Genet. 2005;39:359–407.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Somero GN. Comparative physiology: a “crystal ball” for predicting consequences of global change. Am J Physiol Regul Integr Comp Physiol. 2011;301(1):R1-14.

    Article  CAS  PubMed  Google Scholar 

  9. Somero GN. Thermal physiology and vertical zonation of intertidal animals: Optima, limits, and costs of living. Integr Comp Biol. 2002;42(4):780–9.

    Article  PubMed  Google Scholar 

  10. Portner HO, Peck L, Somero G. Thermal limits and adaptation in marine Antarctic ectotherms: an integrative view. Philos Trans R Soc Lond B Biol Sci. 2007;362(1488):2233–58.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  11. Gibson G, Weir B. The quantitative genetics of transcription. Trends Genet. 2005;21(11):616–23.

    Article  CAS  PubMed  Google Scholar 

  12. Kliebenstein DJ. A role for gene duplication and natural variation of gene expression in the evolution of metabolism. PLoS ONE. 2008;3(3):e1838.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  13. McCairns RJ, Bernatchez L. Adaptive divergence between freshwater and marine sticklebacks: insights into the role of phenotypic plasticity from an integrated analysis of candidate gene expression. Evolution. 2010;64(4):1029–47.

    Article  CAS  PubMed  Google Scholar 

  14. Dayan DI, Crawford DL, Oleksiak MF. Phenotypic plasticity in gene expression contributes to divergence of locally adapted populations of Fundulus heteroclitus. Mol Ecol. 2015;24(13):3345–59.

    Article  PubMed  Google Scholar 

  15. Levine MT, Eckert ML, Begun DJ. Whole-genome expression plasticity across tropical and temperate drosophila melanogaster populations from Eastern Australia. Mol Biol Evol. 2011;28(1):249–56.

    Article  CAS  PubMed  Google Scholar 

  16. Oleksiak MF, Roach JL, Crawford DL. Natural variation in cardiac metabolism and gene expression in Fundulus heteroclitus. Nat Genet. 2005;37(1):67–72.

    Article  CAS  PubMed  Google Scholar 

  17. Kvist J, Wheat CW, Kallioniemi E, Saastamoinen M, Hanski I, Frilander MJ. Temperature treatments during larval development reveal extensive heritable and plastic variation in gene expression and life history traits. Mol Ecol. 2013;22(3):602–19.

    Article  CAS  PubMed  Google Scholar 

  18. Drown MK, DeLiberto AN, Ehrlich MA, Crawford DL, Oleksiak MF. Interindividual plasticity in metabolic and thermal tolerance traits from populations subjected to recent anthropogenic heating. R Soc Open Sci. 2021;8(7):210440.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Yeaman S. Local adaptation by alleles of small effect. Am Nat. 2015;186:S74–89.

    Article  PubMed  Google Scholar 

  20. Barghi N, Tobler R, Nolte V, Jakšić AM, Mallard F, Otte KA, Dolezal M, Taus T, Kofler R, Schlötterer C. Genetic redundancy fuels polygenic adaptation in Drosophila. PLoS Biol. 2019;17(2):e3000128.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Sidell BD, Johnston IA, Moerland TS, Goldspink G. The eurythermal myofibrillar protein complex of the mummichog (Fundulus heteroclitus): adaptation to a fluctuating thermal enviroment. J Comp Physiol. 1983;153(2):167–73.

    Article  CAS  Google Scholar 

  22. Baris TZ, Blier PU, Pichaud N, Crawford DL, Oleksiak MF. Gene by environmental interactions affecting oxidative phosphorylation and thermal sensitivity. Am J Physiol Regul Integr Comp Physiol. 2016;311(1):R157–65.

    Article  PubMed  Google Scholar 

  23. Bryant HJ, Chung DJ, Schulte PM. Subspecies differences in thermal acclimation of mitochondrial function and the role of uncoupling proteins in killifish. J Exp Biol. 2018;221(Pt 24):jeb186320.

    Article  PubMed  Google Scholar 

  24. Chung DJ, Bryant HJ, Schulte PM. Thermal acclimation and subspecies-specific effects on heart and brain mitochondrial performance in a eurythermal teleost (Fundulus heteroclitus). J Exp Biol. 2017;220(8):1459–71.

    PubMed  Google Scholar 

  25. Dhillon RS, Schulte PM. Intraspecific variation in the thermal plasticity of mitochondria in killifish. J Exp Biol. 2011;214(Pt 21):3639–48.

    Article  CAS  PubMed  Google Scholar 

  26. Fangue NA, Hofmeister M, Schulte PM. Intraspecific variation in thermal tolerance and heat shock protein gene expression in common killifish, Fundulus heteroclitus. J Exp Biol. 2006;209(15):2859–72.

    Article  CAS  PubMed  Google Scholar 

  27. Healy TM, Schulte PM. Thermal acclimation is not necessary to maintain a wide thermal breadth of aerobic scope in the common killifish (Fundulus heteroclitus). Physiol Biochem Zool. 2012;85(2):107–19.

    Article  CAS  PubMed  Google Scholar 

  28. Healy TM, Brennan RS, Whitehead A, Schulte PM. Tolerance traits related to climate change resilience are independent and polygenic. Global Change Biol. 2018;24(11):5348–60.

    Article  Google Scholar 

  29. Pierce VA, Crawford DL. Phylogenetic analysis of thermal acclimation of the glycolytic enzymes in the genus Fundulus. Physiol Zool. 1997;70(6):597–609.

    Article  CAS  PubMed  Google Scholar 

  30. Healy TM, Schulte PM. Factors affecting plasticity in whole-organism thermal tolerance in common killifish (Fundulus heteroclitus). J Comp Physiol B. 2012;182(1):49–62.

    Article  PubMed  Google Scholar 

  31. Dayan DI, Du X, Baris TZ, Wagner DN, Crawford DL, Oleksiak MF. Population genomics of rapid evolution in natural populations: polygenic selection in response to power station thermal effluents. BMC Evol Biol. 2019;19(1):61.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Duvernell DD, Lindmeier JB, Faust KE, Whitehead A. Relative influences of historical and contemporary forces shaping the distribution of genetic variation in the Atlantic killifish Fundulus heteroclitus. Mol Ecol. 2008;17(5):1344–60.

    Article  PubMed  Google Scholar 

  33. Wagner DN, Baris TZ, Dayan DI, Du X, Oleksiak MF, Crawford DL. Fine-scale genetic structure due to adaptive divergence among microhabitats. Heredity (Edinb). 2017;118(6):594–604.

    Article  CAS  Google Scholar 

  34. Drown MK, DeLiberto AN, Crawford DL, Oleksiak MF. An innovative setup for high-throughput respirometry of small aquatic animals. Front Mar Sci. 2020;7.

  35. White CR, Alton LA, Frappell PB. Metabolic cold adaptation in fishes occurs at the level of whole animal, mitochondria and enzyme. Proc R Soc B Biol Sci. 2012;279(1734):1740–7.

    Article  CAS  Google Scholar 

  36. Cumillaf JP, Blanc J, Paschke K, Gebauer P, Díaz F, Re D, Chimal ME, Vásquez J, Rosas C. Thermal biology of the sub-polar–temperate estuarine crab Hemigrapsus crenulatus (Crustacea: Decapoda: Varunidae). Biology Open. 2016;5(3):220–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Morgan R, Finnoen MH, Jutfelt F. CTmax is repeatable and doesn’t reduce growth in zebrafish. Sci Rep. 2018;8(1):7099.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  38. Ronning B, Jensen H, Moe B, Bech C. Basal metabolic rate: heritability and genetic correlations with morphological traits in the zebra finch. J Evol Biol. 2007;20(5):1815–22.

    Article  CAS  PubMed  Google Scholar 

  39. Mattila AL, Hanski I. Heritability of flight and resting metabolic rates in the Glanville fritillary butterfly. J Evol Biol. 2014;27(8):1733–43.

    Article  CAS  PubMed  Google Scholar 

  40. Nilsson JA, Akesson M, Nilsson JF. Heritability of resting metabolic rate in a wild population of blue tits. J Evol Biol. 2009;22(9):1867–74.

    Article  PubMed  Google Scholar 

  41. Reemeyer JE, Rees BB. Plasticity, repeatability and phenotypic correlations of aerobic metabolic traits in a small estuarine fish. J Exp Biol. 2020;223(Pt 14):jeb228098.

    Article  PubMed  Google Scholar 

  42. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  43. Orczewska JI, Hartleben G, O’Brien KM. The molecular basis of aerobic metabolic remodeling differs between oxidative muscle and liver of threespine sticklebacks in response to cold acclimation. Am J Physiol Regul Integr Comp Physiol. 2010;299(1):R352–364.

    Article  CAS  PubMed  Google Scholar 

  44. Podrabsky JE, Somero GN. Changes in gene expression associated with acclimation to constant temperatures and fluctuating daily temperatures in an annual killifish Austrofundulus limnaeus. J Exp Biol. 2004;207(Pt 13):2237–54.

    Article  CAS  PubMed  Google Scholar 

  45. Williams LM, Oleksiak MF. Signatures of selection in natural populations adapted to chronic pollution. BMC Evol Biol. 2008;8:282.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  46. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  47. Han L, Shen WJ, Bittner S, Kraemer FB, Azhar S. PPARs: regulators of metabolism and as therapeutic targets in cardiovascular disease Part I: PPAR-alpha. Future Cardiol. 2017;13(3):259–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Visscher PM, Brown MA, McCarthy MI, Yang J. Five years of GWAS discovery. Am J Hum Genet. 2012;90(1):7–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Visscher PM, Goddard ME. From R.A. fisher’s 1918 paper to GWAS a century later. Genetics. 2019;211(4):1125–30.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Hancock AM, Witonsky DB, Ehler E, Alkorta-Aranburu G, Beall C, Gebremedhin A, Sukernik R, Utermann G, Pritchard J, Coop G, et al. Colloquium paper: human adaptations to diet, subsistence, and ecoregion are due to subtle shifts in allele frequency. Proc Natl Acad Sci U S A. 2010;107(Suppl 2):8924–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Koch RE, Buchanan KL, Casagrande S, Crino O, Dowling DK, Hill GE, Hood WR, McKenzie M, Mariette MM, Noble DWA, et al. Integrating mitochondrial aerobic metabolism into ecology and evolution. Trends Ecol Evol. 2021;36(4):321–32.

    Article  PubMed  Google Scholar 

  52. Scheffers BR, De Meester L, Bridge TC, Hoffmann AA, Pandolfi JM, Corlett RT, Butchart SH, Pearce-Kelly P, Kovacs KM, Dudgeon D, et al. The broad footprint of climate change from genes to biomes to people. Science. 2016;354(6313):aaf7671.

    Article  PubMed  CAS  Google Scholar 

  53. Spicer JI, Morley SA, Bozinovic F. Physiological diversity, biodiversity patterns and global climate change: testing key hypotheses involving temperature and oxygen. Philos Trans R Soc Lond B Biol Sci. 2019;374(1778):20190032.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Schulte PM, Healy TM, Fangue NA. Thermal performance curves, phenotypic plasticity, and the time scales of temperature exposure. Integr Comp Biol. 2011;51(5):691–702.

    Article  PubMed  Google Scholar 

  55. Wen X, Zhang XY, Hu YD, Xu JJ, Wang T, Yin SW. iTRAQ-based quantitative proteomic analysis of Takifugu fasciatus liver in response to low-temperature stress. J Proteomics. 2019;201:27–36.

    Article  CAS  PubMed  Google Scholar 

  56. Zhou S, Campbell TG, Stone EA, Mackay TF, Anholt RR. Phenotypic plasticity of the Drosophila transcriptome. Plos Genet. 2012;8(3):e1002593.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Fangue NA, Richards JG, Schulte PM. Do mitochondrial properties explain intraspecific variation in thermal tolerance? J Exp Biol. 2009;212(Pt 4):514–22.

    Article  CAS  PubMed  Google Scholar 

  58. Lucassen M, Schmidt A, Eckerle LG, Portner H-O. Mitochondrial proliferation in the permanent vs temporary cold: enzyme activities and mRNA levels in Antarctic and temperate zoarcid fish. Am J Physiol Regul Integr Comp Physiol. 2003;285(6):R1410–20.

    Article  CAS  PubMed  Google Scholar 

  59. Chapman RW, Mancia A, Beal M, Veloso A, Rathburn C, Blair A, Holland AF, Warr GW, Didinato G, Sokolova IM, et al. The transcriptomic responses of the eastern oyster, Crassostrea virginica, to environmental conditions. Mol Ecol. 2011;20(7):1431–49.

    Article  PubMed  Google Scholar 

  60. Windisch HS, Frickenhaus S, John U, Knust R, Pörtner H-O, Lucassen M. Stress response or beneficial temperature acclimation: transcriptomic signatures in Antarctic fish (Pachycara brachycephalum). Mol Ecol. 2014;23(14):3469–82.

    Article  CAS  PubMed  Google Scholar 

  61. Somero GN. Proteins and temperature. Annu Rev Physiol. 1995;57:43–68.

    Article  CAS  PubMed  Google Scholar 

  62. Gilad Y, Rifkin SA, Pritchard JK. Revealing the architecture of gene regulation: the promise of eQTL studies. Trends Genet. 2008;24(8):408–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Hill MS, Vande Zande P, Wittkopp PJ. Molecular and evolutionary processes generating variation in gene expression. Nat Rev Genet. 2021;22(4):203–15.

    Article  CAS  PubMed  Google Scholar 

  64. Kelley JL, Tobler M, Beck D, Sadler-Riggleman I, Quackenbush CR, Rodriguez LA, Skinner MK. Epigenetic inheritance of DNA methylation changes in fish living in hydrogen sulfide–rich springs. Proc Natl Acad Sci. 2021;118(26):e2014929118.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Hu J, Wuitchik SJS, Barry TN, Jamniczky HA, Rogers SM, Barrett RDH. Heritability of DNA methylation in threespine stickleback (Gasterosteus aculeatus). Genetics. 2021;217(1):1–15.

    Article  PubMed  Google Scholar 

  66. Metzger DCH, Schulte PM. Epigenomics in marine fishes. Mar Genomics. 2016;30:43–54.

    Article  PubMed  Google Scholar 

  67. Tang Y, Ke ZP, Peng YG, Cai PT. Co-expression analysis reveals key gene modules and pathway of human coronary heart disease. J Cell Biochem. 2018;119(2):2102–9.

    Article  CAS  PubMed  Google Scholar 

  68. Liang W, Sun F, Zhao Y, Shan L, Lou H. Identification of susceptibility modules and genes for cardiovascular disease in diabetic patients using WGCNA analysis. J Diabetes Res. 2020;2020:4178639.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  69. Di Y, Chen D, Yu W, Yan L. Bladder cancer stage-associated hub genes revealed by WGCNA co-expression network analysis. Hereditas. 2019;156:7.

    Article  PubMed  PubMed Central  Google Scholar 

  70. Zhai X, Xue Q, Liu Q, Guo Y, Chen Z. Colon cancer recurrenceassociated genes revealed by WGCNA coexpression network analysis. Mol Med Rep. 2017;16(5):6499–505.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Su R, Jin C, Zhou L, Cao Y, Kuang M, Li L, Xiang J. Construction of a ceRNA network of hub genes affecting immune infiltration in ovarian cancer identified by WGCNA. BMC Cancer. 2021;21(1):970.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Yin X, Wang P, Yang T, Li G, Teng X, Huang W, Yu H. Identification of key modules and genes associated with breast cancer prognosis using WGCNA and ceRNA network analysis. Aging (Albany NY). 2020;13(2):2519–38.

    Article  Google Scholar 

  73. Niemira M, Collin F, Szalkowska A, Bielska A, Chwialkowska K, Reszec J, Niklinski J, Kwasniewski M, Kretowski A. Molecular signature of subtypes of non-small-cell lung cancer by large-scale transcriptional profiling: identification of key modules and genes by Weighted Gene Co-Expression Network Analysis (WGCNA). Cancers (Basel). 2019;12(1):37.

    Article  CAS  Google Scholar 

  74. Chen M, Yan J, Han Q, Luo J, Zhang Q. Identification of hub-methylated differentially expressed genes in patients with gestational diabetes mellitus by multi-omic WGCNA basing epigenome-wide and transcriptome-wide profiling. J Cell Biochem. 2020;121(5–6):3173–84.

    Article  CAS  PubMed  Google Scholar 

  75. Huang Z, Ma A, Yang S, Liu X, Zhao T, Zhang J, Wang XA, Sun Z, Liu Z, Xu R. Transcriptome analysis and weighted gene co-expression network reveals potential genes responses to heat stress in turbot Scophthalmus maximus. Comp Biochem Physiol Part D Genomics Proteomics. 2020;33:100632.

    Article  CAS  PubMed  Google Scholar 

  76. Zhang L, Zhang Q, Li W, Zhang S, Xi W. Identification of key genes and regulators associated with carotenoid metabolism in apricot (Prunus armeniaca) fruit using weighted gene coexpression network analysis. BMC Genomics. 2019;20(1):876.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Traylor-Knowles N, Connelly MT, Young BD, Eaton K, Muller EM, Paul VJ, Ushijima B, DeMerlis A, Drown MK, Goncalves A et al: Gene Expression Response to Stony Coral Tissue Loss Disease Transmission in M. cavernosa and O. faveolata From Florida. Front Mar Sci. 2021;8.

  78. Baris TZ, Wagner DN, Dayan DI, Du X, Blier PU, Pichaud N, Oleksiak MF, Crawford DL. Evolved genetic and phenotypic differences due to mitochondrial-nuclear interactions. Plos Genet. 2017;13(3):e1006517.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  79. Campbell-Staton SC, Velotta JP, Winchell KM. Selection on adaptive and maladaptive gene expression plasticity during thermal adaptation to urban heat islands. Nat Commun. 2021;12(1):6195.

    Article  PubMed  PubMed Central  Google Scholar 

  80. Okie JG, Smith VH, Martin-Cereceda M. Major evolutionary transitions of life, metabolic scaling and the number and size of mitochondria and chloroplasts. Proc Biol Sci. 1831;2016:283.

    Google Scholar 

  81. Pettersen AK, Marshall DJ, White CR. Understanding variation in metabolic rate. J Exp Biol. 2018;221(1):jeb166876.

    Article  PubMed  Google Scholar 

  82. Schulte PM. The effects of temperature on aerobic metabolism: towards a mechanistic understanding of the responses of ectotherms to a changing environment. J Exp Biol. 2015;218(Pt 12):1856–66.

    Article  PubMed  Google Scholar 

  83. Clark TD, Jeffries KM, Hinch SG, Farrell AP. Exceptional aerobic scope and cardiovascular performance of pink salmon (Oncorhynchus gorbuscha) may underlie resilience in a warming climate. J Exp Biol. 2011;214(Pt 18):3074–81.

    Article  CAS  PubMed  Google Scholar 

  84. Clark TD, Ryan T, Ingram BA, Woakes AJ, Butler PJ, Frappell PB. Factorial aerobic scope is independent of temperature and primarily modulated by heart rate in exercising Murray cod (Maccullochella peelii peelii). Physiol Biochem Zool. 2005;78(3):347–55.

    Article  CAS  PubMed  Google Scholar 

  85. Jensen DL, Overgaard J, Wang T, Gesser H, Malte H. Temperature effects on aerobic scope and cardiac performance of European perch (Perca fluviatilis). J Therm Biol. 2017;68(Pt B):162–9.

    Article  PubMed  Google Scholar 

  86. Nyboer EA, Chapman LJ. Cardiac plasticity influences aerobic performance and thermal tolerance in a tropical, freshwater fish at elevated temperatures. J Exp Biol. 2018;221(Pt 15):jeb17808.

    Google Scholar 

  87. Chen ZZ, Cheng CHC, Zhang JF, Cao LX, Chen L, Zhou LH, Jin YD, Ye H, Deng C, Dai ZH, et al. Transcrintomic and genomic evolution under constant cold in Antarctic notothenioid fish. P Natl Acad Sci USA. 2008;105(35):12944–9.

    Article  CAS  Google Scholar 

  88. Bost F, Aouadi M, Caron L, Binetruy B. The role of MAPKs in adipocyte differentiation and obesity. Biochimie. 2005;87(1):51–6.

    Article  CAS  PubMed  Google Scholar 

  89. Aramburu J, Ortells MC, Tejedor S, Buxade M, Lopez-Rodriguez C. Transcriptional regulation of the stress response by mTOR. Sci Signal. 2014;7(332):re2.

    Article  PubMed  CAS  Google Scholar 

  90. Sengupta S, Peterson TR, Laplante M, Oh S, Sabatini DM. mTORC1 controls fasting-induced ketogenesis and its modulation by ageing. Nature. 2010;468(7327):1100–U1502.

    Article  CAS  PubMed  Google Scholar 

  91. Reiling JH, Sabatini DM. Stress and mTORture signaling. Oncogene. 2006;25(48):6373–83.

    Article  CAS  PubMed  Google Scholar 

  92. Seebacher F, Little AG. Plasticity of performance curves can buffer reaction rates from body temperature variation in active endotherms. Front Physiol. 2017;8:575.

    Article  PubMed  PubMed Central  Google Scholar 

  93. Wittmann AC, Benrabaa SAM, Lopez-Ceron DA, Chang ES, Mykles DL. Effects of temperature on survival, moulting, and expression of neuropeptide and mTOR signalling genes in juvenile Dungeness crab (Metacarcinus magister). J Exp Biol. 2018;221(21):jeb187492.

    Article  PubMed  Google Scholar 

  94. Frederich M, O’Rourke MR, Furey NB, Jost JA. AMP-activated protein kinase (AMPK) in the rock crab, cancer irroratus: an early indicator of temperature stress. J Exp Biol. 2009;212(5):722–30.

    Article  CAS  PubMed  Google Scholar 

  95. Anttila K, Casselman MT, Schulte PM, Farrell AP. Optimum temperature in Juvenile Salmonids: connecting subcellular indicators to tissue function and whole-organism thermal optimum. Physiol Biochem Zool. 2013;86(2):245–56.

    Article  CAS  PubMed  Google Scholar 

  96. Gross DN, van den Heuvel APJ, Birnbaum MJ. The role of FoxO in the regulation of metabolism. Oncogene. 2008;27(16):2320–36.

    Article  CAS  PubMed  Google Scholar 

  97. Gould SJ. The Panda’s Thumb. New York/London: W. W. Norton & Company, Inc; 1992.

    Google Scholar 

  98. Gould SJ, Lewontin RC. The spandrels of San Marco and the Panglossian paradigm: a critique of the adaptationist programme. Proc R Soc Lond B Biol Sci. 1979;205(1161):581–98.

    Article  CAS  PubMed  Google Scholar 

  99. Linquist S, Doolittle WF, Palazzo AF. Getting clear about the F-word in genomics. Plos Genet. 2020;16(4):e1008702.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  100. Becker CD, Genoway RG. Evaluation of the critical thermal maximum for determining thermal tolerance of freshwater fish. Environ Biol Fishes. 1979;4:245–56.

    Article  Google Scholar 

  101. DeLiberto AN, Drown MK, Oleksiak MF, Crawford DL. Measuring complex phenotypes: a flexible high-throughput design for micro-respirometry. bioRxiv 2020.

  102. Picelli S, Bjorklund AK, Reinius B, Sagasser S, Winberg G, Sandberg R. Tn5 transposase and tagmentation procedures for massively scaled sequencing projects. Genome Res. 2014;24(12):2033–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  103. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  104. Huerta-Cepas J, Forslund K, Coelho LP, Szklarczyk D, Jensen LJ, von Mering C, Bork P. Fast genome-wide functional annotation through orthology assignment by eggNOG-mapper. Mol Biol Evol. 2017;34(8):2115–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  105. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (N Y). 2021;2(3):100141.

    Google Scholar 

  106. Kanehisa M, Sato Y. KEGG mapper for inferring cellular functions from protein sequences. Protein Sci. 2020;29(1):28–35.

    Article  CAS  PubMed  Google Scholar 

  107. Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2016;44(D1):D457–462.

    Article  CAS  PubMed  Google Scholar 

  108. Maere S, Heymans K, Kuiper M. BiNGO: a Cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks. Bioinformatics. 2005;21(16):3448–9.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

Thank you to Amanda DeLiberto and Samantha Sierra-Martinez for assistance in isolation of the Tn5 enzyme used for RNA library preparation and to Dr. Benjamin Young for support in RNA data analysis.

Funding

This research was supported by funding from National Science Foundation award nos. IOS 1556396 and IOS 1754437 to MFO and DLC.

Author information

Authors and Affiliations

  1. University of Miami Rosenstiel School of Marine and Atmospheric Science, 6400 Rickenbacker Causeway, Miami, FL, USA

    Melissa K. Drown, Douglas L. Crawford & Marjorie F. Oleksiak

Authors
  1. Melissa K. Drown
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Douglas L. Crawford
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Marjorie F. Oleksiak
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Data collection, analysis, and visualization MKD. Original manuscript draft MKD. Manuscript editing MFO and DLC. Study design MKD, MFO, and DLC. The author(s) read and approved the final manuscript.

Corresponding author

Correspondence to Melissa K. Drown.

Ethics declarations

Ethics approval and consent to participate

All animal use was approved by the University of Miami Institutional Animal Care and Use Committee (IACUC) and done in accordance with the University of Miami IACUC guidelines. All methods are reported in accordance with ARRIVE guidelines for the reporting of animal experiments.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1:

Figure S1. Trait variance in metabolic rate among ectotherms. A) Standard metabolic rate variance among species from different climates are compared to average variance in metabolic rate within the Fundulus heteroclitus populations used in this study. Variance calculated as mass corrected maximum (in warmer environment) – minimum (in cooler environment)/minimum (in cooler environment) (range spread).All standard metabolic rates are corrected for body mass (residual of metabolic rate vs. body mass + metabolic rate of an average sized fish from this data set). B) Number of species per group. Data from (1). Figure S2. Principal component analysis of all samples. Principal component 1 split heart and brain tissue and explained 86% of variance. Individuals who did not clearly group with the appropriate tissue were removed as outliers. Figure S3. Tissue specific principal component analysis. Heart (N=41, A) first two principal components explain 19% and 7% of variance. Brain (N=45, B) first two principal components explain 11% and 7% of variance. Triangles are 28°C acclimated individuals, circles are 12°C. Individuals from the north reference (N.Ref) are blue, south reference (S.Ref) are purple, and thermal effluent population (TE) are red. Figure S4. Differentially expressed mRNAs among populations within tissue and acclimation temperature. A) Heart at 12°C, B) brain at 12°C, C) heart at 28°C, D) brain at 28°C. Populations are north reference (N.Ref), south reference (S.Ref), and thermal effluent (TE).

Additional file 2:

Table S1. Significant Acclimation effect on mRNA Expression within each Population. Population differential expression.with FDR p<0.05.  Differentially expressed genes among population pairs within 12°C or 28°C with FDR p<0.05. Table S2. Population differential expression for each acclimation temperature with FDR <0.05 . Differentially expressed genes between 12°C and 28°C within populations with FDR p<0.05.Temperature differential expression. Differentially expressed genes between 12°C and 28°C within populations with FDR p<0.05. Table S3. Enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathways within Significant Modules. Module identifier (Module), correlated trait(s), enriched KEGG terms, and KEGG pathways are listed in columns. Enriched KEGG terms with FDR p-value < 0.05 were mapped to KEGG pathways using KEGG Mapper. Terms of interest based on relationship to correlated trait(s) are underlined. Table S4. Enriched Gene Ontology (GO) Terms within Significant Modules. Modules are named based on the tissue they are in and trait they are significantly correlated with. Enriched GO terms with FDR p-value < 0.05 are listed. Terms of interest based on relationship to correlated trait(s) are underlined. Table S5. Primers and barcodes for 3’ mRNA library preparation. First strand synthesis with i7 primers added i7 indicies. Second strand synthesis added i5 index to 96 pooled samples. Table S6. Tn5 loading and PCR amplification primer sequences. The lower-case nucleotides complement the Tn5MErev primer. The upper-case underlined nucleotides match with the i5 (R1) and P7 match with the i7 (part of the NVdT primer). 

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Drown, M.K., Crawford, D.L. & Oleksiak, M.F. Transcriptomic analysis provides insights into molecular mechanisms of thermal physiology. BMC Genomics 23, 421 (2022). https://doi.org/10.1186/s12864-022-08653-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-022-08653-y

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us