TRANS-ACTING GENOTYPES DRIVE M RNA EXPRESSION AFFECTING 1 METABOLIC AND THERMAL TOLERANCE TRAITS 2

Evolutionary processes driving physiological trait variation depend on the underlying genomic mechanisms. Evolution of these mechanisms depends on whether traits are genetically complex (involving many genes) and how gene expression that impact the traits is converted to phenotype. Yet, genomic mechanisms that impact physiological traits are diverse and context dependent (e.g., vary by environment or among tissues), making them difficult to discern. Here we examine the relationships between genotype, mRNA expression, and physiological traits to discern the genetic complexity and whether the gene expression effecting the physiological traits is primarily cis or trans-acting. We use low-coverage whole genome sequencing and tissue specific mRNA expression among individuals to identify polymorphisms directly associated with physiological traits and expressed quantitative trait loci (eQTL) driving variation in six temperature specific physiological traits (standard metabolic rate, thermal tolerance, and four substrate specific cardiac metabolic rates). Not surprisingly, there were few, only five, SNPs directly associated with physiological traits. Yet, by focusing on a select set of mRNAs belonging to co-expression modules that explain up to 82% of temperature specific (12°C or 28°C) metabolism and thermal tolerance, we identified hundreds of significant eQTL for mRNA whose expression affects physiological traits. Surprisingly, most eQTL (97.4% for heart and 96.7% for brain) of eQTL were trans-acting. This could be due to higher effect size or greater importance of trans versus cis acting eQTLs for mRNAs that are central to co-expression modules. That is, we may have enhanced the identification of trans-acting factors by looking for SNPs associated with mRNAs in co-expression modules that are known to be correlated with the expression of 10s or 100s of other genes, and thus have identified eQTLs with widespread effects on broad gene expression patterns. Overall, these data indicate that the genomic mechanism driving physiological variation across environments is driven by trans-acting tissue specific mRNA expression. Author Summary In the salt marsh killifish Fundulus heteroclitus there is amazingly large variation in physiological traits assumed to be under stabilizing selection, which should reduce their variation. To discern the heritability of this physiological variation we took an innovative approach to define the DNA variation that drives mRNA expression linked to physiological variation. This indirect approach revealed many DNA sequence variants associated with physiological variation via their effect on mRNA expression. Surprisingly, these changes were not in the mRNAs themselves, but in unlinked distant genes which regulate mRNA expression. That is, the vast majority (>95%) were trans-acting. This is surprising because trans-acting effects are found less often than DNA variants within or close to mRNA expression genes. Our results are likely related to the select subset of mRNAs across environments that are linked to physiological variation.

