# The adaptive potential of subtropical rainbowfish in the face of climate change: heritability and heritable plasticity for the expression of candidate genes

R. J. Scott Mccairns, Steve Smith, Minami Sasaki, Louis A Bernatchez, Luciano B. Beheregaray

<figure><img src="image_1.png" /></figure>

## To cite this version:

R. J. Scott Mccairns, Steve Smith, Minami Sasaki, Louis A Bernatchez, Luciano B. Beheregaray. The adaptive potential of subtropical rainbowfish in the face of climate change: heritability and heritable plasticity for the expression of candidate genes. Evolutionary Applications, 2016, 9 (4), pp.531–545. 10.1111/eva.12363 . hal-01444059

HAL Id: hal-01444059

https://hal.science/hal-01444059

Submitted on 27 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ORIGINAL ARTICLE

# The adaptive potential of subtropical rainbowfish in the face of climate change: heritability and heritable plasticity for the expression of candidate genes

R. J. Scott McCairns,<sup>1,2</sup> Steve Smith,<sup>1,3</sup> Minami Sasaki,<sup>1</sup> Louis Bernatchez<sup>4</sup> and Luciano B. Beheregaray<sup>1</sup>

1. Molecular Ecology Laboratory, Flinders University, Adelaide, SA, Australia

2. Centre National de la Recherche Scientifique, UMR 7263 (IMBE) Institut Méditerranéen de la Biodiversité et d'Ecologie Marine et Continentale, Marseille, France

3. Department of Integrative Biology and Evolution, University of Veterinary Medicine, Vienna, Austria

4. Institut de Biologie Intégrative et des Systèmes, Université Laval, Québec, QC, Canada

**Keywords**

adaptation, Australia, climate change, ecological genomics, evolutionary physiology, gene expression, teleost, thermal tolerance.

**Correspondence**

Luciano B. Beheregaray, School of Biological Sciences, Flinders University, Adelaide, SA 5001, Australia.

Tel.: +61 8 82015243;

Fax: +61 8 82013015;

e-mail: luciano.beheregaray@flinders.edu.au

Received: 3 September 2015

Accepted: 14 January 2016

doi:10.1111/eva.12363

**Abstract**

Whilst adaptation and phenotypic plasticity might buffer species against habitat degradation associated with global climate change, few studies making such claims also possess the necessary and sufficient data to support them. Doing so requires demonstration of heritable variation in traits affecting fitness under new environmental conditions. We address this issue using an emerging aquatic system to study adaptation to climate change, the crimson-spotted rainbowfish (*Melanotaenia duboulayi*), a freshwater species from a region of eastern Australia projected to be affected by marked temperature increases. Captive born *M. duboulayi* of known pedigree were used to assess the long-term effects of contemporary and 2070-projected summer temperatures on the expression of genes previously identified in a climate change transcriptomics (RNA-Seq) experiment. Nearly all genes responded to increasing temperature. Significant additive genetic variance explained a moderate proportion of transcriptional variation for all genes. Most genes also showed broad-sense genetic variation in transcriptional plasticity. Additionally, molecular pathways of candidate genes co-occur with genes inferred to be under climate-mediated selection in wild *M. duboulayi* populations. Together, these results indicate the presence of existing variation in important physiological traits, and the potential for adaptive responses to a changing thermal environment.

**Introduction**

Global climate change poses many threats to biodiversity, but for ectotherms, the associated impacts affecting temperature niche may be particularly problematic (Pearson et al. 2014). Given pronounced intracontinental temperature increases (Diffenbaugh and Field 2013), and the influence of localized climatic effects on their thermal properties (Caissie 2006), inland and coastal freshwaters may be especially vulnerable to temperature increases. Indeed, freshwater habitats and their biota are among the most threatened by environmental change (Dudgeon et al. 2006), a concern highlighted in the most recent Intergov-

ernmental Panel on Climate Change (IPCC) report (Settele et al. 2014). All freshwater fishes are ectotherms, and thus, metabolically sensitive to environmental temperature. How they are likely to respond to the challenges associated with global climate change is in many ways dependent upon the effects of new environmental inputs on essential physiological processes (Helmuth 2009; Chown et al. 2010; Hofmann and Todgham 2010). Within the context of population-level responses to global climate change, adaptation and phenotypic plasticity are often cited as potential means of species persistence. Yet, few studies present empirical data essential to understanding whether either of these represents viable mechanisms. Critically absent from
the discussion are demonstrations of a genetic basis for
traits likely to influence fitness in the changed environ-
ments (Gienapp et al. 2008; Charmantier and Gienapp
2014; Merilä and Hendry 2014). Moreover, despite a per-
ceived ubiquity of biological adaptation in nature, evidence
of the potential to genetically adapt to climate change is
scarce and ambiguous (Bradshaw and Holzapfel 2006;
Hoffmann and Sgrò 2011; Merilä and Hendry 2014). This
is particularly so in aquatic organisms from the southern
hemisphere, which are largely under-represented in climate
change research (Chambers et al. 2013).

Acclimatization – acclimation within the context of a
controlled, laboratory setting – is a compensatory change
in physiological processes in response to environmental
change (Willmer et al. 2005) and represents a form of phe-
notypic plasticity. Transcriptional mechanisms appear to
play a key role in acclimation to environmentally induced
stressors across broad taxonomic groups (López-Maury
et al. 2008), and in the regulation of metabolism (Des-
vergne et al. 2006). For ectothermic species such as fishes,
the latter is also highly relevant with respect to changes in
the thermal environment. In many fish species, acclimation
can increase threshold temperatures for induction of the
cellular stress response and upper thermal toler-
ance (Fangue et al. 2011; Currie and Schulte 2014). How-
ever, this may come at the cost of increased resting
metabolic rate and related long-term fitness conse-
quences (Pörtner et al. 2006; Pörtner and Knust 2007;
Muñoz et al. 2012). Conversely, interindividual variation
in the thermal response may be reflective of genetic vari-
ability in the underlying physiological processes. Yet, whilst
a genetic basis to functional biology is frequently
assumed in the literature, estimates of heritability for varia-
tion in physiological processes are scarce, especially in
fishes (Chown et al. 2010). This is particularly problematic
with respect to forecasting the effects of climate change, as
any adaptive response is dependent upon the degree of
additive genetic variation underlying organisms’ thermal
physiology. If we are to address these issues in a rigorous
manner, we require ecologically relevant models that are
also amenable to multigenerational experimentation.

Rainbowfishes (Melanotaeniidae) are a diverse (ca. 81
species) and monophyletic group endemic to Australia and
New Guinea (Unmack et al. 2013). They inhabit a diverse
range of freshwater habitats, being found almost every-
where that freshwater is present in Australia (Allen 1989;
McGuigan et al. 2000; Unmack et al. 2013). Despite their
broad geographic distribution and abundance, they are
nonmigratory, typically with small home ranges and strong
microhabitat preferences (Hurwood and Hughes 2001;
Hattori and Warburton 2003). In addition, many rainbow-
fish taxa have radiated relatively recently (McGuigan et al.
2000), especially the lineages endemic to Australia

(Unmack et al. 2013). Given their close evolutionary rela-
tionships, strong metapopulation structures and diverse
habitat preferences, these fishes represent a powerful field
model for the study of local adaptation. Most species also
display key attributes of a valuable laboratory model. They
are small (ca. 3–9 cm), sexually dimorphic, easy to main-
tain and breed in captivity (they are very popular aquarium
fishes), are highly fecund (between 150 and 1600 eggs) and
have a short generation time, reproducing annually in the
wild and at the age of 3 months in the laboratory (Pusey
et al. 2001; Holdway et al. 2008).

The crimson-spotted rainbowfish (*Melanotaenia dubou-
layi*) is an Australian subtropical species found in vegetated
streams over a relatively broad latitudinal range in coastal
drainages east of the Great Dividing Range (ca. 31.4°S to
21.4°S; Allen 1989). Temperatures along this distribution
have been rising steadily, and high CO₂ emissions projec-
tions may produce temperature as high as 10°C warmer by
2070 in the worst affected waterbodies (Meehl et al. 2007;
Reisinger et al. 2014). Contemporary mean summer tem-
peratures (ca. 21°C) are below the preferred temperature
(26°C) for *M. duboulayi* tested in a thermal gradient (King
and Warburton 2007). However, projected 2070 tempera-
tures are expected to exceed the species preference and may
approach levels near its critical thermal temperature (Patra
et al. 2007; Beheregaray, unpublished data). Additionally,
wild-caught *M. duboulayi* exhibit substantial plastic regula-
tory responses to temperature stress. In a short-term tran-
scriptomics (RNA-Seq) experiment, Smith et al. (2013)
identified and annotated 614 upregulated and 349 down-
regulated genes in relation to predicted summer tempera-
tures for 2070. These include a suite of heat-shock genes
that responded sharply to increased temperatures, and
genes related to regulation of metabolic functions and
developmental processes that showed mid-range changes in
expression (Smith et al. 2013). Thus, the availability of
excellent candidate genes for climatic adaptation, coupled
with the species’ life history attributes, makes *M. duboulayi*
an ideal aquatic study system to assess molecular footprints
of adaptive resilience to climate change.

In this study, we assess the adaptive potential to climate
change of the subtropical Australian rainbowfish
*M. duboulayi*, focusing on a genetically healthy population
that is centrally located along the climatically heteroge-
neous latitudinal range of the species. Specifically, we use
the *M. duboulayi* system to examine the effects of long-
term acclimation to increased temperature on the tran-
scription of candidate genes which are differentially
expressed in response to thermal stress, including a test for
sex-specific differences in this response. We hypothesize
that there is a heritable basis underlying the variable tran-
scriptional responses revealed during earlier experiments
that used short-term exposure to temperature treatments
mimicking contemporary and future climate. Critically, we employ a reciprocal common-garden design and incorporate information on the pedigree of experimental animals, to estimate both narrow-sense heritability of transcription and broad-sense genetic variation for its plasticity. Additionally, we test for functional associations with candidate genes inferred to be under temperature-mediated selection in the wild (identified in a companion study) as a preliminary exploration of the potential adaptive value of transcriptional variation.

