Regulatory mechanisms of metabolic flexibility in the dark-eyed junco (Junco hyemalis)

Phenotypic flexibility is an organism's ability to reversibly modify its phenotype with respect to changes in its environment. Because it allows an organism to optimally match its phenotype to prevailing biotic and abiotic conditions, phenotypic flexibility has profound consequences for organismal fitness (Piersma and Drent, 2003). Elucidating the mechanisms that underlie phenotypic flexibility therefore has important implications for understanding how organisms cope with rapid environmental change.

Temporal variation in organismal bioenergetics is a prominent example of phenotypic flexibility. Endotherms reversibly alter their metabolic rates across the annual cycle to meet the seasonal demands associated with migration, reproduction and thermoregulation. This metabolic flexibility is especially important for small endotherms that are resident in seasonally variable environments because metabolic thermoregulatory performance has acute consequences for survival during prolonged cold stress (Hayes and O'Conner, 1999). Birds that winter at high latitudes can provide a valuable model for understanding the physiological mechanisms underlying metabolic flexibility because they must rely largely on physiological adjustments that enhance thermogenic performance for survival. For instance, in winter, temperate-zone birds elevate their capacities for shivering thermogenesis to maintain thermoregulatory homeostasis in the face of falling temperatures and this increased thermogenic performance is associated with elevated summit metabolic rates (e.g. Arens and Cooper, 2005). However, winter is also accompanied by reductions in productivity, decreased foraging time and increased overnight fasting periods, so birds must meet these increased thermoregulatory demands in the face of reduced energetic input (Swanson, 2010). Thus, short-term physiological adjustments that optimally balance competing energetic demands are crucial for survival during the harsh winters in high-latitude environments.

At the biochemical level, variation in aerobic thermogenic performance is manifest through differences in the expression and functional properties of the genes encoding the machinery of cellular metabolism. For instance, deer mice (Peromyscus maniculatus) living at high altitude exhibit greater thermogenic capacities compared with their lowland counterparts, and this enhanced performance has been linked to regulatory changes in the expression of genes that influence the use of metabolic fuel and oxygen (Cheviron et al., 2012,, 2014). Because changes in gene regulation allow the genome to rapidly and reversibly respond to environmental variation, transcriptomic studies can provide key insights into the mechanistic underpinnings of phenotypic flexibility and the adaptive modification of complex physiological traits, both of which require concerted changes in hierarchical and interacting regulatory networks (Cheviron and Brumfield, 2012). The hierarchical nature of complex physiological traits is well illustrated by the following three levels of physiological adjustment that enhance aerobic thermogenic performance (Fig. 1A).

Small birds rely on flight muscles for facultative thermogenesis (Yacoe and Dawson, 1983), and several studies have demonstrated a prominent role for avian pectoralis size on thermogenic capacity (e.g. Swanson et al., 2013). Muscle hypertrophy can be achieved by altering the expression of genes in key growth signaling pathways, such as the rapamycin (mTOR) pathway (regulation of protein synthesis in adult mammalian skeletal muscle) and the insulin-like growth factor 1 (IGF1) signaling pathway (stimulates cell growth and proliferation). Inhibitory growth factors, such as myostatin, can interact with these pathways to impede protein synthesis and induce muscle atrophy; therefore, blocking myostatin activity can result in subsequent growth to myocytes (e.g. McPherron et al., 1997). In principle, modifications could also be made to the relative proportion of fiber types within the muscle; however, numerous studies have shown that the pectoralis of small birds is composed almost exclusively of fast oxidative-glycolytic fibers (e.g. Welch and Altshuler, 2009). Therefore, although it is possible that the avian pectoralis does grow and atrophy seasonally in accordance with changing thermogenic demands, the fiber composition does not seem to be altered to any appreciable degree.

Increased oxygen delivery can be accomplished by a number of cardiovascular changes throughout the body. Within the pectoralis, oxygen diffusion distance can be reduced through increased capillary density in the muscle, as has been shown in altitude-acclimated pigeons (Mathieu-Costello and Agey, 1997). Angiogenesis is primarily stimulated by the vascular endothelial growth factor (VEGF) signaling pathway, which mediates endothelial cell proliferation and migration (Ferrara and Davis-Smyth, 1997). The VEGF pathway in turn can be activated by cellular hypoxia via the hypoxia-inducible factor 1 (HIF-1) signaling cascade (Bruick and McKnight, 2001). In addition to initiating angiogenesis, the HIF-1 cascade regulates oxygen homeostasis in a number of other ways, including eliciting erythropoiesis via hormonal control (Rankin et al., 2007) and inducing vasodilation via nitric oxide production (Beall et al., 2012). Indeed, elevated hematocrit is associated with the winter phenotype in birds and is also positively correlated with thermogenic capacity (Swanson, 1990a; Petit and Vézina, 2014). These many interacting mechanisms can serve to boost oxygen delivery to the muscle, enabling individuals to sustain increased aerobic performance.

Increases in fuel supplies and oxygen utilization can increase the cellular aerobic capacity of thermogenic organs, such as the pectoralis. Energy-rich lipids are the primary fuel source for avian skeletal muscle across a wide range of exercise intensities (Guglielmo, 2010) and adipose stores may contribute 80% of the energy required for avian shivering (Vaillancourt et al., 2005). Fuel availability is not only dependent upon the size of adipose stores, but also upon an individual's ability to rapidly mobilize, catabolize and oxidize these fuels. Mobilization occurs by modifying expression to the many biochemical pathways that comprise the greater fatty acid supply pathway (Fig. 1B). This includes carrier proteins, which are required to move fatty acids to the muscle and across the cellular membrane to the mitochondria. To that effect, birds meet the heightened demands of winter by seasonally increasing fatty acid transporter proteins in the flight muscle (Liknes et al., 2014). Moreover, increased metabolic intensity of the pectoralis can be achieved by increasing mitochondrial density (e.g. Evans et al., 1992) or by elevating activities of catabolic enzymes involved in fatty acid β-oxidation (e.g. Guglielmo et al., 2002). Birds also exhibit corresponding seasonal changes in activities of oxidative enzymes that are involved in the final steps of aerobic metabolism – the citric acid cycle and oxidative phosphorylation pathways (e.g. Liknes and Swanson, 2011; Zheng et al., 2014). Thus, adjustments can be made to fine-tune fuel supplies and oxidative utilization at many points along this extensive cascade.

Despite this range of possibilities, it is still unclear how organisms effectively modify and combine these modular and hierarchical trait components in order to rapidly respond to exogenous factors. Here, we report the results of an experiment designed to gain insight into the fundamental mechanisms of avian metabolic flexibility, integrating new data on genomic transcriptional profiles with previously published data on thermogenic performance and enzyme activities for the same individuals (Swanson et al., 2014a). This work focuses on the dark-eyed junco (Junco hyemalis), a small passerine that winters from Alaska to Arizona (Nolan et al., 2002). Wintering juncos exhibit increases in thermogenic performance and cold tolerance compared with summer-acclimatized birds (Swanson, 1991,, 2001). To induce changes in thermogenic performance, we exposed wild-caught juncos to synthetic photoperiod and temperature regimes, two seasonal cues that have been previously implicated in altering avian metabolic rate and thermogenic performance (Carey and Dawson, 1999; Swanson, 2010). We then related variation in thermogenic performance to differences in underlying gene expression profiles of a plastic, thermogenically active tissue, the pectoralis. We explicitly surveyed the genes associated with the cell-, tissue- and organ-level responses outlined above to test whether juncos made regulatory changes at each of the three levels of physiological organization. By concurrently examining the effects of photoperiod and temperature, we address how each environmental stimulus induces regulatory changes in the pectoralis and how these transcriptomic responses, in turn, affect whole-organism thermogenic performance.

We exposed juncos to synthetic temperature and photoperiod cues using a two-by-two experimental design in order to induce changes in thermogenic performance. As previously reported (Swanson et al., 2014a), these acclimation treatments effectively altered junco summit metabolic rate (ΔM_sum; Fig. 2). Cold- and warm-acclimated individuals exhibited M_sum similar to wild-caught winter and summer juncos, respectively (Swanson, 1990b). Cold-acclimated birds increased M_sum by 1.19 ml O₂ min⁻¹ compared with warm-acclimated birds (ANOVA F_1,35=33.998, P=1.29×10⁻⁶), whereas photoperiod treatments did not result in changes to M_sum (F_1,35=0.168, P=0.684), nor was the interaction of the two treatments significant (F_1,35=0.994, P=0.326). Mass-adjusted pectoralis size (M_p) did not vary among acclimation treatments (ANOVA temperature, F_1,35=2.352, P=0.134; photoperiod, F_1,35=1.592, P=0.215; interaction, F_1,35=0.221, P=0.641; supplementary material Fig. S1).

We used a general linear model approach to identify genes that exhibited a significant difference in transcript abundance due to acclimation treatment. A total of 1325 genes (12.2% of total transcripts) were differentially expressed (DE) among photoperiod treatments [false discovery rate (FDR)<0.05; supplementary material Table S2]. Within this subset of photoperiod-sensitive DE genes (those with significant photoperiod effects), we identified 51 transcriptional modules of coregulated genes using modulated modularity clustering (MMC; Stone and Ayroles, 2009; Fig. 3A). We then used principal component and regression analyses to determine which of these transcriptional modules were associated with ΔM_sum. Overall module expression (summarized as PC1 scores) correlated with individual ΔM_sum for three transcriptional modules (Table 1). These modules were significantly enriched for gene ontology (GO) terms related to glycolysis (P36) and angiogenesis (P45; supplementary material Table S3). Similarly, five transcriptional modules showed strong associations with M_p, and these modules were significantly enriched for terms involved in ion transport (P8, P21) and lipid metabolism (P32, P51).

A total of 756 genes (7.0% of total transcripts) were differentially expressed across the temperature treatments (supplementary material Table S4). Of these, 196 DE genes were also differentially expressed in response to photoperiod manipulation (supplementary material Fig. S2). Temperature-sensitive DE genes (those with significant temperature effects) clustered into 38 transcriptional modules (Fig. 3B) and 23 of these temperature-sensitive transcriptional modules exhibited strong associations with ΔM_sum. These modules were enriched for genes annotated with terms related to lipid metabolism (N=5 modules), oxidative phosphorylation (N=5), glycolysis (N=2), angiogenesis and oxygen transport (N=4), and growth (N=7; Table 2). Additionally, three transcriptional modules showed strong associations with M_p. These modules were enriched for terms involved in muscle contraction (T5), and growth and angiogenesis (T38).

Our dataset included a majority of genes (84–100%) from each of nine candidate regulatory pathways that we examined in entirety (Table 3). A high proportion of the genes involved in fatty acid metabolism (27%), glycerolipid metabolism (27%), glycolysis (33%) and pyruvate metabolism (29%) pathways exhibited significant differential expression among temperature treatments. Photoperiod also had a strong effect on the expression of genes involved in fatty acid metabolism (27% of genes differentially expressed). We used a binomial test to examine large-scale changes in gene expression across each pathway. This analysis revealed evidence of concerted expression changes in response to photoperiod across three integrated pathways (fatty acid metabolism, citric acid cycle, and oxidative phosphorylation; Table 3). A disproportionate number of genes exhibited trends towards upregulation in short-day treatments (SD) across both the citric acid cycle and oxidative phosphorylation, whereas the reverse pattern [upregulation in long-day treatments (LD)] was exhibited for the fatty acid oxidation pathways. Additionally, the fatty acid metabolism pathway exhibited concerted expression in response to temperature (upregulation in cold) and the citric acid cycle also exhibited a non-significant trend in the same direction (P=0.057).

Expression patterns were largely paralleled by changes in the activities of key enzymes in relevant metabolic pathways (Fig. 4). We quantified enzyme activities for enzymes involved in fatty acid transport capacity (carnitine palmitoyl transferase; CPT), fatty-acid oxidation capacity (β-hydroxyacyl coenzyme-A dehydrogenase; HOAD), and cellular metabolic intensity (citrate synthase; CS). Enhanced capacity for aerobic thermogenesis was associated with elevated activities of CPT alone (P=0.044; supplementary material Fig. S3) and CPT was also the only enzyme to reveal temperature-sensitive activity patterns (P=0.039; supplementary material Fig. S4). In contrast, CS and HOAD enzyme activities exhibited photoperiod effects, such that LD individuals demonstrated increased activity of both enzymes (P_CS=0.039, P_HOAD=0.001), but neither was associated with ΔM_sum.

Because transcriptomic changes do not necessarily reflect corresponding variation at the proteomic level (e.g. Feder and Walser, 2005), we then tested for correlations within our own dataset by examining the relationship between transcript abundance and enzyme activity using linear regression. Transcript abundance for the gene encoding CS correlated with CS activity (Fig. 4). If the outlying CS activity datum (Grubbs' test, P=4.03×10⁻⁵) is removed, this relationship is no longer statistically significant (R²=−0.023, P=0.659); however, we have no obvious reason to exclude this result. Similarly, each of the three genes encoding HOAD subunits correlated with HOAD activity (both with and without the outlying individual). However, CPT activity was not significantly correlated with abundance of either the CPT1a or CPT2 transcript.

Phenotypic flexibility allows an organism to optimally match its phenotype to prevailing environmental conditions. This reversible modification of complex physiological traits can be driven by changes in the regulation of genes at several levels of hierarchical organization (Cheviron and Brumfield, 2012). Our acclimation treatments successfully altered junco thermogenic performance, and we found that changes in performance were associated with underlying transcriptomic changes in interacting hierarchical pathways. Seasonal cues associated with winter induced distinct transcriptomic responses, highlighting the importance of temperature and photoperiod cues in mediating distinct, seasonal responses in avian physiology. Cold-acclimated individuals exhibited elevated thermogenic capacities and transcriptomic signatures indicative of enhanced O₂ delivery, as well as enhanced transport and β-oxidation of fatty acids. These differences in transcript abundances were largely associated with the activities of enzymes that influence flux through the fatty acid oxidation and citric acid cycle pathways, suggesting that changes in gene regulation in these pathways translate into functional variation at the proteomic level. Conversely, changes in photoperiod were associated with differential regulation of pathways that influence cellular oxidative capacities. However, photoperiod manipulation was not associated with significant changes in thermogenic performance. Collectively, these results provide important clues as to the differential regulation of hierarchal, multistep pathways that influence complex physiological adjustments to abiotic stressors.

As in other bird species (e.g. Swanson, 2010), winter increases in junco thermogenic performance are often accompanied by pectoralis muscle hypertrophy (Swanson, 1991). However, we found that M_p was not affected by acclimation treatment in this study, thus increases in thermogenic performance were not accompanied by corresponding increases in pectoralis size (supplementary material Fig. S1). The lack of changes in M_p as a driver of improved thermogenic performance further emphasizes the importance of cellular adjustments that enhance lipid oxidation.

Despite the absence of changes in M_p across the acclimation treatments, we still found transcriptomic responses, suggesting either the induction of muscle growth or muscle repair. Both ΔM_sum and M_p were correlated with the expression of transcriptomic modules containing genes related to muscle development. Temperature-sensitive modules associated with thermogenic capacity were enriched for genes annotated with collagen (two in T28), myoblast proliferation (T22), skeletal muscle cell differentiation (T23), transforming growth factor (T29, T33) and regulation of growth (T38) terms. Some of these same terms, in combination with that of inflammatory response (T20, T22, T24, T27), are implicated in mammalian muscle repair (Karalaki et al., 2009). Considering that pectoralis sizes were similar among treatments, differential expression among these genes could reflect a potential lag between gene expression signaling and ensuing changes in muscle growth. Alternatively, differential levels of muscle repair among temperature treatments could have resulted from varying rates of either acute injury to the muscle or protein metabolism. Indeed, prolonged periods of intense shivering are associated with elevated levels of muscle damage in some small birds, although not in juncos (Swanson and Thomas, 2007). Enhanced tissue carbon turnover during cold exposure has also been documented in birds (Bauchinger et al., 2010) and higher protein metabolism in the cold could elevate expression of muscle growth or repair elements.

Very few genes (6%) in the mTOR signaling pathway were differentially regulated among temperature treatments, and those that were differentially expressed were counterintuitively upregulated in warm-acclimated birds. Upstream regulators of the mTOR pathway, myostatin (MSTN) and IGF1, are central in regulating avian embryonic muscle development (e.g. Guernec et al., 2003), but their role in seasonal regulation of avian adult muscle size is unclear. Swanson et al. (2009, 2014b) found that increases in pectoralis size of wintering passerine birds were associated with decreased expression of MSTN or its metalloproteinase activators (TLL-1, TLL-2), but increases in flight muscle size of photo-stimulated white-throated sparrows (Zonotrichia albicollis) were conflictingly accompanied by increased expression of both IGF1 and MSTN (Price et al., 2011). Although IGF1 expression was not detected in our dataset, MSTN expression did not differ significantly among either temperature or photoperiod treatments. Taken as a whole, this indicates that cold-acclimated individuals were not using any of the major muscle growth elements to enhance thermogenic performance.

We found some evidence to indicate that juncos may be enhancing oxygen delivery to the muscle via angiogenesis (Table 3). Although the temperature treatments were not associated with differential expression of either hypoxia inducible factor, and were only associated with differential expression of a single gene (paxillin, PXN) across the VEGF pathway, VEGFC was upregulated in cold-acclimated individuals (included in M_sum-associated module T28). VEGFC is a VEGF family member that has also been shown to increase angiogenesis in avian embryonic development (Eichmann et al., 1998). Similar function in cold-acclimated individuals could lead to decreased oxygen diffusion distance at the muscles, enabling increased oxygen delivery for use in oxidative phosphorylation and ultimately aerobic thermogenesis.

In addition to enhancements of O₂ delivery, increased thermogenic performance could also result from regulatory changes that enhance the supplies of metabolic fuels to thermogenic organs. Price (2010) provides a thorough review of the avian fatty acid supply pathway (summarized in Fig. 1B). We targeted genes in this integrated pathway to determine whether seasonal increases in thermogenic capacity could, in part, be driven by increased fuel supplies to the pectoralis mitochondria. Juncos seasonally increase muscle lipoprotein lipase during migration (Ramenofsky, 1990), but its role in enhancing circulatory lipid transport capacity during winter is not clear (Swanson, 2010). We did not find evidence of differential regulation in lipoprotein lipases in response to our acclimation treatments, suggesting that these lipases did not contribute to elevated thermogenic performance. Similarly, increases in the expression of sarcolemmal fatty acid translocase (FAT) in migrating white-throated sparrows suggest that these could also mediate seasonal variation in fatty acid transport capacities (McFarlan, 2007), but FAT was not differentially expressed in juncos across treatments. Seasonal increases in the abundance of cytosolic fatty acid binding protein (FABP) have been documented in the pectoralis of a number of avian species during winter (e.g. Liknes et al., 2014) and cytosolic FABP was upregulated in the pectoralis of cold-acclimated individuals, suggesting that it plays an essential role in fatty acid muscular uptake from the circulation in juncos. Finally, seasonal increases in the enzymes facilitating fat transport across the mitochondrial membrane have also been documented in migrating sandpipers (Guglielmo et al., 2002) but their roles in avian metabolic performance for winter birds are unknown (Swanson, 2010). Although we did not find differential expression of any of the genes encoding the three enzymes that facilitate transport across the mitochondrial membrane (CPT1, CPT2 and carnitine-acylcarnitine translocase), we did find a significant trend towards increased CPT enzyme activity in cold-acclimated individuals. Collectively, these results indicate that increased junco thermogenic performance was not enabled by wholesale changes in gene expression across the lipid fuel supply pathway, but upregulation of certain steps (e.g. cytosolic FABP) with cold acclimation suggests that these steps may be crucial modulators of metabolic flexibility.

Regulatory changes to genes involved in lipid oxidation and oxidative phosphorylation are known to be important for thermogenic capacity in mammals (Cheviron et al., 2012,, 2014). We found similar patterns in juncos, although our results emphasize the importance of the former, suggesting that regulatory changes that enhance thermogenic performance in juncos are made to upstream portions of these integrated pathways. Module T32 (negatively correlated with ΔM_sum, Table 2) included one subunit of NADH dehydrogenase (complex I, NDUFB9), which was downregulated in cold-acclimated individuals, but this is the only member of the oxidative phosphorylation pathway to exhibit temperature-sensitive differential expression. It is unclear how differential regulation of a single component of NADH dehydrogenase – a large complex composed of 41 avian subunits – could influence the function of complex I and subsequently the electron transport chain as a whole. Increased metabolic performance in cold-acclimated juncos was instead associated with heightened expression of genes related to the β-oxidation of fatty acids. For instance, temperature-sensitive transcriptional module T31 was enriched for genes involved in fatty-acyl-CoA biosynthesis and HOAD activity, whereas the expression of T33 was negatively correlated with thermogenic capacity and included terms related to fatty acid elongation. Roughly one-third of the genes involved in the fatty acid metabolism pathway were significantly differentially expressed among temperature treatments. In combination with trends in CPT activity, temperature-sensitive expression patterns in the genes involved in fatty acid oxidation are largely indicative of enhanced lipid oxidative capacity in cold-acclimated individuals.

Cold-acclimated birds also downregulated 25% of the genes involved in glycolysis, indicating a preference for lipid over carbohydrate fuels to enhance thermogenic capacity. This preference has similarly been shown in birds and rodents under conditions of prolonged shivering (Vaillancourt et al., 2005,, 2009). Recent in vitro studies of avian pectoralis mitochondria suggest that selection of fatty acids over glycolytic fuels is due to suppression of pyruvate oxidation rather than relative lipid catalytic potentials (Kuzmiak-Glancy and Willis, 2014). Although we found temperature-sensitive expression patterns in 29% of the genes in the pyruvate metabolism pathway, cold-acclimated birds downregulated the expression of four of these genes while upregulating the expression of the remaining two genes. Thus, these regulatory changes were not wholly indicative of pyruvate oxidation suppression in cold-acclimated individuals.

For a number of birds, winter increases in M_sum are associated with increased cellular metabolic and lipid oxidation capacities, as inferred from metabolic enzyme activities of the pectoralis (e.g. Liknes and Swanson, 2011). Of the three muscle metabolic enzyme activities that we measured, enhanced capacity for aerobic thermogenesis was only associated with elevated activities of CPT. CPT was also the only enzyme to reveal temperature-sensitive activity patterns, whereas CS and HOAD activities exhibited photoperiod-sensitive patterns. Enzyme activities for CS and HOAD were also correlated with transcript abundances of the enzyme-encoding genes, indicating that transcriptomic changes translated into functional effects at the enzyme level.

Nonetheless, a lack of correspondence between CPT and corresponding transcript abundances could be an artifact of the complex relationship between mRNA and protein abundances. This includes factors such as post-transcriptional regulation, targeted protein degradation and differences in half-lives among mRNA and proteins, all of which can obscure the true association between these two entities (Suarez and Moyes, 2012). The correlation between activity and transcript abundance for CPT was further complicated by two factors. First, CPT assays may include both CPT1 and CPT2 enzymes, although their relative contributions to our measurements are not known. We think that we primarily measured CPT1 given that the assay we used followed the production of free CoA after the addition of carnitine as substrate, which is the in vivo role of CPT1, and did not use detergent, usually required to solubilize CPT2. We did not find a relationship between CPT activity and transcript abundance for either gene separately (Fig. 4C) or in combination (not shown). Second, the CPT1b paralog (muscle-specific) is, as of yet, unidentified in the zebra finch (Taeniopygia guttata) genome, thus we had little power to detect it within our dataset. Expression of the CPT1a paralog (liver-specific) was quite low; therefore, if our enzyme assays were indeed biased towards CPT1 activity, we may have missed variation at the transcript level. This lack of stoichiometry is then not surprising, especially in combination with the potential action of other mechanisms, such as post-transcriptional modification.

Enzyme activity patterns also largely paralleled concomitant changes in expression across the corresponding metabolic pathways. Increased thermogenic capacity among populations of deer mice is associated with increased HOAD and cytochrome c oxidase activities, as well as concerted expression patterns across both the fatty acid β-oxidation and oxidative phosphorylation pathways (Cheviron et al., 2012,, 2014). Similarly, mitochondrial transport is considered to be the rate-limiting step in fatty acid oxidation (Holloway et al., 2006), and we found that increased CPT activity in cold-acclimated individuals was accompanied by concerted expression across the fatty acid oxidation pathway in response to temperature. Likewise, in response to photoperiod we found concomitant changes across fatty acid oxidation, the citric acid cycle, and oxidative phosphorylation pathways, and these patterns corresponded with increased CS and HOAD activities among the same treatments. Together, these results indicate differential effects of temperature and photoperiod on the two ends of this integrated pathway: supply is manipulated in response to temperature with seemingly little effect further downstream, whereas photoperiod has less of an effect to supply steps and largely modifies downstream oxidative capacity. In winter, photoperiod and temperature change in unison, suggesting that these concerted environmental cues may help to match supply and demand.

The mechanisms underlying seasonal phenotypic flexibility remain uncertain, but recent advances in genomic technologies have enabled high-throughput characterization of the distinct and concurrent physiological adjustments an organism employs when responding to changes in abiotic stressors. Previous studies that surveyed a limited number of parameters suggest that individuals can increase aerobic performance by making modifications to three broad levels of physiological organization – enhancements to thermogenic organ size, oxygen delivery and cellular aerobic capacity (Fig. 1). We found evidence that juncos are using only a small subset of these potential mechanisms to enhance thermogenic performance in response to experimental manipulations of their abiotic environment. This limited response could be due, in part, to the severity of the stressor that we used and the subsequent strength of the physiological response that was required to offset it. Although the magnitude of increase by cold-acclimated individuals is representative of seasonal changes to M_sum documented in free-living passerines (10–50% from summer to winter; Swanson, 2010), changes in thermogenic performance shown here are likely less dramatic than those juncos sustain in the wild while wintering at higher latitudes. M_sum varies with environmental temperature (Swanson and Olmstead, 1999) and the temperature used for cold-acclimation (3°C) was milder than what juncos probably experience throughout much of their winter range. We suggest that, in the face of maintaining energetically expensive metabolic machinery, juncos may only engage the minimum degree of change that is necessary to meet their specific aerobic demands. Given the modular structure of trait components that influence aerobic performance, efficient regulation of metabolic flexibility may involve differential regulation of trait components according to the relative costs and benefits of particular adjustments.

Changes in thermogenic performance were elicited within a relatively short time period, with cold-acclimated individuals already showing incremental increases in thermogenic capacity when measured midway through the acclimation period (at 3 weeks; Swanson et al., 2014a), indicating that temperature is a proximal cue that can be integrated to make rapid physiological adjustments. The proximate influence of winter temperature in regulating junco metabolic rate has been shown previously (Swanson and Olmstead, 1999), and is not surprising as temperatures can be highly variable in the temperate zone – in contrast to daylength – such that birds must quickly respond to precipitous drops in temperature that could serve as selective events (e.g. Bumpus, 1898). It therefore follows that these changes would occur at the level of fuel mobilization and oxygen utilization, representing a rapid vehicle for response. In comparison, relatively immediate (and temporary) responses in the oxygen cascade come in the form of accelerated breath rate and enhanced cardiac output, which require no physiological reorganization, whereas changes that increase vasculature and mitochondrial abundance represent more persistent changes that are achieved over longer time periods. We found that temperature-sensitive transcriptional modules associated with increased thermogenic performance showed upregulation of genes related to β-oxidation and oxygen utilization. These modules also contain genes related to muscle development and repair, possibly reflecting longer-term changes in muscle size and structure; however, individuals did not exhibit corresponding pectoralis hypertrophy. Collectively, these results suggest differential regulation of pathway components that are associated with short-term and long-term phenotypic adjustments. It remains to be seen how quickly birds can mount responses to temperature and photoperiod manipulations, and which pathway components respond most rapidly to these changes.

Finally, many previous investigations of avian seasonal metabolic flexibility have focused on the migratory period. Avian migration is an extreme athletic feat associated with dramatic reorganization of multiple physiological systems (Piersma and van Gils, 2011), yet it represents an acute stressor that lasts for a relatively short period of time (of the order of days). This is in stark contrast to the chronic stress of winter, which lasts 6 months or more in the temperate north. Therefore, while it may be possible for birds to concurrently maintain many physiological responses for the duration of a migratory event, it could prove too costly to sustain these same responses for an entire winter season. Further investigations are therefore needed to determine how free-living birds use these mechanisms of metabolic flexibility amidst winter conditions.

Details of the sampling methods have been previously described by Swanson et al. (2014a). Briefly, we captured dark-eyed juncos Junco hyemalis Linnaeus 1758 overwintering near Vermillion, SD, USA (∼43°N, 97°W) in mid-December. We housed birds individually and, after a 2 week adjustment period, subjected individuals to one of four experimental temperature-photoperiod regimes, each lasting 6 weeks. The treatments were 16 h light:8 h dark at 24°C (warm LD, N=10), 8 h light:16 h dark at 24°C (warm SD, N=10), 16 h light:8 h dark at 3°C (cold LD, N=10), or 8 h light:16 h dark at 3°C (cold SD, N=9). To obtain appropriate sample sizes, procedures were performed over the course of two consecutive winter seasons (2011–2012 and 2012–2013). We tested for and found no differences in phenotypes among years and therefore pooled samples for subsequent analysis (supplementary material Table S1).

We measured individual thermogenic performance before and after the treatments. We quantified summit metabolic rate (M_sum; the maximum metabolic rate attained under cold exposure; Swanson, 2010), a proxy for thermogenic capacity, using open-flow respirometry in a heliox atmosphere (21% oxygen/79% helium) with sliding cold exposure (as per Swanson et al., 1996). Following the M_sum trial, we measured cloacal temperature to verify that individuals were hypothermic (body temperature ≤37°C). We identified M_sum as the period of peak oxygen consumption over a 5 min window and report it as ml O₂ min⁻¹. To control for existing differences in individual M_sum before acclimation treatments, we calculated the change in individual M_sum before and after treatments (ΔM_sum) for use in all subsequent analyses. We tested for differences among treatment groups after the acclimation period using a two-way ANOVA on ΔM_sum.

At the end of the acclimation treatment, immediately following the final M_sum measurement, we euthanized juncos via cervical dislocation, measured body mass and rapidly excised pectoralis tissue on ice. We weighed pectoralis tissue to the nearest 0.1 mg, then flash froze it in liquid N₂ for use in enzyme activity assays or stored it in RNAlater and froze it at −60°C for transcriptomic analyses. We collected birds under appropriate state (11-7, 12-2) and federal (MB758442) scientific collecting permits and all procedures were approved by the University of South Dakota Institutional Animal Care and Use Committee (Protocol 79-01-11-14C). To test for differences in pectoralis size among treatment groups after the acclimation period, we first expressed pectoralis mass as a proportion of body mass (M_p) and then used a two-way ANOVA on M_p.

We used massively parallel sequencing of pectoralis muscle transcriptomes (RNA-seq; Wang et al., 2009) to quantify genome-wide patterns of gene expression for each of the 39 individuals. We isolated mRNA from pectoralis using Trizol reagent and prepared cDNA libraries for sequencing using the Illumina TruSeq RNA Sample Preparation Kit following the manufacturer's protocols. Libraries were sequenced as 100 nt single-end reads on the Illumina HiSeq2000 platform. We multiplexed 10 individuals in a single flow cell lane using Illumina index primers. Image analysis and base calling were performed using Illumina pipeline software. The raw sequence reads have been deposited in the NCBI sequence read archive (SRA PRJNA256328).

We performed a series of sequential filtering steps to remove low quality reads and index primer sequences. First, we removed all reads with mean Phred quality scores less than 30. Second, we trimmed low-quality bases (Phred score <30) from these remaining high-quality sequences. Finally, we scanned reads for adaptor sequences, and if detected, they were trimmed from the sequence read. The filtering steps resulted in a final dataset of nearly 743 million sequence reads, with an average of 19 million reads per individual (range, 3.5–63.0 million reads), and the trimming steps resulted in an average read length of 97.7 nt.

We estimated transcript abundance by mapping sequence reads to the Zebra Finch genome (build taeGut3.2.4) using CLC Genomics workbench (v.6.0.4). In total, sequence reads mapped to 15,362 unique genes. However, we excluded genes with less than an average of 5 reads per individual because genes with low count values are typically subject to increased measurement error (Robinson and Smyth, 2007). This filtering step resulted in a final dataset of 10,868 detected genes.

We performed normalization and differential expression analyses using the package edgeR (Robinson et al., 2010) in Program R (ver. 3.0, R Development Core Team, 2013). We first normalized read counts and controlled for differences in total library size (number of total reads) among individuals using the function ‘calcNormFactors’ (Robinson and Oshlack, 2010). Following normalization, we tested for differences in transcript abundance among treatment groups using a generalized linear model approach. We first specified a paired design in order to test the effects of photoperiod treatment while adjusting for baseline differences among temperature treatments. We then estimated model dispersion for each gene separately using the function ‘estimateTagwiseDisp’ (McCarthy et al., 2012). We tested for genes that exhibited significant expression differences among photoperiod treatments using a GLM likelihood ratio test and controlled for multiple tests by enforcing a genome-wide FDR of 0.05 (Benjamini and Hochberg, 1995). We repeated this procedure to test the effects of temperature treatment while adjusting for baseline differences among photoperiod treatments.

By detecting suites of genes that exhibit correlated transcriptional patterns, previously undetected modules of ostensibly co-regulated genes can be defined. These modules may expose biologically meaningful associations between surprising assemblages of genes, and thus their use in subsequent functional analyses can provide insight into the mechanistic underpinnings of complex traits (Ayroles et al., 2009). To identify coordinated expression changes among interacting genes, we assessed the degree of correlation in transcript abundance among genes with significant treatment effects (FDR <0.05) using modulated modularity clustering (MMC; Stone and Ayroles, 2009). This method calculates Pearson correlation coefficients among transcript abundances for all pair-wise combinations and then identifies modules of highly intercorrelated genes (Stone and Ayroles, 2009). We denote modules composed of photoperiod-sensitive genes with P and those of temperature-sensitive genes with T.

For each of the modules, we then summarized overall module expression using a principal component analysis and retained the first principal component axis (PC1), which represented on average 84% of within module variation (49–100%; Table 1). To determine whether suites of DE genes were related to whole-organism thermogenic performance, we performed linear regressions of individual ΔM_sum on module expression (PC1) and identified correlations at a significance level of P≤0.05. We repeated this procedure to identify correlations between module expression and pectoralis size (M_p). For modules that exhibited significant associations with either phenotype, we performed functional enrichment analyses using the CORNA algorithm (Wu and Watson, 2009), applied on the zebra finch GO Analysis web interface by Michael Watson, to identify functional categories that were over-represented among the genes that comprise a given module.

We performed analyses on targeted genes and pathways to test their influences on metabolic flexibility. These included all genes in the major muscle growth (mTOR), angiogenesis (VEGF), fuel supply [glycerolipid metabolism, peroxisome proliferator-activated receptor (PPAR) signaling, pyruvate metabolism] and metabolic pathways (glycolysis, fatty acid oxidation, citric acid cycle, oxidative phosphorylation). We identified these genes using T. guttata KEGG Database pathway maps (Kanehisa and Goto, 2000) downloaded from CytoKEGG (v.0.0.5). We also included selected genes in the IGF1 and HIF-1 signaling pathways for which complete T. guttata maps were not available. We identified the presence of these genes among those differentially expressed (either among photoperiod or temperature treatments) and noted fold changes with respect to the environmental stimuli.

Because phenotypic responses could manifest from large-scale changes in expression across the entire pathway, we tested for concerted expression across these same pathways in response to both photoperiod and temperature treatments. To do this, we mapped expression values onto T. guttata pathway maps and counted the number of genes exhibiting positive and negative log fold-change values. We performed binomial tests to test for concerted expression across each pathway with the null expectation that an equal number of genes would be up- and down-regulated in given pathway simply by chance. We identified concerted expression at a significance level of P≤0.05.

Variation in transcript abundance does not consistently produce a proportional response in corresponding protein abundance (Abreu et al., 2009). To investigate the relationship between the junco transcriptome and its corresponding gene products, we quantified the activities of three metabolic enzymes that serve as biomarkers for the flux capacity of important metabolic pathways. Enzyme activity (V_max) is a product of both the concentration and kinetic efficiency of an enzyme, and consequently is not strictly equivalent to protein abundance. Nonetheless, because we expect any potential variation in kinetic efficiency to be randomly distributed among individuals, we assumed that differences in enzyme activities largely reflected differences in enzyme concentration, which in turn should be related to differences in the transcript abundance of enzyme-encoding genes. The enzymes we surveyed included: carnitine palmitoyl transferase (CPT), an indicator of fatty acid transport capacity across the mitochondrial membrane; β-hydroxyacyl coenzyme-A dehydrogenase (HOAD), an indicator of fatty acid oxidation capacity and citrate synthase (CS), an indicator of mitochondrial abundance and cellular metabolic intensity. We utilized pectoralis tissue to quantify enzymes using standard assays conducted at 39±2°C with a Beckman DU 7400 spectrophotometer. Further details have been previously described in Swanson et al. (2014a). We report all enzyme activities as mean mass specific activity (μmoles min⁻¹ g⁻¹). To test whether differences in transcript abundance translated into variation in enzyme activity, we performed linear regressions of individual metabolic enzyme activity on the corresponding transcript abundances for each of the enzyme-encoding genes. This included a single gene for CS, three genes (EHHADH, HADHA, HADHB) representing each of the HOAD subunits, and – because our measures of CPT activity may include both CPT1 and CPT2 enzymes – two genes (CPT1a and CPT2) representing two separate CPT enzymes. We identified associations between enzyme activity and transcript abundance at a significance level of P≤0.05 and identified potential outliers using Grubbs' test after verifying normality (Grubbs, 1950). For measures of CS activity, we first log-transformed the data, which resulted in a normal distribution, before applying Grubbs’ test.

We are very appreciative of Jennifer Jones, for her patience and guidance in the laboratory, and Jeff Brawn, Nathan Senner, Nicholas Sly and Cory Suski for comments made to an earlier version of the manuscript. We also thank Marisa King, Jin-Song Liu, Christopher Merkord and Yufeng Zhang for their contributions quantifying metabolic rates and enzyme activities; Ming Liu, Kyle Kirby, Stephanie Owens, Tyler Tordsen, Jake Johnson and Travis Carter for technical assistance in the laboratory and field; Sol Redlin and Kathie Rasmussen for access to their property for collection of juncos; and Brian Allan for use of his tissuelyzer.