ABSTRACT 10 Evolutionary processes driving physiological trait variation depend on the underlying 11 genomic mechanisms. Evolution of these mechanisms depends on whether traits are genetically 12 complex (involving many genes) and how gene expression that impact the traits is converted to 13 phenotype. Yet, genomic mechanisms that impact physiological traits are diverse and context 14 dependent (e.g., vary by environment or among tissues), making them difficult to discern. Here 15 we examine the relationships between genotype, mRNA expression, and physiological traits to 16 discern the genetic complexity and whether the gene expression effecting the physiological traits 17 is primarily cis or trans-acting. We use low-coverage whole genome sequencing and tissue 18 specific mRNA expression among individuals to identify polymorphisms directly associated with 19 physiological traits and expressed quantitative trait loci (eQTL) driving variation in six 20 temperature specific physiological traits (standard metabolic rate, thermal tolerance, and four 21 substrate specific cardiac metabolic rates). Not surprisingly, there were few, only five, SNPs 22 directly associated with physiological traits. Yet, by focusing on a select set of mRNAs 23 belonging to co-expression modules that explain up to 82% of temperature specific (12°C or 24 28°C) metabolism and thermal tolerance, we identified hundreds of significant eQTL for mRNA 25 whose expression affects physiological traits. Surprisingly, most eQTL (97.4% for heart and 26 96.7% for brain) of eQTL were trans-acting. This could be due to higher effect size or greater 27 importance of trans versus cis acting eQTLs for mRNAs that are central to co-expression 28 modules. That is, we may have enhanced the identification of trans-acting factors by looking for 29 SNPs associated with mRNAs in co-expression modules that are known to be correlated with the 30 expression of 10s or 100s of other genes, and thus have identified eQTLs with widespread effects 31 on broad gene expression patterns. Overall, these data indicate that the genomic mechanism 32 driving physiological variation across environments is driven by trans-acting tissue specific 33 mRNA expression. These attributes make it difficult to predict or identify the genetic variation driving physiological 71 variation. One approach to simplify this multifaceted complexity is to identify the genomic 72 mechanism affecting mRNA expression that drives phenotypic variation. 73 mRNA expression variation is often biologically important in that complex or multivariate 74 mRNA expression can explain variation in a diverse suite of traits including thermal tolerance, 75 disease response, and metabolism (Zhang et  physiologically induced; yet mRNA expression also is heritable (Gibson and Weir 2005) and has 78 large variation among common gardened individuals (Oleksiak et al. 2002). Thus, it is possible 79 to identify heritable genetic loci associated with mRNA expression variation. These associations 80 between genetic loci and mRNA expression are identified as expression quantitative trait loci 81 (eQTL), where eQTLs mediate the expression of one or many genes. Furthermore, when eQTL 82 controlled mRNA expression impacts physiological traits, eQTLs may be evolutionary targets 83 for adaptation (Whitehead and Crawford 2006). 84 Here, we apply association studies to identify genetic loci directly (genome wide 85 association, GWAS) and indirectly (eQTL) driving physiological variation. We use three 86 common gardened wild populations to capture natural variation in six complex physiological 87 traits: whole animal metabolism (standard metabolic rate, SMR), critical thermal maximum 88 (CTmax), and four substrate specific cardiac metabolic rates (MRcardiac) ). 89 Physiological traits and heart and brain mRNA expression were quantified under two 90 ecologically relevant acclimation temperatures (12°C and 28°C) in the same individuals. We 91 found little evidence of population divergence in physiological traits  or 92 differentially expressed mRNAs (Drown et al. 2022). Therefore, we treated all individuals as 93 belonging to a single population and, using a tissue specific (heart and brain) weighted gene co-94 expression network analysis (WGCNA, (Langfelder and

Population structure 119
To determine the genetic structure among populations, we conducted an admixture 120 analysis. Using NGSadmix we tested seven K values where K is the number of ancestral 121 populations. We found K=4 to be the best fit based on log-likelihood probability with no clear 122 structure among populations (Fig. S1). These analyses, as well as physiological ( Horvath 2008)) to identify co-expression mRNAs and group them into modules (MEs). Single 141 mRNA expression was limited to the top 10 mRNAs from each WGCNA co-expression module. 142 There were 80 MEs: 39 heart modules and 41 brain modules (Table S1), 390 single heart  143 mRNAs, and 410 single brain mRNAs (Drown et al. 2022). 144 All association results are reported from ANGSD score test (-doAsso 2) and after multiple 145 test correction of p-values (Benjamini and Hochberg 1995). Physiological traits and mRNA 146 expression were measured at two acclimation temperatures (12°C and 28°C). Acclimation 147 temperature was included as a variable in the association tests and yielded similar results when 148 compared to within temperature associations. Sample size varied among association tests based 149 on phenotypic and genotypic data availability and are reported in Table S2. 150 1. Direct genotypic association to whole animal and tissue level metabolic and 151 thermal tolerance traits 152 A total of five independent SNPs were significantly associated with three traits: GLU 153 MRcardiac, LKA MRcardiac, and END MRcardiac (Table 1, Fig. 2). Of these five significant 154 associations, two were associated with GLU MRcardiac, one was associated with LKA MRcardiac, 155 and the remaining two were associated with END MRcardiac (Table 1, Fig. 2). Only one of these 156 SNPs (erfl1, directly associated with END MRcardiac) is also a significant eQTL. We did not find 157 significantly associated SNPs for SMR or CTmax, however, these traits were previously 158 associated with at least one mRNA co-expression module (Drown et al. 2022), described below. 159 The low number of SNPs directly associated with physiological traits is most likely due to the 160 small sample size: average sample size was 29.60  6.23 (mean  one standard deviation) 161 individuals per association test (Table S2). 162 We previously identified tissue specific mRNA co-expression modules for heart and brain 166 that are associated with variation in the six physiological traits (Drown et al. 2022). Co-167 expression modules included 39 significant heart modules with 90-554 mRNAs per module and 168 41 significant brain modules with 142-393 mRNAs per module (Drown et al. 2022). Each 169 module was assigned a module eigengene (ME, the first principal component of multivariate 170 mRNA expression), and each mRNA in the module was assigned a module membership defined 171 7 as the correlation coefficient between that mRNA and the ME. From each module we choose the 172 top 10 single mRNAs (based on module membership) and used these single mRNA expression 173 values in a series of association tests to identify eQTL. Additionally, we used the ME for each 174 module as a phenotypic value in a separate set of association studies to find SNPs associated 175 with multivariate mRNA expression (eQTLME). This allowed us to identify SNPs that explain 176 mRNA expression patterns previously correlated to acclimation temperature specific 177 physiological traits (SMR, CTmax, and MRcardiac measured at 12°C and 28°C) (Drown et al. 178 2022). Notably, we did not test all possible mRNAs (~10,000 mRNAs per tissue); instead, we 179 were interested in discrete relationships between SNPs and specific mRNAs and MEs. 180

Single mRNA expression associations 181
For hearts, the 390 single heart mRNAs tested had 79 significant independent genetic 182 associations among 56 unique mRNAs with 52 unique SNPs (Fig. 3). These 79 significant 183 associations with 56 mRNAs and 52 SNPs occur because a SNP tends to be associated with 184 expression of multiple mRNAs (average 1.58  0.95, max=5) and an mRNA tends to have more 185 than one eQTL (average 1.48  0.97, max=4). For brains, the 410 single brain mRNAs tested had 186 245 significant independent associations among 152 unique mRNAs with 146 unique SNPs (Fig.  187 3). Again, many SNPs were correlated to more than one mRNA (average 1.68  1.08, max=8), 188 and most mRNAs had more than one significant eQTL (average 2.94  2.49, max=14). Despite 189 testing a similar number of heart and brain mRNAs, brain had 3.1x more total significant eQTL 190 associations than heart (245 brain vs. 79 heart). There was also a greater proportion of mRNAs 191 with at least one significant eQTL in brain (152/410 [37.10%] compared to heart mRNAs 192 (56/390 [14.36%]) (Fig. 3). This difference is not explained by a difference in sample size (i.e., 193 power), which was similar between tissues. 194 One explanation for a single SNP being associated with multiple mRNAs is that the 195 mRNA affected by the SNP regulates the expression of many genes. This could occur when an 196 eQTL is found in a transcription factor protein or its regulatory region (e.g., promoter). To 197 determine if an eQTL impacted many mRNAs through a transcription factor, we annotated SNPs 198 and identified those within 5kb of a transcription factor (Table S4). Most of the changes were not 199 in transcription factors. Instead, of the seventeen heart eQTLs (33.7%) associated with more than 200 one mRNA (hereafter identified as "hotspots"), sixteen were in known protein coding regions 201 (within an intron or exon) but only four of these sixteen (25%) were within an annotated 202 transcription factor (erfl1, tada1, atf1, and zbtb3). The one eQTL hotspot not within a known 203 protein coding region was intergenic and not within 5kb of an annotated transcription factor. In 204 brains, 63 eQTLs (43.2%) were identified as hotspots, and of those, 32 were in known protein 205 coding regions (within an intron or exon); and only 2 of the 32 (6.25%) were within an annotated 206 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. ; https://doi.org/10.1101/2023.01.15.524165 doi: bioRxiv preprint transcription factor (erfl1 and zbtb3) while 1 intergenic SNP was also within 5kb of a zinc finger 207 protein (oocyte zinc finger protein XlCOF6-like). 208

Multivariate mRNA expression associations 209
Among the 39 heart co-expression modules, there were a total of 189 significant 210 associations between 102 eQTLME and 33 MEs (Fig. 4). These MEs are independent (not 211 correlated to each other) and do not share any mRNAs. As with the single mRNA associations, 212 the relationship between SNPs and MEs was not one to one. A given SNP was correlated with up 213 to six heart MEs with an average of 1.85  1.18 associations per SNP. Of the 102 unique SNPs, 214 54 were associated with a single heart ME, and the other 48 associated with two or more heart 215 MEs. Among those SNPs correlated with two or more heart MEs, the average number of MEs 216 correlations per SNP was 2.81  1.10. 217 Among the 41 brain co-expression modules, we found a total of 218 significant 218 associations between 133 eQTLME and 37 MEs. Single SNPs were correlated with up to six brain 219 MEs with an average of 1.64  1.06 associations per SNP (Fig. 4). Of the 133 unique SNPs, 87 220 were associated with a single brain ME, and the other 46 associated with two or more brain MEs. 221 Among those SNPs correlated with two or more brain MEs, the average number of ME 222 correlations per SNP was 2.85  1.01. Sample size was similar between heart and brain ME 223 associations. 224 Similar to the single mRNA eQTL results, several eQTLME were identified as potential 225 eQTLME hotspots because they were associated with expression of two or more MEs (Table S5). 226 Out of 48 heart eQTLME hotspots, 26 were in known protein coding regions (within an intron or 227 exon) with two in annotated transcription factors (jdp2 and atf1). None of the 22 intergenic 228 eQTLME heart hotspots were within 5kb of an annotated transcription factor. Out of 46 eQTLME 229 brain hotspots, 21 were in known protein coding regions (within an intron or exon) with two in 230 transcription factors (brf1b, trrap). Additionally, four intergenic eQTLME were within 5kb of an 231 annotated transcription factor (npas3, XlCGF7.1, and two SNPs in zmym3). One heart and one 232 brain hotspot were also within a gene involved in chromatin remodeling (chd3). 233 The eQTL studies for single and multivariate mRNAs indicated significant heritable 234 variation underlying tissue specific mRNA expression. Previously, multivariate mRNAs were 235 correlated to physiological traits (Drown et al. 2022). Here, we find that these multivariate 236 mRNA expression patterns are controlled by heritable variation. Thus, many eQTLME are 237 indirect drivers of physiological trait variation. Specifically, eQTLME underlying heart mRNA 238 expression are linked to 28°C CTmax, 12°C FA MRcardiac, 12°C SMR, 12°C LKA MRcardiac, and 239 12°C heart mass. Multivariate brain mRNA expression linked to 12°C body mass, 28°C CTmax, 240 and 28°C SMR is also under indirect genetic control via eQTLME. 241 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. ; https://doi.org/10.1101/2023.01.15.524165 doi: bioRxiv preprint

Trans-acting effects on mRNA expression 242
Although many (61.3% of eQTL and 50% of eQTLME) eQTL and eQTLME were within 243 genic regions (introns or exons), few (7.5% of eQTL and 6.4% of eQTLME) were found in or 244 near (within 5kb) transcription factors. This is similar to the percentage of transcription factors 245 found globally in eukaryotic genomes and indicates that there is not an enrichment for eQTL or 246 eQTLME within or near transcription factors (Riechmann 2002 To better understand the genomic context of eQTL and eQTLME we determined their 252 proximity to the mRNA(s) they effect. Interestingly, we found that >95% of single mRNA 253 eQTLs (97.4% for heart and 96.7% for brain) were trans-acting (defined as SNPs found on a 254 different chromosome or scaffold than that of the mRNA with which they were associated). 255 Specifically, for heart 94.1% (16/17) of eQTL hotspots, and 97.1% (34/35) of eQTL impacting 256 single mRNAs were trans-acting. For brain 95.2% (60/63) of eQTL hotspots, and 96.4% (80/83) 257 of eQTL impacting single mRNAs were trans-acting. While eQTLME could not be classified as 258 cis or trans-acting because they effect co-expression modules containing many disparately 259 located mRNAs, we determined whether the eQTLME were found within 5kb of mRNAs that 260 were part of the module. All heart and 98% of brain eQTLME were in genes not containing 261 mRNAs from that module. Thus, eQTLME are not simply cis-acting on a high-ranking mRNA 262 within the module but have broad effects on the expression of many module mRNAs. 263 264

Patterns of shared association 265
For both heart and brain, we looked for eQTL association for 10 mRNAs from each co-266 expression module. Within an ME, the 10 mRNAs have some degree of correlation with each 267 other causing them to be grouped into the same module. Thus, we examined whether mRNAs 268 from the same module had more shared eQTL than those from different modules. Surprisingly, 269 within 83% of heart and 67% of brain modules the top ten mRNAs all had unique eQTLs. This 270 "uniqueness" is demonstrated in the frequency of shared eQTLs within versus among modules: 271 for both hearts and brains, there was fewer shared eQTLs within modules than among modules 272 (T-test, heart p-value=0.027, brain p-value=0.008). This eQTL uniqueness would be expected if 273 they were cis acting, for example, a SNP in the mRNA promoter. Yet, many eQTLs are not near 274 or within the genes encoding the mRNA they effect. This indicates that these eQTL are trans-275 acting regulatory SNPs in the sense that they are found outside of the genes that code for 276 modular mRNAs. Similarly, for eQTLME, nearly all (98% of brain and 100% of heart) were in 277 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. ; https://doi.org/10.1101/2023.01.15.524165 doi: bioRxiv preprint genes not containing mRNAs from that module. Instead, they were found in different parts of the 278 genome and in genes unassociated with the physiological trait(s) correlated to the module. 279 Finally, we looked for overlap among SNP sets associated with a given phenotype 280 (physiological traits, tissue specific single mRNAs, or tissue specific MEs, Fig. 5). First, to 281 address whether the eQTLME were impacting module expression through action on a single high-282 ranking mRNA in that module, we assessed whether any SNP was associated with an mRNA and 283 the module containing that same mRNA (overlap in eQTL and eQTLME). There were no shared 284 mRNAs between eQTL and eQTLME, suggesting that the eQTLME are not acting on a single 285 mRNAs but represent a more complex mechanism of multivariate mRNA expression control. 286 Second, we examined overlap between tissues for SNPs associated with both single mRNAs and 287 modules. There was only one mRNA (atp7a, 0.13% of total) explained by at least one heart and 288 at least one brain eQTL. However, the eQTLs for this mRNA were not the same for heart and 289 brain expression. These data suggest that in the different tissues the variation in atp7a mRNA is 290 affected by different SNPs. There is a limited generality of this finding because we only tested 291 single heart and brain mRNAs that were high-ranking in modules, with only 3% of tested 292 mRNAs shared between heart and brain. Thus, while our results suggest that genetic control of 293 mRNA expression is tissue specific, examining a larger set of mRNAs expressed in several 294 tissues would likely be more informative about the role of conserved eQTL among tissues. 295

Genetic diversity 296
To compare genetic diversity among SNP sets, heterozygosity (He) was quantified. Among 297 all 1,406,282 high-probability variant sites, the average heterozygosity was 0.231  0.164 (mean 298  standard deviation). Heterozygosity for eQTL SNPs for heart and brain were significantly 299 higher than for all SNPs (heart = 0.328  0.074, brain = 0.325  0.068) and was highest for 300 eQTLME SNPs (heart = 0.368  0.091, brain = 0.361  0.099). Tissue specific eQTL and eQTLME 301 SNPs did not significantly differ in heterozygosity (Fig. S3). 302 303 We examined associations among SNPs and specific mRNAs that were previously 304 identified as biologically important based on their membership in co-expression modules. This 305 reduced our test from ~1.5 million SNPs across ~10,000 mRNAs in each tissue to <500 mRNAs 306 in each tissue. Additionally, we use WGCNA (Langfelder and Horvath 2008) to summarize 307 multivariate mRNA expression for co-expressed mRNAs. Using this approach, we summarized 308 ~10,000 mRNAs per tissue into 39 heart and 41 brain modules that could be used in our 309 association tests. In comparison to testing each of these mRNAs individually this may have 310 reduced the signal to noise ratio in our mRNA expression data and increased our power by 311 reducing the number of tests (Westra and Franke 2014). 312 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

Discussion
The copyright holder for this preprint this version posted January 18, 2023.
The data suggest that there are many more significant indirect eQTL than direct SNP 313 associations for the physiological traits we examined. Most of these indirect eQTL were trans-314 acting and in known protein coding regions (within introns or exons). Many (33.7% of heart and 315 43.2% of brain) eQTL were hotspots--associated with more than one mRNA. These results are 2021), we acknowledge that this may limit our findings in two ways. First, small sample sizes 321 could limit findings to only large effect loci. Yet, the important insights are that nearly all 322 significant eQTL are trans-acting and affect multiple mRNAs that are linked to physiological 323 variation. The observation that nearly all significant eQTL are trans-acting suggests that trans-324 acting eQTL have larger effect size than cis acting eQTL, which allowed us to detect them here. 325 Second, small sample sizes may increase the risk of detecting spurious associations. Yet, the 326 observation that eQTL were trans-acting and affecting multiple mRNA indicates that the eQTL 327 are not spurious. This is, both because it is unlikely that multiple independent associations would 328 spuriously be significant for a given eQTL and because it is unlikely that trans-acting factors 329 (defined as being on one of the 24 separate chromosomes) would represent more than 90% of 330 significant associations if they were spurious. Thus, while all our conclusions are based on the 331 limits of detections, these limits indicate the importance of trans-acting factors affecting 332 physiologically important mRNA expression. Still, we caution against making assumptions about 333 the role of specific SNPs that have been identified here in explaining physiological traits or gene 334 expression in other populations. This is both due to the limited sample size and the lack of direct 335 evidence (e.g., a traditional QTL study) for potential quantitative trait loci and eQTLs. Instead, 336 we focus the discussion on general patterns that speak to the genetic and molecular control of 337 physiological traits. 338 Our first question asked if there were SNPs directly associated with physiological trait 339 variation. We found only five SNPs directly association with physiological trait and all five were 340 to substrate specific MRcardiac. Second, we asked if mRNA expression patterns were under 341 genetic control. We found many significant eQTL associated with single mRNAs or MEs in 342 heart or brain tissue suggesting that mRNA expression patterns are influenced by genotype and 343 likely heritable (Table S2). Third, we leveraged data from our prior study that found significant 344 correlations between mRNA expression modules and physiological traits to ask if genetic control 345 of mRNA expression was likely to impact physiological trait variation. Nine heart modules and 346 four brain modules were previously identified as biologically important (correlated to 347 physiological traits). Of these, seven out of nine heart and all four brain modules were associated 348 with at least one eQTLME. Our results suggest that physiological traits are heritable via genetic 349 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made Only one direct associated SNP was also an eQTL (erfl1), the remaining eQTLs suggest unique 354 genetic and molecular mechanisms underlying the physiological traits we examined. 355 We treated individuals as belonging to a single population based on the admixture analysis, 356 which revealed no structure among the three sampled populations (Fig S1) In both heart and brain tissue, we found eQTLs for 14.6% and 20.7% of the tested heart 405 and brain mRNAs respectively, suggesting that heritable genetic variation is important in 406 explaining these mRNA expression patterns (Fig. 2). Interestingly, the majority (97.4% for heart 407 and 96.7% for brain) of eQTL that we detected were found on different chromosomes than the 408 affected mRNA (trans-acting). This could reflect a higher effect size or greater importance of 409 trans versus cis acting eQTLs for the selected mRNAs expression. While we did not classify SNPs associated with MEs as cis or trans because there is no 430 single genomic location for a module, we did examine whether SNPs associated with modules 431 were found in genes that produced mRNAs that were part of that module. We found no 432 association between SNP identity and module membership. That is, nearly all the eQTLs for 433 MEs were not found in or within 5,000bp of the genes for the mRNAs in our MEs. Additionally, 434 the SNPs associated with single mRNAs were not the same as those associated with ME 435 expression, suggesting that eQTLME were functionally independent from any single mRNA. This 436 demonstrates that ME expression is not driven by the effect of an eQTL on a single high-ranking 437 mRNA within the module. No studies to our knowledge have used multivariate mRNAs defined 438 by a WGCNA analysis as a phenotypic outcome in a genetic association study. 439 The approximately 25 individuals for any one SNP limits conclusion on the diversity of 440 mRNA and molecular mechanisms affecting the six physiological traits. Yet, we find many 441 significant eQTLs for both single and multivariate mRNA expression where a vast majority were 442 trans-acting or distant to the variable mRNA loci. These findings suggest that our approach may 443 allow for better detection of trans-acting elements, which may impact expression of dozens or 444 hundreds of co-expressed mRNAs. 445 Finally, we find that eQTL and eQTLME have greater heterozygosity when compared to all 446 1,406,282 variant sites used in this study. This may be explained by the association approach as 447 less variable sites have less variance among individuals that can be used to explain variance in 448 the physiological traits and mRNA expression patterns. It is also possible that heterozygosity 449 differences are biologically relevant. If eQTL and eQTLME sites are under directional selection, 450 we might expect a loss of genetic diversity in these sites. Yet, we find that heterozygosity is 451 higher in eQTL and eQTLME when compared to all SNPs. Higher heterozygosity may be due to 452 genetic redundancy, as indicated by the presence of different SNPs underlying eQTL and 453 eQTLME and the association of correlated physiological traits with different underlying loci. 454 These patterns could be driven by spatial or temporal variation in selection preventing allelic 455 fixation and aiding in the maintenance of biologically important variation at the SNP and mRNA 456 level. 457 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made with 12°C SMR and that these traits were both associated with two MEs (heart ME4 and heart 506 ME5). If these traits also had a shared genetic basis, we expected eQTLME for heart ME4 and 507 heart ME5 to be shared. Despite the shared association with these two traits, these MEs were 508 associated with two different eQTLME each. Similar patterns were found in brain among ME2, 509 ME3, and ME4, which were all correlated with 28°C CTmax but did not share any eQTLME. In 510 general, eQTLME were correlated to more than one ME (46.2% of heart and 35% of brain); 511 however, only four SNPs (1.7%) were shared among >1 module that was also associated with a 512 physiological trait (Fig. S4). Additionally, the five SNPs associated directly with physiological 513 traits were not shared among traits and do not overlap with ane eQTL, although only one gene 514 contains multiple SNPs that are either directly associated with a trait or are eQTLs (erfl1). 515 Overall, the limited overlap among SNPs associated directly or indirectly (eQTL) with 516 physiological traits within a tissue was surprising. This may be due to our limited ability to 517 detect small effect QTL, and we might expect greater overlap between QTL and eQTL if more 518 mRNAs were tested across more individuals. However, within the power of our data, we 519 detected diverse and complex molecular mechanisms contributing to physiological trait 520 variation. Tissue specific expression patterns appear to be under unique genetic control and 521 multivariate mRNA expression is not explained by a single eQTL impacting mRNA expression 522 of a gene highly correlation to a given module. 523 Thus, while different traits are correlated to the same ME, the nucleotide polymorphism, or 524 genetic control, of mRNA expression is distinct. This suggests that there is substantial genetic 525 variation underlying the physiological traits we have measured, with a diversity of molecular and 526 genetic mechanisms contributing to trait variation. The paradoxical genetic independence of 527 physiologically related traits (here metabolism and thermal tolerance) is not uncommon (see 528 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. ; https://doi.org/10.1101/2023.01.15.524165 doi: bioRxiv preprint (Herrewege 1980;Baker et al. 2015;Healy et al. 2018)) and emphasizes that these traits may still 529 be evolutionarily distinct, although they are linked at the molecular or physiological level. 530

531
Relationships among genotype, gene expression, and physiological traits explain 532 biologically important natural variation found in wild populations. In particular, substantial and 533 diverse genetic variation impacts these traits through direct and indirect (eQTL and eQTLME) 534 mechanisms. Demonstrated here, much of the mRNA expression variation is associated with a 535 diverse set of trans-acting eQTLs. Surprisingly, these trans-acting eQTLs are unique even for 536 mRNAs that affect multiple traits. Under a simpler genetic architecture, we may expect mRNAs 537 that have a shared association with cardiac and whole animal metabolism to also share the same 538 trans-acting eQTL, but this does not occur. Instead, the mRNA expression changes that affect 539 multiple physiological traits are regulated by different trans-acting SNPs. Finally, the SNPs 540 directly or indirectly associated with physiological traits have greater heterozygosity (genetic 541 variation) compared to all SNPs, and this greater genetic variation likely contributes to F. GuHCl buffer. DNA was extracted using carboxyl coated magnetic beads. The DNA quality was 559 assayed using gel electrophoresis and spectrophotometry to ensure high molecular weight and 560 low contamination. 561 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. ; https://doi.org/10.1101/2023.01.15.524165 doi: bioRxiv preprint

Library preparation 562
The analysis presented here uses a subset of samples that were part of a larger sequencing 563 run. A total of 1121 individuals were sequenced (Table S3). All samples were quantified in 564 triplicate using spectrophotometry and normalized to 100ng for sequencing library preparation. 565 The whole genome sequencing library was prepared using a tagmentation approach. Briefly, 566 DNA was digested with an in-house purified Tn5 transposase (as in (Picelli et al. 2014)) loaded 567 with partial adapter sequences. After tagmentation, the fragmented DNA was amplified using 568 barcoded primers such that each individual sample would contain a unique i7 and a plate level for F18 samples was added to the pool to achieve higher coverage because whole animal, tissue, 575 and molecular (mRNA expression) level phenotypic data was available for these individuals. 576

Raw sequence analysis 577
We followed best practices for low coverage whole genome sequencing (lc-WGS) data 578 processing as in (Lou et al. 2021). Briefly, adapter sequences and low-quality bases were 579 trimmed using Trimmomatic (v0.39) (Bolger et al. 2014). Flash (v1.2.11)(Magoč and Salzberg 580 2011) was used to combine overlapping reads and to parse singletons and paired reads. 581 Singletons and paired reads were mapped separately using BWA mem (v0.7.17) and resulting 582 sam files converted to bam files using samtools (v1.3.1) (Danecek et al. 2021). The first and 583 second sequencing run were processed separately until BAM files were produced and found to 584 be of similar quality assessed by comparing total percentage of mapped reads and levels of 585 dually mapped reads before combining for the remaining file processing steps. Picard (v2.26.4) 586 was used to add read group information, which is needed for duplicate marking downstream. 587 After combining all mapped reads for a single individual, BAM files were further filtered for 588 mapping quality using samtools and overlapping reads softclipped using bamutil (v1.0.15). Heart and brain tissues were collected after measuring MRcardiac and stored in chaotropic 610 buffer for mRNA expression analysis. The mRNA data includes tissue specific expression counts 611 for single mRNAs and tissue specific module mRNA expression (ME) from a whole genome co- (ME) as a multivariate molecular level phenotype that may be predicted using genotype 618 likelihoods. In addition, we examined association between genotype likelihoods and the top 10 619 mRNAs for each module (based on module membership). 620

Association studies 621
All results are reported from the score test conducted in ANGSD using -doAsso 2 with 622 filters -minHigh 4, -minCount 1. These filters are relaxed from the default filters (default: -623 minHigh 10, -minCount 10) and allowed us to examine associations for all traits with reduced 624 sample size requirements. The sample size for each association can be found in Table S2. For all 625 phenotypes, acclimation temperature was included as a covariate. For SMR, CTmax, and 626 MRcardiac, acclimation order (individuals were acclimated to 12°C then 28°C or 28°C then 12°C) 627 was included as an additional covariate. For SNP associations to mRNA expression, heart and 628 brain mRNA expression were examined as separate phenotypes. P-values for genotype to 629 phenotype associations were corrected for multiple testing using the Benjamini-Hochberg 630 approach (Benjamini and Hochberg 1995) and significant associations identified as those with a 631 corrected p-value <0.05. To examine patterns among independent SNPs, in cases where SNPs 632 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  Physiological trait variation can be driven directly or indirectly (through gene expression) by genotype. To understand the molecular and genetic basis of physiological trait variation, comprehensive datasets can be used to investigate 1) direct associations between genotype and physiological traits, 2) direct correlations between gene expression and physiological traits, for example using weighted gene co-expression analysis (WGCNA), and 3) indirect effects of genotype on physiological traits, which may occur when expression quantitative trait loci (eQTL) impact expression of genes underlying physiological traits.
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023.  . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. ; https://doi.org/10.1101/2023.01.15.524165 doi: bioRxiv preprint 668 669 670 671 Figure 4: Distribution of multivariate mRNA co-expression associated eQTLs among modules that explain physiological traits. Number of significant SNPs (y-axis, count) associated with each mRNA coexpression module in A) heart and B) brain tissue. Modules are tissue specific and were identified using a weighted gene co-expression network analysis. Modules previously found to be associated with physiological traits have bright colors and are labeled. Most modules were associated with >1 SNP with an average of 5.72  3.84 SNPs per module for heart (max = 13) and 5.74  4.35 SNPs per module for brain (max = 18).
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. ; https://doi.org/10.1101/2023.01.15.524165 doi: bioRxiv preprint 672 673 674 Figure 5: Unique eQTLs explain single mRNA and multivariate mRNA co-expression between heart and brain tissues. Overlap among eQTLs associated with single mRNAs and co-expressed mRNA modules (Heart and Brain Modules) and SNPs associated with physiological traits. In both tissues, SNPs uniquely explain either single or co-expressed mRNAs with few shared among SNP sets (five rightmost vertical bars). Twelve (6.2%) eQTLs are correlated with single mRNA expression in both heart and brain tissues and seven (3.0%) eQTLs are correlated with module expression in both heart and brain tissues. Finally, a single SNP (chr NC046363.1 location 31037486) was correlated to expression of lmnb1 (a protein involved in cellular architecture) in brain and endogenous cardiac metabolism. Horizontal bars (far right) show the number of SNPs associated with each variable (single mRNA expression, module mRNA expression, or physiological traits), vertical bars (top) show number of SNPs unique (five leftmost vertical bars) or shared (five rightmost vertical bars) among datasets.
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. ; https://doi.org/10.1101/2023.01.15.524165 doi: bioRxiv preprint Figure S3: Heterozygosity among all single nucleotide polymorphisms and eQTL sets. Heterozygosity was significantly higher for SNPs identified as eQTLME and eQTLs when compared to all SNPs. Heterozygosity was significantly higher for eQTLME compared to eQTLs with no difference between tissues.
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. ; https://doi.org/10.1101/2023.01.15.524165 doi: bioRxiv preprint 681 682 Figure S4: Overlap among eQTLME for modules correlated to physiological traits. Genetic control of modules that explain physiological trait variation appears independent with few shared SNPs in any eQTLME set. The exceptions are a single shared eQTLME between brain ME1 and brain ME4 (no shared trait correlation), between brain ME1 and brain ME2 (no shared trait correlation), and among heart ME3, heart ME4, and heart ME8 (shared correlation of ME3 and ME4 with 12°C FA MRcardiac and shared correlation of ME4 and ME8 with 12°C SMR).
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. ; https://doi.org/10.1101/2023.01.15.524165 doi: bioRxiv preprint . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted January 18, 2023. ; https://doi.org/10.1101/2023.01.15.524165 doi: bioRxiv preprint