## Materials and methods

### Origin of broodstock, spawning & husbandry

Sixty adult *M. duboulayi* were collected using bait traps from a location in the upper reaches of the Brisbane River near the township of Fernvale (27°26′37.39″S, 152°40′12.76″E). The Brisbane River population was targeted because it is abundant compared with several *M. duboulayi* found in smaller rivers and is centrally located along the north-south range of the species in eastern Australia. In addition, a landscape genomic study of *M. duboulayi* based on over 17 000 SNPs screened in 22 populations from across the species range (Smith et al. unpublished) indicates that Brisbane River rainbowfish are part of a large representative group of populations locally adapted to climate-related variables (hottest maximum and coldest minimum temperatures, run-off, quarterly rainfall and annual radiation).

Fish were transported by air to the Animal House and Aquaculture facility at Flinders University in Adelaide and allowed to acclimate for 6 months. During the acclimation period, fish were maintained in single sex tanks (12 fish/100L) at 21°C under conditions of 12 h light/12 h dark. They were fed once daily on a mixture of blood worms and fish pellets. After 6 months, photoperiod was changed to a cycle of 14 h light/10 h dark to induce breeding behaviour. When fish began exhibiting nuptial coloration and gravidity, 11 mating groups of a single male and two females were selected at random and transferred to breeding tanks. Breeding tanks were 100L volume and contained 'spawning mops' as artificial substrate for egg deposition. Mops were removed daily in the morning and were transferred to corresponding rearing tanks. Breeding was conducted over a 4-month period from June to October 2012. Larvae were then reared in separate, 100-L family-rearing tanks at 21°C under a photoperiod of 12 h light/12 h dark and fed daily, initially with rotifers and later with blood worms and fish pellets.

### Pedigree reconstruction

DNA was extracted from all parents and offspring (n = 199) using a modified salting-out procedure (Sun-

nucks and Hales 1996). DNA quality and quantity were measured with a NanoDrop 1000 spectrophotometer. All individuals were genotyped using microsatellite markers developed for *M. australis* (Ma03, Ma09, Ma10 & Ma11; Young et al. 2009) and *M. splendida* (Ms22, Ms24 & Ms40; Zhu et al. 1998). Loci were amplified in 5 µL polymerase chain reactions (PCRs) using procedures and touch-down cycles following (Beheregaray and Sunnucks 2000). Microsatellite fragments were scored using GENEMAPPER 4.0 (Applied Biosystems, Foster City, CA, USA), and scoring error evaluated using MICRO-CHECKER (van Oosterhout et al. 2004). As breeding tanks contained only a single male, offspring paternity was known; maternity was inferred from resulting genotype of offspring and parents using COLONY (version 2.0.5.8; Jones and Wang 2010). Individuals were separated into full-sib groups and kept in separate tanks until thermal treatments began.

### Thermal treatment & tissue sampling

At approximately 6 months of age, individuals from each family were assigned randomly, and in equal proportions, to one of two groups: a control representing current mean summer temperatures (21°C), and the projected mean 2070 summer temperature (31°C) for the Brisbane River region (Meehl et al. 2007; Reisinger et al. 2014). It should be noted that this projected summer temperature is below the upper thermal critical maximum (CT<sub>MAX</sub> = 38°C; SD = 0.39) of *M. duboulayi* (Beheregaray, unpublished data), empirically estimated for Brisbane River rainbowfish following the method of Becker and Genoway (1979). Individuals from each temperature group were also randomly assigned to one of 19 replicate tanks. For the 31°C treatment, temperature was increased by 2°C per day until the experimental temperature was achieved. Tanks were monitored daily for signs of stress or disease associated with the higher temperature. Fish were then maintained for 80 days to approximate the typical duration of summer temperatures in the Brisbane River region.

At the end of the treatment period, fish were sacrificed using AQUI-S® solution and dissected to remove their livers for immediate RNA extraction. Although increased temperature has been shown to differentially induce transcription changes in different tissue types, we selected liver due to previous research specifically linking this tissue type to heat stress responses (Smith et al. 2013). The final data set contained 108 individuals (60 females & 48 males), assigned in approximately equal proportions to control (21°C: 26 females; 21 males) and treatment (31°C: 34 females; 27 males) groups. Nineteen families were present, 16 of which were related as paternal half-sibs, spawned from 19 dams and 11 sires. Family size ranged from one to 16 individuals (median = 5). All families comprised of six
or fewer individuals had multiple, related individuals in a
corresponding half-sib family; the three families with no
half-sib relations contained from six to 14 individuals.

RNA extractions & qRT-PCR

Total RNA was extracted from liver tissue using the Ambion
Magmax™-96 total RNA isolation kit (Thermo Fisher Scien-
tific, Waltham, MA, USA), following Smith et al. (2013).
RNA quality and concentration were measured using a Nan-
drop 1000 spectrophotometer, then diluted to a standard-
ized concentration (10 ng/µL). Each dilution was checked
for DNA contamination by PCR amplification of a house-
keeping gene using the SensiFAST SYBR Hi-ROX kit (Bio-
line, London, UK), followed by electrophoresis on agarose
gel. Samples exhibiting amplified DNA were treated with
DNase MAX Kit (MO BIO) for a second time. Finally, we
used random primers and the High Capacity cDNA Reverse
Transcription Kit (Applied Biosystems) for cDNA synthesis,
following manufacturer's protocols.

Primers were designed to amplify twelve candidate genes that are putatively differentially expressed in response to thermal stress (Smith et al. 2013) – primer sequences and related information can be found in Table S1. Most importantly, candidate genes are believed to play a role in biological processes associated with a response to thermal challenge (e.g. oxidation reduction, metabolism, regulation of transcription; see Table S2). Candidates included homologues for the following zebrafish (Danio rerio) genes: alanine: glyoxylate aminotransferase a (agxta); cytochrome P450, family 1, subfamily A (cyp1a); glutathione S-transferase rho (gstr, a.k.a. zgc:162356); 3-hydroxy-3-methylglutaryl-coenzyme A synthase 1 (hmgcs1); heat-shock protein 90, alpha (cytosolic), class A member 1, tandem duplicate 2 (hsp90aa1.2); nuclear receptor subfamily 1, group D, member 4b (nr1d4b); NAD(P)-dependent steroid, dehydrogenase-like (nsdhl); peroxisome proliferator-activated receptor alpha b (pparab); thyrotrophic embryonic factor a (tefa); thioredoxin-related transmembrane protein 2b (tmx2b); UDP glucuronosyltransferase 2 family, polypeptide A4 (ugt2a4); and ubiquinol-cytochrome c reductase core protein IIB (uqcrc2b). Reference genes, to control for differences in template quantity due to variation among extractions and/or dilutions, were identified from a previous study as nondifferentially expressed between treatments (Smith et al. 2013) and included ba1 globin (ba1, a.k.a. hbb) and phosphoglycerate mutase 1a (pgam1a).

Quantitative PCR was performed in 96-well plates on a
StepOnePlus Real-Time PCR System (Applied Biosystems),
using a SYBR Green reporter. Reactions contained 1 µL of
diluted cDNA (10 ng/µL), 5 µL of 2× SensiFAST SYBR
Hi-ROX (Bioline), 0.05 µM of each forward and reverse
primer and 3.8 µL of RNase-free water. Thermal cycles

consisted of an initial denaturation step (95°C; 2 min), and
40 cycles of denaturation (95°C; 3 s), annealing (60°C;
10 s) and extension (72°C; 10 s). This was followed by a
melt curve cycle of 95°C for 15 s, 60°C for 1 min and 95°C
for 15 s. A single gene was amplified for a given 96-well
plate that also comprised: three negative controls; a four-
fold serial dilution of a positive control, also in triplicate;
three technical replicates for 30 experimental animals.

Analyses

For each gene, raw fluorescence values of all individuals,
including all technical replicates and controls, were com-
piled into a single data file. Initial mRNA quantities
(N₀) were estimated using a regression-based procedure
implemented in LinRegPCR (Ramakers et al. 2003; Rui-
jter et al. 2013). Replicate N₀ values of each candidate
gene were standardized relative to the geometric mean of
N₀ values for all technical replicates of reference genes
for a given individual (Vandesompele et al. 2002). Stan-
dardized data were then screened for outliers: technical
replicates greater than two standard deviations from
their corresponding family-by-treatment mean were
flagged and removed – on average, outlier removal
resulted in the loss of a single data point (i.e. individual)
per gene in the final data set (range 0–3). Following
standardization and outlier removal, technical replicates
were averaged by individual.

To evaluate effects of temperature and/or sex on
candidate gene transcription, and to estimate additive
genetic variance (V_A) for transcription, we analysed
data under a pedigree-weighted, mixed effects model (i.e.
the animal model) implemented in the R package
‘MCMCglmm’ (Hadfield 2010). In order to evaluate the
significance of additive genetic variance for transcription,
we ran two nested models: a minimal/null model including
only random variation among rearing tanks, and a second
model which also included ‘animal’ variance (V_A); both
models included coefficients describing fixed effects of tem-
perature treatment, sex and their interaction. To ensure
convergence and sufficient mixing, models were run for
25 000 000 iterations comprising an initial burn-in of
20 000 000 iterations and a sampling chain of 5 000 000
iterations. To reduce autocorrelation of estimates, each
5000th point of the Markov chain was sampled, yielding a
total of 1000 posterior estimates of model parameters.
Model (i.e. V_A) significance was evaluated based on the
deviance information criterion (DIC); significance of fixed
effects was determined by profiling the posterior distribu-
tions of estimated coefficients.

As mRNA extraction necessitated lethal sampling, indi-
viduals could not be sampled in multiple environments,
thereby precluding the possibility of estimating additive
genetic variance for transcriptional plasticity. However, to
evaluate the potential for genetic variance in plasticity
($V_G \times E$), we ran two additional nested models: a model
including broad-sense genetic variance in mean transcrip-
tion (i.e. variation among families in the intercept coeffi-
cient; $V_G$), and a model including broad-sense (i.e. among
family) variation in the temperature (i.e. slope) coefficient.
Both models also included random variation among rear-
ing tanks and were fully parameterized for all fixed effects.
Models were run with a burn-in of 500 000 iterations, fol-
lowed by an additional 500 000 iterations from which each
500th point was sampled. Significance of $V_G \times E$ was deter-
mined by contrasting the DIC score against both the $V_G$
and previous ‘null’ models.

Inferring functional/adaptive value

As a tentative first inference into the potential adaptive
value of variation in the candidate genes, we explored
their functional relationships with a subset of genes pre-
viously identified as putatively under selection in a natu-
ral thermal gradient. In a companion paper, Smith et al.
(unpublished) explored genomewide SNP variation along
the latitudinal range of *M. duboulayi* and observed 1141
loci putatively under climate-mediated selection, as iden-
tified via outlier analysis using Bayesian (Foll and Gag-
giotti 2008). We first conducted a distant homology
BLAST scan (gap opening and extension penalties 5 and
2; match/mismatch scores of 2 and -3; minimum per-
centage sequence identity 75) of the 93-bp sequences of
these outliers against the zebrafish genome (GRCz10;
accessed July 28, 2015), reducing the data set to 411
annotated genes. We then used the BioMart package to
retrieve identifiers (e.g. ENTREZ/NCBI accessions, offi-
cial gene symbols) for use in subsequent analyses (Dur-
inck et al. 2005).

We queried the KEGG database to identify all known
molecular pathways in zebrafish comprising a co-occur-
rence of at least one of the 12 candidate genes assayed for
transcription and one of the annotated outlier
genes (Kanehisa and Goto 2000; Kanehisa et al. 2014).
Additionally, we performed a functional enrichment analy-
sis and clustering of the two gene lists based on biological
process and molecular function ontologies using DAVID
(Huang et al. 2007; Sherman et al. 2007). Finally, we iden-
tified potential protein-protein interactions between mem-
bers of each gene set by querying the STRING
database (v10; accessed August 3, 2015; Snel et al. 2000;
Szkarczyk et al. 2015). For each of the 12 candidate genes,
we identified all known and predicted protein-protein
interactions with the 'highest' confidence score (0.900).
From this list, we identified 23 from the 411 genes puta-
tively under selection and plotted their interaction network
with the 12 candidate genes.

Results

Effects of temperature and sex on transcription

With the exception of *tmx2b* (Table 1; Fig. 1A), all candidate genes exhibited a significant response to increasing temperature. Nearly all genes were significantly upregulated at 31°C (Fig. 1B–J), with the exception of *hsp90aa1.2*, for which female transcription was decreased (Fig. 1K), and *nr1d4b*, which was generally downregulated in response to increasing temperature (Fig. 1L). Five genes were also differentially expressed between sexes: *agxta* (Fig. 1F), *gstr* (Fig. 1G) and *ugt2a4* (Fig. 1I) were significantly higher in males than in female; in contrast, transcription of *tmx2b* (Fig. 1A) and *tefa* (Fig. 1H) was significantly greater in females than in males. Significant interaction effects also revealed environment dependent sex differences. Transcription of *cyp1a* was significantly higher in males at 21°C; however, this pattern was reversed at 31°C, where female levels were higher (Table 1; Fig. 1J). For *hsp90aa1.2* and *nr1d4b*, sex differences were observed only at 21°C (Table 1), with significantly greater levels of transcription for females (Fig. 1K) and males (Fig. 1L), respectively.

Heritability and G × E for transcriptional variation

DIC for models containing both random variation due to
tank effects ($V_{\text{tank}}$) and $V_A$ were lower than those of corre-
sponding ‘null’ models (Table 2), indicating significant
additive genetic variance for transcription of all candidate
genes. For most genes, point estimates of $V_{\text{tank}}$ were an
order of magnitude less than those of $V_A$ (Table 1). With
the exception of *agxta* and *hsp90aa1.2*, residual/error vari-
ance ($V_{\text{resid}}$) tended to be highest. Consequently, narrow-
sense heritability for transcription was moderate, being
0.25 on average (range = 0.11–0.54; Table 2).

Most candidate genes also exhibited evidence of signifi-
cant broad-sense genetic variation for transcriptional plas-
ticity (Fig. S1), based on DIC comparisons with ‘null’ and
$V_G$ models (Tables 2 & S3). $V_{G \times E}$ was deemed nonsignif-
icant for only *agxta* and *tmx2b* (Fig. S2) and was question-
able for *hsp90aa1.2* given the relatively minor difference in
DIC scores among models – point estimates of random
slope variance ($V_{\beta\text{Treat|Fam}}$) for these genes were also sub-
stantially lower than those of $V_{\text{resid}}$ (Table S3). Of the genes
with putatively significant $V_{G \times E}$, family-specific reaction
norms were largely in the same direction: none exhibited a
crossing pattern indicative of both up- and downregulation
in response to temperature (Fig. S1). Although based
strictly on mean estimates *hmgcs1* showed a trend in this
direction for three families (Fig. 1B), 95%PDIs clearly
overlapped, indicating no significant plasticity in transcrip-
tion associated with temperature (i.e. slope estimates not
exclusive of zero, Fig. S1E). This was the pattern observed
<table><caption>**Table 1.** Statistical evaluation of model coefficients describing effects of temperature treatment (Treat) and differences between sexes in the rainbowfish *Melanotaenia duboulayi*. Letters in parentheses beside gene symbols correspond to the respective panels of Fig. 1. The average fold-change in response to temperature treatment is also reported in parentheses. For genes with a significant sex-by-treatment interaction, results of reduced models contrasting sexes within each temperature and describing sex-specific temperature effects are also presented. Variance estimates are based on the fully parameterized model including all fixed effects, additive genetic variance (V<sub>A</sub>) and random variance attributable to tank effects (V<sub>tank</sub>). Estimates are based on the posterior mode (Post. Mode) and are bounded by 95% posterior density intervals (PDIs).</caption><thead><tr><th rowspan="2">Gene</th><th rowspan="2">Fixed effects<br>Model term</th><th rowspan="2">P-Value</th><th colspan="3">Variance estimates</th></tr><tr><th>Post. Mode</th><th>95% PDIs</th></tr></thead><tbody><tr><td rowspan="3">agxta (F)</td><td>Treat (1.24x)</td><td>0.001</td><td>V<sub>A</sub></td><td>0.472</td><td>0.078 0.825</td></tr><tr><td>Sex</td><td>0.002</td><td>V<sub>tank</sub></td><td>0.069</td><td>0.025 0.309</td></tr><tr><td>Sex x Treat</td><td>0.774</td><td>V<sub>resid</sub></td><td>0.196</td><td>0.045 0.505</td></tr><tr><td rowspan="5">cyp1a (J)</td><td>Treat (3.05x)</td><td>0.001</td><td>V<sub>A</sub></td><td>0.105</td><td>0.024 0.383</td></tr><tr><td>Sex</td><td>0.014</td><td>V<sub>tank</sub></td><td>0.084</td><td>0.032 0.295</td></tr><tr><td>Sex x Treat</td><td>0.001</td><td>V<sub>resid</sub></td><td>0.256</td><td>0.091 0.368</td></tr><tr><td>Interactions</td><td></td><td></td><td></td><td></td></tr><tr><td>21°C: M vs F</td><td>0.004</td><td></td><td></td><td></td></tr><tr><td>31°C: M vs F</td><td>0.002</td><td></td><td></td><td></td></tr><tr><td>F: 21°C vs 31°C</td><td>0.001</td><td></td><td></td><td></td></tr><tr><td>M: 21°C vs 31°C</td><td>0.001</td><td></td><td></td><td></td></tr><tr><td rowspan="3">gstr (G)</td><td>Treat (4.05x)</td><td>0.001</td><td>V<sub>A</sub></td><td>0.011</td><td>0.005 0.051</td></tr><tr><td>Sex</td><td>0.006</td><td>V<sub>tank</sub></td><td>0.009</td><td>0.003 0.030</td></tr><tr><td>Sex x Treat</td><td>0.502</td><td>V<sub>resid</sub></td><td>0.029</td><td>0.010 0.045</td></tr><tr><td rowspan="3">hmgcs1 (B)</td><td>Treat (1.92x)</td><td>0.001</td><td>V<sub>A</sub></td><td>0.153</td><td>0.026 0.503</td></tr><tr><td>Sex</td><td>0.554</td><td>V<sub>tank</sub></td><td>0.050</td><td>0.016 0.203</td></tr><tr><td>Sex x Treat</td><td>0.182</td><td>V<sub>resid</sub></td><td>0.255</td><td>0.058 0.410</td></tr><tr><td rowspan="5">hsp90aa1.2 (K)</td><td>Treat (n.a.)</td><td>0.001</td><td>V<sub>A</sub></td><td>1.0E-03</td><td>1.7E-04 0.002</td></tr><tr><td>Sex</td><td>0.001</td><td>V<sub>tank</sub></td><td>2.7E-04</td><td>5.7E-05 0.001</td></tr><tr><td>Sex x Treat</td><td>0.001</td><td>V<sub>resid</sub></td><td>3.4E-04</td><td>3.9E-05 0.001</td></tr><tr><td>Interactions</td><td></td><td></td><td></td><td></td></tr><tr><td>21°C: M vs F</td><td>0.001</td><td></td><td></td><td></td></tr><tr><td>31°C: M vs F</td><td>0.074</td><td></td><td></td><td></td></tr><tr><td>F: 21°C vs 31°C</td><td>0.001</td><td></td><td></td><td></td></tr><tr><td>M: 21°C vs 31°C</td><td>0.448</td><td></td><td></td><td></td></tr><tr><td rowspan="5">nr1d4b (L)</td><td>Treat (–2.44x)</td><td>0.001</td><td>V<sub>A</sub></td><td>1.5E-05</td><td>4.2E-06 6.2E-05</td></tr><tr><td>Sex</td><td>0.001</td><td>V<sub>tank</sub></td><td>1.2E-05</td><td>3.4E-06 3.6E-05</td></tr><tr><td>Sex x Treat</td><td>0.008</td><td>V<sub>resid</sub></td><td>4.1E-05</td><td>1.5E-05 6.4E-05</td></tr><tr><td>Interactions</td><td></td><td></td><td></td><td></td></tr><tr><td>21°C: M vs F</td><td>0.014</td><td></td><td></td><td></td></tr><tr><td>31°C: M vs F</td><td>0.192</td><td></td><td></td><td></td></tr><tr><td>F: 21°C vs 31°C</td><td>0.001</td><td></td><td></td><td></td></tr><tr><td>M: 21°C vs 31°C</td><td>0.001</td><td></td><td></td><td></td></tr><tr><td rowspan="3">nsdhl (C)</td><td>Treat (2.38x)</td><td>0.001</td><td>V<sub>A</sub></td><td>1.4E-03</td><td>2.4E-04 0.004</td></tr><tr><td>Sex</td><td>0.134</td><td>V<sub>tank</sub></td><td>5.4E-04</td><td>2.0E-04 0.002</td></tr><tr><td>Sex x Treat</td><td>0.534</td><td>V<sub>resid</sub></td><td>2.4E-03</td><td>8.2E-04 0.004</td></tr><tr><td rowspan="3">pparab (D)</td><td>Treat (2.51x)</td><td>0.001</td><td>V<sub>A</sub></td><td>0.005</td><td>0.001 0.032</td></tr><tr><td>Sex</td><td>0.390</td><td>V<sub>tank</sub></td><td>0.005</td><td>0.002 0.018</td></tr><tr><td>Sex x Treat</td><td>0.066</td><td>V<sub>resid</sub></td><td>0.019</td><td>0.004 0.029</td></tr><tr><td rowspan="3">tefa (H)</td><td>Treat (1.35x)</td><td>0.001</td><td>V<sub>A</sub></td><td>0.008</td><td>0.001 0.035</td></tr><tr><td>Sex</td><td>0.001</td><td>V<sub>tank</sub></td><td>0.005</td><td>0.001 0.018</td></tr><tr><td>Sex x Treat</td><td>0.090</td><td>V<sub>resid</sub></td><td>0.021</td><td>0.005 0.032</td></tr><tr><td rowspan="3">tmx2b (A)</td><td>Treat (ns)</td><td>0.186</td><td>V<sub>A</sub></td><td>2.0E-06</td><td>7.8E-07 1.9E-05</td></tr><tr><td>Sex</td><td>0.010</td><td>V<sub>tank</sub></td><td>4.4E-06</td><td>9.8E-07 1.4E-05</td></tr><tr><td>Sex x Treat</td><td>0.442</td><td>V<sub>resid</sub></td><td>1.9E-05</td><td>9.8E-06 2.9E-05</td></tr><tr><td rowspan="3">ugt2a4 (I)</td><td>Treat (1.38x)</td><td>0.001</td><td>V<sub>A</sub></td><td>0.005</td><td>0.001 0.033</td></tr><tr><td>Sex</td><td>0.001</td><td>V<sub>tank</sub></td><td>0.005</td><td>0.001 0.021</td></tr><tr><td>Sex x Treat</td><td>0.130</td><td>V<sub>resid</sub></td><td>0.024</td><td>0.008 0.035</td></tr></tbody></table>

(continued)
<table><caption>**Table 1.** (continued)</caption><thead><tr><th rowspan="2">Gene</th><th rowspan="2">Fixed effects<br>Model term</th><th rowspan="2">P-Value</th><th colspan="4">Variance estimates</th></tr><tr><th>Post. Mode</th><th>95% PDIs</th></tr></thead><tbody><tr><td rowspan="3">uqcrc2b (E)</td><td>Treat (1.84x)</td><td>0.001</td><td>V<sub>A</sub></td><td>2.3E-04</td><td>6.3E-05</td><td>0.001</td></tr><tr><td>Sex</td><td>0.292</td><td>V<sub>tank</sub></td><td>3.3E-04</td><td>1.1E-04</td><td>0.001</td></tr><tr><td>Sex × Treat</td><td>0.732</td><td>V<sub>resid</sub></td><td>8.8E-04</td><td>4.2E-04</td><td>0.001</td></tr></tbody></table>

for the majority of genes: one to few families exhibited no effect of thermal treatment, whereas others significantly altered transcription in response to the temperature increase. In three genes (*cyp1a*, *gstr* & *uqcrc2b*), all families exhibited significant upregulation in response to increased temperature, although the magnitude of this effect differed.

## Functional annotations

We identified ten molecular pathways in which candidate genes (6 of 12) were found to co-occur with genes putatively under selection (Table 3) in the natural temperature gradient of Smith et al. (unpublished). These were largely dominated by metabolic processes, but also included generation of an innate immune response (NOD-like signalling pathway), and oocyte maturation. Enrichment analyses of overlapping GO categories revealed 28 functional clusters, with 193 of the 411 selected genes found to cluster with 8 candidate genes. Seven functional clusters comprised individual ontologies/categories with a mean fold enrichment of at least one (Table 4). Genes for transcriptional regulation and nucleotide binding were enriched at significantly higher rates than other categories, which included clusters associated with ion homeostasis and proteolysis.

Although six candidate genes mapped to a direct protein-protein interaction with at least one putatively selected gene, two major subnetworks were most prevalent: one in which hsp90aa1.2 appeared to occupy a central node, another implicating both hmgcs1 and pparab (Fig. 2). The hsp90aa1.2 subnetwork contained the most informative functional annotations (Fig. S3), with the majority of catalysis links implicated in the deacetylation of heat-shock factor 1 (*hsf1*) by sirtuin 1 (sirt1), and biochemical reactions associated with *hsf1* acetylation – in Fig. S3, these links are found mediated via the other heat-shock proteins and cognates (e.g. hsc70). The remaining catalysis reaction of this subnetwork, mediated via dynactin 1a (dctn1a), mapped to the phosphorylation of polo-like kinase 1 (*plk1*) by aurora kinase A (aurka). Annotations for the hmgcs1/pparab subnetwork were dominated by links of expression activation/induction (Fig. S4).

## Discussion

Clarifying the roles of phenotypic plasticity (environmentally induced response) and microevolution (genetic

response) in scenarios of environmental change is a major research priority (Hoffmann and Willi 2008; Merilä and Hendry 2014). Whilst under certain conditions plasticity can be adaptive (Gotthard and Nylin 1995; Ghalambor et al. 2007), genetic responses are generally a requisite for long-term persistence to predicted climate changes (Bradshaw and Holzapfel 2006). Moreover, the potential for an evolved response requires demonstration of heritable variation in traits affecting individual performance under new environmental conditions (Gienapp et al. 2008; Charmantier and Gienapp 2014). We can draw two conclusions from our long-term common-garden experiment of functional responses to predicted climate change with the Australian rainbowfish *M. duboulayi*. The first is that increases in seasonal austral summer temperature projected under a high emissions scenario are likely to have a significant effect on the transcriptional profiles of this subtropical rainbowfish. The second is that even in the face of significant temperature effects, there is a variation in both transcription and the transcriptional response to increased temperatures. Most importantly, this variation has an underlying genetic basis and, as such, might allow for an evolved response.

## Additive genetic variance and adaptive potential

Perhaps more important than confirmation of the transcriptional effects of projected temperature increases, this study adds to the growing body of evidence that transcriptional variance has an underlying heritable basis. Much of the previous work in the field has been dominated by studies designed to identify quantitative trait loci responsible for transcriptional variation (eQTL, Gibson and Weir 2005; Skelly et al. 2009). Whilst this line of enquiry has provided many valuable insights into the regulatory nature of gene expression, classical pedigree-based analyses of variance components may provide more reliable estimates of the fraction of transcriptional variance due to additive genetic effects (discussed in Leder et al. 2015). Thus, not only is this a highly cost-effective approach to understanding the genetics of transcriptional variance, but also it has the potential to yield more accurate estimates of the parameters required to predict a response to changing selection pressures. Yet to date surprisingly few studies investigating the transcriptional effects of changing environmental stressors
<figure><img src="image_3.png" /><figcaption>**Figure 1** Mean initial RNA concentrations ($N_0$), conditioned on additive genetic variance and random variation due to tank effects, as a function of sex and temperature treatment in the rainbowfish *Melanotaenia duboulayi*. Females are plotted in white and males in grey. Whiskers denote 95% PDIs; boxes indicate the quartile range of posterior estimates. (A) *tmx2b*; (B) *hmgcs1*; (C) *nsdhl*; (D) *pparab*; (E) *uqcrc2b*; (F) *agxta*; (G) *gstr*; (H) *tefa*; (I) *ugt2a4*; (J) *cyp1a*; (K) *hsp90aa1.2*; (L) *nr1d4b*. Coloured lines denote family-specific reaction norms for 6 families represented in each temperature treatment. Note that for clarity, only female-specific reaction norms are presented; plots of random intercept and slopes for both males and females can be found as online supplementary material (Figs S1 and S2).</figcaption></figure>
<table><caption>**Table 2.** Summary of gene effects in the rainbowfish *Melanotaenia duboulayi*. Significance of heritability estimates is evaluated by contrasting the deviance information criterion (DIC) for a fully parameterized model ($V_A$) with a simpler model (null) containing only random variance due to tank effects. Estimates of narrow-sense heritability ($h^2$) for candidate gene expression are bounded by 95% posterior density interval estimates (PDIs). Genes whose transcription was significantly affected by rearing temperature (Temp.) are indicated by T; note that for three genes, marked *, temperature effects were sex specific. Genes with significant heritable variation in expression and those with significant family-specific variation in temperature reaction norms ($G \times E$) are also marked T. Adaptive and/or functional significance was inferred on the basis of candidate gene co-occurrence with outlier genes in an annotated molecular pathway (Path.), a functional enrichment category (Enrich.), or in a predicted protein-protein interaction network (Prot.).</caption><thead><tr><th rowspan="2">Gene</th><th colspan="2">Model DIC</th><th colspan="3">Heritability (h<sup>2</sup>)</th><th colspan="3">Statistical significance</th><th colspan="3">Adaptive inference</th></tr><tr><th>Null</th><th>V<sub>A</sub></th><th>Est.</th><th>95% PDIs</th><th>Temp.</th><th>h<sup>2</sup></th><th>G &times; E</th><th>Path.</th><th>Enrich.</th><th>Prot.</th></tr></thead><tbody><tr><td>agxta</td><td>252.08</td><td>193.65</td><td>0.540</td><td>0.137</td><td>0.829</td><td>T</td><td>T</td><td>-</td><td>T</td><td>T</td><td>T</td></tr><tr><td>cyp1a</td><td>199.01</td><td>179.38</td><td>0.246</td><td>0.054</td><td>0.599</td><td>*</td><td>T</td><td>T</td><td>T</td><td>T</td><td>T</td></tr><tr><td>gstr</td><td>-21.84</td><td>-46.00</td><td>0.178</td><td>0.068</td><td>0.651</td><td>T</td><td>T</td><td>T</td><td>T</td><td>T</td><td>T</td></tr><tr><td>hmgcs1</td><td>223.59</td><td>182.54</td><td>0.205</td><td>0.083</td><td>0.767</td><td>T</td><td>T</td><td>T</td><td>T</td><td>-</td><td>T</td></tr><tr><td>hsp90aa1.2</td><td>-399.63</td><td>-493.37</td><td>0.518</td><td>0.197</td><td>0.876</td><td>*</td><td>T</td><td>-</td><td>T</td><td>T</td><td>T</td></tr><tr><td>nr1d4b</td><td>-724.48</td><td>-754.21</td><td>0.169</td><td>0.056</td><td>0.647</td><td>*</td><td>T</td><td>T</td><td>-</td><td>T</td><td>-</td></tr><tr><td>nsdhl</td><td>-287.79</td><td>-314.24</td><td>0.281</td><td>0.063</td><td>0.716</td><td>T</td><td>T</td><td>T</td><td>-</td><td>-</td><td>-</td></tr><tr><td>pparab</td><td>-71.46</td><td>-94.00</td><td>0.112</td><td>0.043</td><td>0.683</td><td>T</td><td>T</td><td>T</td><td>T</td><td>T</td><td>T</td></tr><tr><td>tefa</td><td>-55.05</td><td>-91.82</td><td>0.322</td><td>0.046</td><td>0.720</td><td>T</td><td>T</td><td>T</td><td>-</td><td>T</td><td>-</td></tr><tr><td>tx2b</td><td>-815.84</td><td>-827.23</td><td>0.106</td><td>0.029</td><td>0.516</td><td>-</td><td>T</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>ugt2a4</td><td>-53.16</td><td>-79.46</td><td>0.157</td><td>0.044</td><td>0.654</td><td>T</td><td>T</td><td>T</td><td>-</td><td>-</td><td>T</td></tr><tr><td>uqcrc2b</td><td>-420.25</td><td>-428.61</td><td>0.153</td><td>0.040</td><td>0.512</td><td>T</td><td>T</td><td>T</td><td>-</td><td>T</td><td>-</td></tr></tbody></table>

<table><caption>**Table 3.** Molecular pathways comprising at least one candidate gene assayed for transcriptional variation in the rainbowfish *Melanotaenia duboulayi*, and genes putatively under selection in a natural temperature gradient (Selected Partners; Smith et al. unpublished). Accession numbers for the Kyoto Encyclopedia of Genes and Genomes (KEGG) database are provided for pathways associated with zebrafish homologues.</caption><thead><tr><th>KEGG ID</th><th>Pathway</th><th>Candidates</th><th>Selected partners</th></tr></thead><tbody><tr><td>dre00140</td><td>Steroid hormone biosynthesis</td><td>cyp1a</td><td>cyp11a1</td></tr><tr><td>dre00250</td><td>Alanine, aspartate and glutamate metabolism</td><td>agxta</td><td>aldh4a1</td></tr><tr><td>dre00260</td><td>Glycine, serine and threonine metabolism</td><td>agxta</td><td>alas2</td></tr><tr><td>dre00280</td><td>Valine, leucine and isoleucine degradation</td><td>hmgcs1</td><td>acaa1</td></tr><tr><td>dre00480</td><td>Glutathione metabolism</td><td>gstr</td><td>mgst3, ggctb</td></tr><tr><td>dre00980</td><td>Metabolism of xenobiotics by cytochrome P450</td><td>cyp1a, gstr</td><td>mgst3</td></tr><tr><td>dre03320</td><td>PPAR signalling pathway</td><td>pparab</td><td>acaa1</td></tr><tr><td>dre04621</td><td>NOD-like receptor signalling pathway</td><td>hsp90aa1.2</td><td>erbb2ip</td></tr><tr><td>dre04914</td><td>Progesterone-mediated oocyte maturation</td><td>hsp90aa1.2</td><td>map2k1, rps6ka3b</td></tr><tr><td>dre04920</td><td>Adipocytokine signalling pathway</td><td>pparab</td><td>ptpn11a</td></tr></tbody></table>

have taken such a pedigree-based approach (though see Côté et al. 2007; Roberge et al. 2007; Debes et al. 2012).

Although this study does provide requisite proof of genetic variation, at present we do not know how transcriptional variance in these candidate genes might affect fitness. In fact, the general importance of transcriptional variance as a contributor to fitness remains a point of contention (Feder and Walser 2005; Evans 2015). Given the difficulty of quantifying fitness, it is perhaps unsurprising that few studies have attempted to integrate fitness metrics into the analysis of transcriptional variation in wild vertebrate models (but see McCairns and Bernatchez 2010). Nevertheless, we may infer that some contributions to fitness must be conferred given the signatures of selection observed in transcriptional variation in natural populations. For example, stickleback populations inhabiting a broad latitudinal gradient in northern Europe show signatures of directional selection in nearly 16% of all transcripts, with many enriched for functional groups associated with temperature stress (Leder et al. 2015). Likewise, 22% of transcripts assayed in killifish adapted to a steep temperature gradient along the coast of eastern North America also show evidence of selection (Whitehead and Crawford 2006). These studies demonstrate that transcriptional variation can clearly respond to selection and aid in the adaptation to divergent thermal environments. However, establishing a direct link to fitness remains an important next step if we are to understand its role in *M. duboulayi*'s adaptive potential in a changing thermal environment.
<table><caption>**Table 4.** Functionally enriched clusters based on overlapping gene ontology (GO) annotations of biological process and molecular function for candidate genes and genes identified as putatively under selection in a natural temperature gradient. Overall enrichment score (Enrich.), total number of candidate and selected genes with overlapping annotations (Genes), and the number of gene ontology biological processes and molecular functions comprising the cluster (No. GO) are presented. Fold Enrichment describes the average enrichment score for all GO terms comprising the cluster.</caption><thead><tr><th>Functional cluster</th><th>Enrich.</th><th>Candidates present</th><th>Genes</th><th>No. GO</th><th>Fold enrichment</th></tr></thead><tbody><tr><td>Transcriptional regulation</td><td>2.27</td><td>nr1d4b, pparab, tefa</td><td>45</td><td>4</td><td>1.57 (1.49–1.65)</td></tr><tr><td>Nucleotide binding</td><td>1.81</td><td>hsp90aa1.2</td><td>56</td><td>9</td><td>1.41 (1.26–1.52)</td></tr><tr><td>Nucleotide metabolism</td><td>0.92</td><td>gstr</td><td>8</td><td>33</td><td>3.27 (1.64–9.65)</td></tr><tr><td>Vitamin & cofactor binding</td><td>0.62</td><td>agxta</td><td>6</td><td>4</td><td>2.31 (1.11–3.19)</td></tr><tr><td>Cellular ion homeostasis</td><td>0.47</td><td>gstr</td><td>3</td><td>8</td><td>2.52 (1.76–3.40)</td></tr><tr><td>Metal ion binding</td><td>0.40</td><td>cyp1a, gstr, pparab, uqcrc2b</td><td>61</td><td>4</td><td>1.05 (1.01–1.12)</td></tr><tr><td>Peptidase proteolysis</td><td>0.18</td><td>uqcrc2b</td><td>13</td><td>3</td><td>0.99 (0.87–1.08)</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 2** Protein-protein interaction network involving candidate genes assayed for transcriptional variation in the rainbowfish *Melanotaenia duboulayi* (bold red text), and genes putatively under selection in a natural temperature gradient (Smith et al. unpublished). Note that detailed descriptions of network edges for the two major subnetworks (hsp90aa1.2 & pparab/hmgcs1), indicating the types of interactions can be found as online supplementary material (Figs S3 and S4).</figcaption></figure>

## Inference from predicted functional associations

Although function alone does not imply fitness, it does contribute to whole organism performance, and so may help to inform our understanding of the targets of selection (Dalziel et al. 2009). Many of the candidate genes are known to be directly involved in pathways associated with oxidative stress (e.g. *gstr*; Hermes-Lima and Zen-teno-Savín 2002) and the regulation of metabolism (e.g. *nr1d4b*; Desvergne et al. 2006), both physiological processes likely to be affected by a warming thermal environment. And whilst new analyses suggest that mRNA levels might be a far better predictor of protein abundance than previously believed (Li and Biggin 2015),

increased protein production is not necessarily equivalent to increased functional activity (Feder and Walser 2005). However, in evolved/steady-state comparisons, transcriptional variation can correspond remarkably well to differences in enzyme activity (Nikinmaa et al. 2013). Given the extensive acclimation period of this study, we are likely safe to assume that current transcriptional data are reflective of steady-state levels typical of summer conditions, both current and projected. As such, it is plausible that a degree of functional significance is captured in these data and that the most likely candidates in the physiologically mediated function/fitness nexus should be found at the functional intersection with genes under selection in the wild.
## A role for plasticity?

All genes exhibited significant heritability for transcription. Although the proportion of total phenotypic variance explained by additive genetic effects was moderate (Table 2), estimates are similar to that observed in various taxa (reviewed in Leder et al. 2015). Whether these are too low to effectively respond to selection exerted by increasing temperature remains an open question. For example, in UK populations of the great tit (*Parus major*), heritability for breeding time is similarly modest ($h^2$ = 0.16), and insufficient to account for shifts which appear to track environmentally induced changes in an important prey item's phenology (Charmantier et al. 2008). In this example, plasticity in breeding time appears to have facilitated the shift. Given the labile nature of transcription, this raises the possibility of plasticity also facilitating a shift towards new optima beyond the limits of existing genetic variation.

Plasticity in the form of acclimation seems an unlikely long-term strategy. Indeed, species with an evolutionary history in highly variable thermal environments (i.e. those frequently adjusting/acclimating to changing conditions) may currently exist within their tolerance limits for physiological stress (Tomanek 2010), and so may have limited plastic acclimation potential beyond this to meet additional physiological challenges imposed by climate change. The view with respect to transcriptional plasticity is also ambiguous. For example, a comparison of acute versus chronic (3 weeks) thermal stress in the bald notothen (*Pagothenia borchgrevinki*) revealed that mRNA levels of genes associated with ROS scavenging were increased in acute exposure but did not differ between control and chronic exposure; most critically, chronically exposed fish showed significantly higher levels of oxidative damage (Almroth et al. 2015). It should be noted, however, that this observation comes from a highly stenothermal species, a group which may lack a typical heat-shock response (Tomanek 2010). In contrast, brown trout (*Salmo trutta*) sampled from different thermal environments exhibit patterns of G × E for transcription suggestive of either a role for plasticity during adaptation to different temperatures, or the evolution of plasticity within them (Meier et al. 2014).

One possible caveat to the role of transcriptional plasticity in adaptation, however, may be seen in the yeast model. Although studied in the context of nutritional environmental quality and not temperature, significant G × E for transcription was largely restricted to genes unlikely to have unique essential functions, for example paralogous genes (Landry et al. 2006). If this observation also holds true for other taxa, this might limit the number and/or types of genes with adequate heritable plasticity. One likely candidate, however, remains the PPAR pathway. Transcriptional regulation via PPARs (e.g. *pparab*) may play a role in

acclimation to both temperature and nutritional stress (Seebacher et al. 2010). In freshwater sticklebacks, plasticity in the transcription of *pparaa* (a *pparab* paralogue) is thought to have played a role in the colonization of a broader range of thermal environments than experienced by the ancestral population (Morris et al. 2014).

In all cases, for plasticity to play a meaningful role in adaptation to a changing environment, it must move the phenotype towards the new fitness peak and it should be heritable (Gotthard and Nylin 1995; Ghalambor et al. 2007). In subtropical rainbowfish, differences in temperature reaction norms indicate considerable variation in the potential response to increasing temperature (Fig. S1). Likewise, the relatively wide-ranging posterior density intervals for certain families can be interpreted as a proxy for plasticity in the classical sense (i.e. the ability of a given genotype to produce a range of phenotypes in different environments). Moreover, we suggest that significant variation among families in the slope coefficient for nine of the twelve candidate genes is representative of broad-sense heritability for plasticity. As such, if plasticity itself becomes an important determinant of fitness, there may be sufficient genetic variation upon which selection might act. Given that increased temperature variations are also a feature of most projections for the region (Reisinger et al. 2014), such potential for adaptive plasticity may play an important role in the persistence of *M. duboulayi* along the east coast of Australia.

## Acclimation time influences transcription pattern

For many candidate genes, the transcriptional response to a season-length duration of increased temperature showed patterns similar to those exhibited following short-term (i.e. 14d) exposure (Smith et al. 2013). Three in particular (*hmgcs1*, *nr1d4b* & *nsdhl*) were not only consistent between studies of *M. duboulayi*, but also with respect to observations in threespine stickleback (*Gasterosteus aculeatus*) subjected to an acute temperature stress (Leder et al. 2015). That the range of fold-change values observed among these studies was also remarkably similar would suggest that these genes may play a general role in the thermal stress response and/or thermal acclimation in teleost fishes – indeed increased liver transcription of *hmgcs1* has also been reported in zebrafish (*Danio rerio*) exhibiting metabolic and oxidative stress (Hugo et al. 2012; Zhang et al. 2012; Tsedensodnom et al. 2013), both conditions associated with thermal stress. Similarly, two other candidate genes (*gstr* & *uqcrc2b*) were consistently upregulated in this and the previous study of *M. duboulayi*. However, seven of the twelve genes studied exhibited patterns of transcription opposite of that inferred by transfrag quantitation of RNA-Seq data.
Given the preponderance of transcription differences between sexes, and the fact that only males were assayed in the previous RNA-Seq study, sex-specific patterns are one potential explanation for observed discrepancies. However, only *txmx2b* and *hsp90aa1.2* show sufficient variability (i.e. overlapping posterior density intervals; Fig. 1A and K) from which biased sampling might potentially suggest increased transcription in response to temperature, as observed previously. Moreover, although males show a trend towards increased *hsp90aa1.2* transcription, this is neither significant (Table 1) nor within an order of magnitude of the degree of fold change previously reported (ca. 140×). Conversely, threespine sticklebacks have been shown to exhibit a ca. 10-fold increase of *hsp90aa1.2* in liver tissue immediately following a temperature increase (Leder et al. 2015). Likewise, killifish (*Austrofundulus limnaeus*) show levels of *hsp90α* induction that appear to increase over the first 24 h of thermal treatment, but maintain a relatively constant level of upregulation over a 14-d period (Podrabsky and Somero 2004). The common element between these experiments and previous observations in *M. duboulayi* is the relatively short duration of exposure, likely occurring within a window of physiological stress/disruption of homeostasis. In contrast, current observations were made on fish likely acclimated to the increased temperature.

Acclimation time is known to affect the expression of a suite of physiological traits (Schulte et al. 2011; Healy and Schulte 2012), and molecular responses to thermal challenge are no exception. With respect to mRNA transcription, expression is dependent not only upon acclimation time, but also the degree of thermal change experienced (Currie and Schulte 2014). Whilst the common, immediate response to increased temperature involves initiation of a stress response, there appears to be considerable variability in the length of time required for homeostasis to be re-established. For many temperate-zone fishes, including both stenothermal and eurythermal species, transcription-wide patterns of expression have been observed at steady state within a 4-week acclimation period (Vornanen et al. 2005; Logan and Somero 2010). In contrast, transcriptional patterns consistent with those of acute thermal stress persist in zebrafish following 4 weeks of exposure to elevated temperatures (Vergauwen et al. 2010). Little is known about the duration and/or periodicity of the molecular response to thermal stress in *M. duboulayi* specifically; as such, we can only assume that transcript levels following 80 days of acclimation to experimental temperatures are at steady state. Irrespective, given the choice of experimental temperatures, these observations are likely to capture some elements of the transcriptional mechanisms underlying acclimatization to future summer seasons.

It should also be noted that previous observations were made on wild-caught fish. Although acclimated to control temperature for 1 month prior to experimentation, it is unknown to what extent experimental individuals may have experienced different temperatures when inhabiting a naturally heterogeneous thermal environment. Temperature variation during embryonic development can have a significant impact on thermal acclimation potential and the transcriptomic response of adult fish to thermal stress (Schaefer and Ryan 2006; Scott and Johnston 2012). Thus, discrepancies could merely reflect differences between fish reared under fluctuating versus constant temperatures during a critical developmental stage. Whilst this may open reasonable debate regarding the extent to which observations on model animals reared under artificially homogeneous conditions can be extrapolated to wild populations, we would re-iterate that this is outweighed by a significant advantage: the ability to disentangle environmental and genetic sources of variation.

In conclusion, this study provides a tangible example that can be added to the debate surrounding adaptive evolution in response to climate change. Evidence of heritability in candidate gene expression, in addition to broad-sense genetic variation in transcriptional plasticity, speaks to the potential for an adaptive response to future increases in environmental temperature. Additionally, molecular pathways of candidate genes were found to co-occur with genes putatively under climate-mediated selection in wild *M. duboulayi* populations. Such findings are essential to meaningful discussions of adaptive change, but are sorely lacking in many studies of climate change adaptation. Although establishing a direct link between transcriptional variation and fitness remains essential if we are to understand its role in *M. duboulayi*'s adaptive potential in a changing thermal environment, our study represents important first steps towards that end.

## Acknowledgement

This study was funded by the Discovery Program of the Australian Research Council (ARC grant DP110101207 and DP150102903 to L. Beheregaray and L. Bernatchez). L. Beheregaray also acknowledges support from ARC FT130101068. We thank Leo O'Reilly for assistance with sampling and Hillary Mahon, Leslie Morrison and Simon Westergaard for contributing to fish husbandry. Animal ethical approval was received from Flinders University (AWC E342).

## Conflict of interest

The authors declare no conflict of interest.
## Data archiving

Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.h40d8

## Literature cited

Allen, G. R. 1989. Freshwater Fishes of Australia. THF Publications, Neptune, NJ.

Almroth, B. C., N. Asker, B. Wassmur, M. Rosengren, F. Jutfelt, A. Gräns, K. Sundell et al. 2015. Warmer water temperature results in oxidative damage in an Antarctic fish, the bald notothen. Journal of Experimental Marine Biology and Ecology **468**:130–137.

Becker, C. D., and R. G. Genoway 1979. Evaluation of the critical thermal maximum for determining thermal tolerance of freshwater fish. Environmental Biology of Fishes **4**:245–256.

Beheregaray, L. B., and P. Sunnucks 2000. Microsatellite loci isolated from *Odontesthes argentinensis* and the *O. perugiae* species group and their use in other South American silverside fish. Molecular Ecology **9**:629–631.

Bradshaw, W. E., and C. M. Holzapfel 2006. Evolutionary response to rapid climate change. Science **312**:1477–1478.

Caissie, D. 2006. The thermal regime of rivers: a review. Freshwater Biology **51**:1389–1406.

Chambers, L. E., R. Altwegg, C. Barbraud, P. Barnard, L. J. Beaumont, R. J. M. Crawford, J. M. Durant et al. 2013. Phenological changes in the southern hemisphere. PLoS ONE **8**:e75514.

Charmantier, A., and P. Gienapp 2014. Climate change and timing of avian breeding and migration: evolutionary versus plastic changes. Evolutionary Applications **7**:15–28.

Charmantier, A., R. H. McCleery, L. R. Cole, C. Perrins, L. E. B. Kruuk, and B. C. Sheldon 2008. Adaptive phenotypic plasticity in response to climate change in a wild bird population. Science **320**:800–803.

Chown, S. L., A. A. Hoffmann, T. N. Kristensen, M. J. Angilletta, N. C. Stenseth, and C. Pertoldi 2010. Adapting to climate change: a perspective from evolutionary physiology. Climate Research **43**:3–15.

Côté, G., G. Perry, P. Blier, and L. Bernatchez 2007. The influence of gene-environment interactions on GHR and IGF-I expression and their association with growth in brook charr, *Salvelinus fontinalis* (Mitchell). BMC Genetics **8**:87.

Currie, S., and P. M. Schulte 2014. Thermal stress. In: D. H. Evans, J. B. Claiborne, and S. Currie, eds. The Physiology of Fishes, 4th edn, pp. 257–287. CRC Press, Boca Raton, FL.

Dalziel, A. C., S. M. Rogers, and P. M. Schulte 2009. Linking genotypes to phenotypes and fitness: how mechanistic biology can inform molecular ecology. Molecular Ecology **18**:4997–5017.

Debes, P. V., E. Normandeau, D. J. Fraser, L. Bernatchez, and J. A. Hutchings 2012. Differences in transcription levels among wild, domesticated, and hybrid Atlantic salmon (*Salmo salar*) from two environments. Molecular Ecology **21**:2574–2587.

Desvergne, B., L. Michalik, and W. Wahli 2006. Transcriptional regulation of metabolism. Physiological Reviews **86**:465–514.

Diffenbaugh, N. S., and C. B. Field 2013. Changes in ecologically critical terrestrial climate conditions. Science **341**:486–492.

Dudgeon, D., A. H. Arthington, M. O. Gessner, Z. I. Kawabata, D. J. Knowler, C. Lévèque, R. J. Naiman et al. 2006. Freshwater biodiversity: importance, threats, status and conservation challenges. Biological Reviews **81**:163–182.

Durinck, S., Y. Moreau, A. Kasprzyk, S. Davis, B. De Moor, A. Brazma, and W. Huber 2005. BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis. Bioinformatics **21**:3439–3440.

Evans, T. G. 2015. Considerations for the use of transcriptomics in identifying the 'genes that matter' for environmental adaptation. Journal of Experimental Biology **218**:1925–1935.

Fangue, N. A., E. J. Osborne, A. E. Todgham, and P. M. Schulte 2011. The onset temperature of the heat-shock response and whole-organism thermal tolerance are tightly correlated in both laboratory-acclimated and field-acclimatized tidepool sculpins (*Oligocottus maculosus*). Physiological and Biochemical Zoology **84**:341–352.

Feder, M. E., and J. C. Walser 2005. The biological limitations of transcriptomics in elucidating stress and stress responses. Journal of Evolutionary Biology **18**:901–910.

Foll, M., and O. Gaggiotti 2008. A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics **180**:977–993.

Ghalambor, C. K., J. K. McKay, S. P. Carroll, and D. N. Reznick 2007. Adaptive versus non-adaptive phenotypic plasticity and the potential for contemporary adaptation in new environments. Functional Ecology **21**:394–407.

Gibson, G., and B. Weir 2005. The quantitative genetics of transcription. Trends in Genetics **21**:616–623.

Gienapp, P., C. Teplitsky, J. S. Alho, J. A. Mills, and J. Merilä 2008. Climate change and evolution: disentangling environmental and genetic responses. Molecular Ecology **17**:167–178.

Gotthard, K., and S. Nylin 1995. Adaptive plasticity and plasticity as an adaptation: a selective review of plasticity in animal morphology and life history. Oikos **74**:3–17.

Hadfield, J. D. 2010. MCMC methods for multi-response generalized linear mixed models: the MCMCglmm R package. Journal of Statistical Software **33**:1–22.

Hattori, A., and K. Warburton 2003. Microhabitat use by the rainbow-fish *Melanotaenia duboulayi* in a subtropical Australian stream. Journal of Ethology **21**:15–22.

Healy, T. M., and P. M. Schulte 2012. Thermal acclimation is not necessary to maintain a wide thermal breadth of aerobic scope in the common killifish (*Fundulus heteroclitus*). Physiological and Biochemical Zoology **85**:107–119.

Helmuth, B. 2009. From cells to coastlines: how can we use physiology to forecast the impacts of climate change? Journal of Experimental Biology **212**:753–760.

Hermes-Lima, M., and T. Zenteno-Savín. 2002. Animal response to drastic changes in oxygen availability and physiological oxidative stress. Comparative Biochemistry and Physiology C – Toxicology & Pharmacology **133**:537–556.

Hoffmann, A. A., and C. M. Sgrò 2011. Climate change and evolutionary adaptation. Nature **470**:479–485.

Hoffmann, A. A., and Y. Willi 2008. Detecting genetic responses to environmental change. Nature Reviews Genetics **9**:421–432.

Hofmann, G. E., and A. E. Todgham 2010. Living in the now: physiological mechanisms to tolerate a rapidly changing environment. Annual Review of Physiology **72**:127–145.

Holdway, D. A., J. Hefferman, and A. Smith 2008. Multigeneration assessment of nonylphenol and endosulfan using a model Australian freshwater fish, *Melanotaenia fluviatilis*. Environmental Toxicology **23**:253–262.

Huang, D. W., B. T. Sherman, Q. Tan, J. R. Collins, W. G. Alvord, J. Roayaei, R. Stephens et al. 2007. The DAVID gene functional classification
cation tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biology **8**:R183.

Hugo, S. E., L. Cruz-Garcia, S. Karanth, R. M. Anderson, D. Y. R. Stainier, and A. Schlegel 2012. A monocarboxylate transporter required for hepatocyte secretion of ketone bodies during fasting. *Genes & Development* **26**:282–293.

Hurwood, D. A., and J. M. Hughes 2001. Historical interdrainage dispersal of eastern rainbowfish from the Atherton Tableland, north-eastern Australia. *Journal of Fish Biology* **58**:1125–1136.

Jones, O. R., and J. L. Wang 2010. COLONY: a program for parentage and sibship inference from multilocus genotype data. *Molecular Ecology Resources* **10**:551–555.

Kanehisa, M., and S. Goto 2000. KEGG: Kyoto Encyclopedia of Genes and Genomes. *Nucleic Acids Research* **28**:27–30.

Kanehisa, M., S. Goto, Y. Sato, M. Kawashima, M. Furumichi, and M. Tanabe 2014. Data, information, knowledge and principle: back to metabolism in KEGG. *Nucleic Acids Research* **42**:D199–D205.

King, S., and K. Warburton 2007. The environmental preferences of three species of Australian freshwater fish in relation to the effects of riparian degradation. *Environmental Biology of Fishes* **78**:307–316.

Landry, C. R., J. Oh, D. L. Hartl, and D. Cavalieri 2006. Genome-wide scan reveals that genetic variation for transcriptional plasticity in yeast is biased towards multi-copy and dispensable genes. *Gene* **366**:343–351.

Leder, E. H., R. J. S. McCairns, T. Leinonen, J. M. Cano, H. M. Viitaniemi, M. Nikinmaa, C. R. Primmer et al. 2015. The evolution and adaptive potential of transcriptional variation in sticklebacks – signatures of selection and widespread heritability. *Molecular Biology and Evolution* **32**:674–689.

Li, J. J., and M. D. Biggin 2015. Statistics requantitates the central dogma. *Science* **347**:1065–1066.

Logan, C. A., and G. N. Somero 2010. Transcriptional responses to thermal acclimation in the eurythermal fish *Gillichthys mirabilis* (Cooper 1864). American Journal of Physiology – Regulatory Integrative and Comparative Physiology **299**:R843–R852.

López-Maury, L., S. Marguerat, and J. Bähler 2008. Tuning gene expression to changing environments: from rapid responses to evolutionary adaptation. *Nature Reviews Genetics* **9**:583–593.

McCairns, R. J. S., and L. Bernatchez 2010. Adaptive divergence between freshwater and marine sticklebacks: insights into the role of phenotypic plasticity from an integrated analysis of candidate gene expression. *Evolution* **64**:1029–1047.

McGuigan, K., D. Zhu, G. R. Allen, and C. Moritz 2000. Phylogenetic relationships and historical biogeography of melanotaeniid fishes in Australia and New Guinea. *Marine and Freshwater Research* **51**:713–723.

Meehl, G. A., T. F. Stocker, W. D. Collins, P. Friedlingstein, A. T. Gaye, J. M. Gregory, A. Kitoh et al. 2007. Global climate projections. In S. Solomon, D. Qin, M. Manning, M. Marquis, K. Averyt, M. M. B. Tignor, H. L. Miller, and Z. L. Chen, eds. *Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change*, pp. 747–845. Cambridge University Press, Cambridge, UK.

Meier, K., M. M. Hansen, E. Normandeau, K. L. D. Mensberg, J. Frydenberg, P. F. Larsen, D. Bekkevold et al. 2014. Local adaptation at the transcriptome level in brown trout: evidence from early life history temperature genomic reaction norms. *PLoS ONE* **9**:e85171.

Merilä, J., and A. P. Hendry 2014. Climate change, adaptation, and phenotypic plasticity: the problem and the evidence. *Evolutionary Applications* **7**:1–14.

Morris, M. R. J., R. Richard, E. H. Leder, R. D. H. Barrett, N. Aubin-Horth, and S. M. Rogers 2014. Gene expression plasticity evolves in response to colonization of freshwater lakes in threespine stickleback. *Molecular Ecology* **23**:3226–3240.

Muñoz, N. J., R. D. Breckels, and B. D. Neff 2012. The metabolic, locomotor and sex-dependent effects of elevated temperature on Trinidadian guppies: limited capacity for acclimation. *Journal of Experimental Biology* **215**:3436–3441.

Nikinmaa, M., R. J. S. McCairns, M. W. Nikinmaa, K. A. Vuori, M.kanerva, T. Leinonen, C. R. Primmer et al. 2013. Transcription and redox enzyme activities: comparison of equilibrium and disequilibrium levels in the three-spined stickleback. *Proceedings of the Royal Society B-Biological Sciences* **280**:20122974.

van Oosterhout, C., W. F. Hutchinson, D. P. M. Wills, and P. Shipley 2004. MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. *Molecular Ecology Notes* **4**:535–538.

Patra, R. W., J. C. Chapman, R. P. Lim, and P. C. Gehrke 2007. The effects of three organic chemicals on the upper thermal tolerances of four freshwater fishes. *Environmental Toxicology and Chemistry* **26**:1454–1459.

Pearson, R. G., J. C. Stanton, K. T. Shoemaker, M. E. Aiello-Lammens, P. J. Ersts, N. Horning, D. A. Fordham et al. 2014. Life history and spatial traits predict extinction risk due to climate change. *Nature Climate Change* **4**:217–221.

Podrabsky, J. E., and G. N. Somero 2004. Changes in gene expression associated with acclimation to constant temperatures and fluctuating daily temperatures in an annual killifish *Austrofundulus limnaeus*. *Journal of Experimental Biology* **207**:2237–2254.

Pörtner, H. O., and R. Knust 2007. Climate change affects marine fishes through the oxygen limitation of thermal tolerance. *Science* **315**:95–97.

Pörtner, H. O., A. F. Bennett, F. Bozinovic, A. Clarke, M. A. Lardies, M. Lucassen, B. Pelster et al. 2006. Trade-offs in thermal adaptation: the need for a molecular to ecological integration. *Physiological and Biochemical Zoology* **79**:295–313.

Pusey, B. J., A. H. Arthington, J. R. Bird, and P. G. Close 2001. Reproduction in three species of rainbowfish (Melanotaeniidae) from rain-forest streams in northern Queensland, Australia. *Ecology of Freshwater Fish* **10**:75–87.

Ramakers, C., J. M. Ruijter, R. H. L. Deprez, and A. F. M. Moorman 2003. Assumption-free analysis of quantitative real-time polymerase chain reaction (PCR) data. *Neuroscience Letters* **339**:62–66.

Reisinger, A., R. L. Kitching, F. Chiew, L. Hughes, P. C. D. Newton, S. S. Schuster, A. Tait et al. 2014. Australasia. In V. R. Barros, C. B. Field, D. J. Dokken, M. D. Mastrandrea, K. J. Mach, T. E. Bilir, M. Chatterjee, K. L. Ebi, Y. O. Estrada, R. C. Genova, B. Girma, E. S. Kissel, A. N. Levy, S. MacCracken, P. R. Mastrandrea, and L. L. White, eds. *Climate Change 2014: Impacts, Adaptation, and Vulnerability*. Part B: Regional Aspects. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel of Climate Change, pp. 1371–1438. Cambridge University Press, Cambridge, UK.

Roberge, C., H. Guderley, and L. Bernatchez 2007. Genomewide identification of genes under directional selection: gene transcription *QST* scan in diverging Atlantic salmon subpopulations. *Genetics* **177**:1011–1022.

Ruijter, J. M., M. W. Pfaffl, S. Zhao, A. N. Spiess, G. Boggy, J. Blom, R. G. Rutledge et al. 2013. Evaluation of qPCR curve analysis methods for reliable biomarker discovery: bias, resolution, precision, and implications. *Methods* **59**:32–46.
Schaefer, J., and A. Ryan 2006. Developmental plasticity in the thermal tolerance of zebrafish *Danio rerio*. *Journal of Fish Biology* **69**:722–734.

Schulte, P. M., T. M. Healy, and N. A. Fangue 2011. Thermal performance curves, phenotypic plasticity, and the time scales of temperature exposure. *Integrative and Comparative Biology* **51**:691–702.

Scott, G. R., and I. A. Johnston 2012. Temperature during embryonic development has persistent effects on thermal acclimation capacity in zebrafish. *Proceedings of the National Academy of Sciences of the United States of America* **109**:14247–14252.

Seebacher, F., M. D. Brand, P. L. Else, H. Guderley, A. J. Hulbert, and C. D. Moyes 2010. Plasticity of oxidative metabolism in variable climates: molecular mechanisms. *Physiological and Biochemical Zoology* **83**:721–732.

Settele, J., R. Scholes, R. Betts, S. E. Bunn, P. Leadley, D. Nepstad, J. T. Overpeck et al. 2014. Terrestrial and inland water systems. In C. B. Field, V. R. Barros, D. J. Dokken, K. J. Mach, M. D. Mastrandrea, T. E. Bilir, M. Chatterjee, K. L. Ebi, Y. O. Estrada, R. C. Genova, B. Girma, E. S. Kissel, A. N. Levy, S. MacCracken, P. R. Mastrandrea, and L. L. White, eds. *Climate Change 2014: Impacts, Adaptation, and Vulnerability*. Part A: Global and Sectoral Aspects. *Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel of Climate Change*, pp. 271–359. Cambridge University Press, Cambridge, UK.

Sherman, B. T., D. W. Huang, Q. N. Tan, Y. J. Guo, S. Bour, D. Liu, R. Stephens et al. 2007. DAVID knowledgebase: a gene-centered database integrating heterogeneous gene annotation resources to facilitate high-throughput gene functional analysis. *BMC Bioinformatics* **8**:426.

Skelly, D. A., J. Ronald, and J. M. Akey 2009. Inherited variation in gene expression. *Annual Review of Genomics and Human Genetics* **10**:313–332.

Smith, S., L. Bernatchez, and L. B. Beheregaray 2013. RNA-seq analysis reveals extensive transcriptional plasticity to temperature stress in a freshwater fish species. *BMC Genomics* **14**:375.

Snel, B., G. Lehmann, P. Bork, and M. A. Huynen 2000. STRING: a web-server to retrieve and display the repeatedly occurring neighbourhood of a gene. *Nucleic Acids Research* **28**:3442–3444.

Sunnucks, P., and D. F. Hales 1996. Numerous transposed sequences of mitochondrial cytochrome oxidase I-II in aphids of the genus *Sitobion* (Hemiptera: Aphididae). *Molecular Biology and Evolution* **13**:510–524.

Szklarczyk, D., A. Franceschini, S. Wyder, K. Forslund, D. Heller, J. Huerta-Cepas, M. Simonovic et al. 2015. STRING v10: protein-protein interaction networks, integrated over the tree of life. *Nucleic Acids Research* **43**:D447–D452.

Tomanek, L. 2010. Variation in the heat shock response and its implication for predicting the effect of global climate change on species' biogeographical distribution ranges and metabolic costs. *Journal of Experimental Biology* **213**:971–979.

Tsedensodnom, O., A. M. Vacaru, D. L. Howarth, C. Y. Yin, and K. C. Sadler 2013. Ethanol metabolism and oxidative stress are required for unfolded protein response activation and steatosis in

zebrafish with alcoholic liver disease. *Disease Models & Mechanisms* **6**:1213–1226.

Unmack, P. J., G. R. Allen, and J. B. Johnson 2013. Phylogeny and bio-geography of rainbowfishes (*Melanotaeniidae*) from Australia and New Guinea. *Molecular Phylogenetics and Evolution* **67**:15–27.

Vandesompele, J., K. De Preter, F. Pattyn, B. Poppe, N. Van Roy, A. De Paepe, and F. Speleman. 2002. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. *Genome Biology* **3**:research0034.1–0034.11. doi: 10.1186/gb-2002-3-7-research0034.

Vergauwen, L., D. Benoot, R. Blust, and D. Knapen 2010. Long-term warm or cold acclimation elicits a specific transcriptional response and affects energy metabolism in zebrafish. *Comparative Biochemistry and Physiology A - Molecular & Integrative Physiology* **157**:149–157.

Vornanen, M., M. Hassinen, H. Koskinen, and A. Krasnov 2005. Steady-state effects of temperature acclimation on the transcriptome of the rainbow trout heart. *American Journal of Physiology – Regulatory Integrative and Comparative Physiology* **289**:R1177–R1184.

Whitehead, A., and D. L. Crawford 2006. Neutral and adaptive variation in gene expression. *Proceedings of the National Academy of Sciences of the United States of America* **103**:5425–5430.

Willmer, P., G. Stone, and I. Johnston 2005. *Environmental Physiology of Animals*. Blackwell Publishing, Malden, MA.

Young, M. J., L. W. Simmons, and J. P. Evans 2009. Isolation and characterization of 12 novel DNA microsatellites in the western rainbowfish, *Melanotaenia australis*. *Molecular Ecology Resources* **9**:1252–1254.

Zhang, W., Y. Liu, H. X. Zhang, and J. Y. Dai 2012. Proteomic analysis of male zebrafish livers chronically exposed to perfluorononanoic acid. *Environment International* **42**:20–30.

Zhu, D. Q., S. Degnan, and C. Moritz 1998. Evolutionary distinctiveness and status of the endangered Lake Eacham rainbowfish (*Melanotaenia eachamensis*). *Conservation Biology* **12**:80–93.

## Supporting Information

Additional Supporting Information may be found in the online version of this article:

<caption>**Table S1.** Primer sequences of reference and candidate genes.</caption>

**Table S2** Putative protein family identity and functional role(s) of candidate genes in the rainbowfish *Melanotaenia duboulayi*.

**Table S3.** Variance estimates from nested mixed-effects models in the rainbowfish *Melanotaenia duboulayi*.

**Figure S1.** Family-specific reaction norms for genes with significant $G \times E$ for expression in the rainbowfish *Melanotaenia duboulayi*.

**Figure S2.** Family-specific reaction norms for genes with non-significant or marginal/questionable $G \times E$ for expression in the rainbowfish *Melanotaenia duboulayi*.

**Figure S3.** Interactions between proteins for the hsp90aa1.2 sub-network.

**Figure S4.** Interactions between proteins for the pparab/hmgcs1 sub-network.