# Gene expression shifts in emperor penguin adaptation to the extreme antarctic environment

Josephine R. Paris, Flávia Akemi Nitta Fernandes, Federica Pirri, Samuele Greco, Marco Gerdol, Alberto Pallavicini, Marine Benoiste, Clément Cornec, Lorenzo Zane, Brian Haas, et al.

► To cite this version:

Josephine R. Paris, Flávia Akemi Nitta Fernandes, Federica Pirri, Samuele Greco, Marco Gerdol, et al.. Gene expression shifts in emperor penguin adaptation to the extreme antarctic environment. Molecular Ecology, 2025, 34 (18), pp.e17552. 10.1111/mec.17552. hal-04787153

HAL Id: hal-04787153

https://hal.science/hal-04787153v1

Submitted on 17 Nov 2024

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ORIGINAL ARTICLE

OPEN ACCESS

# Gene Expression Shifts in Emperor Penguin Adaptation to the Extreme Antarctic Environment

Josephine R. Paris<fnref n="1" /> | Flávia A. Nitta Fernandes<fnref n="1" />,² | Federica Pirri<fnref n="1" />,³ | Samuele Greco⁴ | Marco Gerdol⁴ |
Alberto Pallavicini⁴ | Marine Benoiste² | Clément Cornec²,⁵ | Lorenzo Zane³ | Brian Haas⁶ | Céline Le Bohec²,⁷,⁸ |
Emiliano Trucchi<fnref n="1" />

**Correspondence:** Céline Le Bohec (celine.lebohec@iphc.cnrs.fr) | Emiliano Trucchi (e.trucchi@univpm.it)

**Received:** 29 November 2023 | **Revised:** 17 September 2024 | **Accepted:** 26 September 2024

**Handling Editor:** Sean Rogers

**Funding:** This study was supported by the Institut Polaire Français Paul-Emile Victor (IPEV) within the framework of the Project 137-ANTAVIA (PI: CLB), by the Centre Scientifique de Monaco with additional support from the LIA-647 and RTPI-NUTRESS (CSM/CNRS-UNISTRA), by the Centre National de la Recherche Scientifique (CNRS) through the long-term Studies in Ecology and Evolution (SEE-Life) and Zone Atelier Antarctique et Subantarctique (ZATA) programs. The study was approved by the French ethics committee (last: APAFIS#29338-2020070210516365) and the French Polar Environmental Committee, and permits to handle animals and access breeding sites were delivered by the "Terres Australes et Antarctiques Françaises" (TAAF). JRP is supported by funding from the European Union's Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101068395-"Poly2Adapt." SG was supported by PNRA16_00099 and PNRA16_00234. ET was supported by PNRA_16 00164 ("Programma Nazionale di Ricerca in Antartide." Bando PNRA 5 aprile 2016, n. 651. - Linea B "Genomica degli adattamenti estremi alla vita in Antartide"). Bioinformatic analyses were performed on the HPC clusters at the Department of Life and Environmental Sciences ("HappyComputing@DiSVA"), Marche Polytechnic University, and the Department of Life Sciences and Biotechnology, University of Ferrara.

**Keywords:** birds | cold adaptation | rapid adaptation | speciation | transcriptomics

## ABSTRACT

Gene expression can accelerate ecological divergence by rapidly tweaking the response of an organism to novel environments, with more divergent environments exerting stronger selection and supposedly, requiring faster adaptive responses. Organisms adapted to extreme environments provide ideal systems to test this hypothesis, particularly when compared to related species with milder ecological niches. The Emperor penguin (*Aptenodytes forsteri*) is the only endothermic vertebrate breeding in the harsh Antarctic winter, in stark contrast with the less cold-adapted sister species, the King penguin (*A. patagonicus*). Assembling the first *de novo* transcriptomes and analysing multi-tissue (brain, kidney, liver, muscle, skin) RNA-Seq data from natural populations of both species, we quantified the shifts in tissue-enhanced genes, co-expression gene networks, and differentially expressed genes characterising Emperor penguin adaptation to the extreme Antarctic. Our analyses revealed the crucial role played by muscle and liver in temperature homeostasis, fasting, and whole-body energy metabolism (glucose/insulin regulation, lipid metabolism, fatty acid beta-oxidation, and blood coagulation). Repatterning at the regulatory level appears as more important in the brain of the Emperor penguin, showing the lowest signature of differential gene expression, but the largest co-expression gene network shift. Nevertheless, over-expressed genes related to mTOR signalling in the brain and the liver support their central

<footnote n="1">¹Department of Life and Environmental Sciences, Marche Polytechnic University, Ancona, Italy | ²Université de Strasbourg, CNRS, IPHC UMR 7178, Strasbourg, France | ³Department of Biology, University of Padova, Padova, Italy | ⁴Department of Life Sciences, University of Trieste, Trieste, Italy | ⁵ENES Bioacoustics Research Laboratory, CRNL, CNRS, Inserm, University of Lyon, Saint-Etienne, France | ⁶Broad Institute of MIT and Harvard, Cambridge, Massachusetts, USA | ⁷CEFE, Université de Montpellier, CNRS, EPHE, IRD, Montpellier, France | ⁸Département de Biologie Polaire, Centre Scientifique de Monaco, Monaco, Monaco</footnote>

<footnote>Josephine R. Paris and Flávia A. Nitta Fernandes joint first co-authors.</footnote>

<footnote>This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>

<footnote>© 2024 The Author(s). *Molecular Ecology* published by John Wiley & Sons Ltd.</footnote>
role in cold and fasting responses. Besides contributing to understanding the genetics underlying complex traits, like body energy reservoir management, our results provide a first insight into the role of gene expression in adaptation to one of the most extreme environmental conditions endured by an endotherm.

# 1 | Introduction

Organisms adapted to extreme environments provide ideal systems to study the evolutionary mechanisms that drive adaptation and speciation. Extreme environments are characterised by severe abiotic or biotic conditions that are challenging, or lethal, for most organisms to function in (Bell and Callaghan 2012; Rothschild and Mancinelli 2001). From an evolutionary perspective, the strong selective regimes imposed by such environments allow us to probe the fringe of biological resilience and evolutionary innovation, and more generally, the causes and consequences of adaptive evolution (Grant et al. 2017; Hoffmann and Hercus 2000; Hoffmann and Parsons 1997). Comparison of these highly specialised organisms to ancestral lineages, or sister species, occupying milder environments, provides a powerful framework to understand the evolutionary processes that mediate organismal function under adverse environmental conditions (Tobler et al. 2018).

A predictably extreme environment (sensu Wingfield, Kelley, and Angelier 2011), such as a desert, a polar or an alpine region, is one outside an ancestral species' centre of distribution or environmental envelope (Bell and Callaghan 2012), which represents a vacant niche that the ancestral species must colonise (Chapin, Autumn, and Pugnaire 1993; Chesson and Huntly 1997; Jackson et al. 2009; Lekevičius 2009). In such a scenario, we can predict strong selection and rapid adaptation (Chevin and Hoffmann 2017; Grant et al. 2017). Compelling examples include rapid adaptation to temperature extremes in *Drosophila* (Hoffmann, Sørensen, and Loeschcke 2003), and microbial adaptation to UV irradiation, extreme pH, and salinity (Wani et al. 2022). Amongst vertebrates, support for rapid adaptation to extreme environments has been evidenced in sheep following domestication to plateaus and deserts (8–9kya; Yang et al. 2016), and polar bears to the extreme cold of the Arctic (400kya; Liu et al. 2014). Characterising the evolutionary processes contributing to the physiological, metabolic, morphological, and life cycle adaptations to extreme environments could therefore also enhance our predictions for describing the capacity of organisms to respond to the rapidly changing conditions of the Anthropocene (Botero et al. 2015; Chevin, Lande, and Mace 2010).

Rapid adaptive responses to extreme environments can be achieved through plastic or genetic evolutionary changes (or a combination of the two). The former can be accounted for by processes regulating gene expression in response to environmental change (Schlichting and Smith 2002). But how much does the regulation of gene expression favour rapid adaptation to an extreme environment? Gene expression serves as a functional link between the genotype and the cellular and organisomal physiology and can potentially lead to adaptive phenotypes. Expression patterns can evolve rapidly in response to differential selection between environments (Fay and Wittkopp 2008;

Whitehead and Crawford 2006a), where the diversity of cellular functions is achieved through the expression of individual genes orchestrated by large, layered regulatory networks (Davidson and Erwin 2006; Wittkopp 2007). If similar environmental conditions have been encountered previously, such responses may be “hard-wired,” involving the under- or over-expression of pre-selected genes and networks regulating specific biological processes (Jacob and Monod 1961; Romero, Ruvinsky, and Gilad 2012). Gene expression networks can also contain a considerable degree of freedom, allowing an organism to stochastically tune its gene expression and cellular functions to cope with any novel challenges, like those not encountered during its evolutionary history (López-Maury, Marguerat, and Bähler 2008; Tirosh et al. 2006; Xiao et al. 2019). Alongside conventional gene regulation, adaptive reinforcement of initially stochastic gene expression changes may be a key option for organisms to adapt to new and challenging environments (Freddolino et al. 2018; López-Maury, Marguerat, and Bähler 2008). Therefore, if gene expression is indeed important in shaping adaptation to extreme environments (see Hoffmann and Parsons 1997), we expect the signal of gene expression divergence at both individual genes and gene expression networks to be more related to the newly adaptive phenotypes of the focal species than to broader processes governing the biology of both species.

A model that can be used to test these expectations is the case of the Antarctic Emperor penguin (*Aptenodytes forsteri*) and the sub-Antarctic King penguin (*Aptenodytes patagonicus*). Despite having diverged ca. 100,000 generations ago (Bird et al. 2020; Gavryushkina et al. 2017), the two species show clear morphological and behavioural differentiation as a consequence of ecological adaptation to totally divergent temperature regimes. The Emperor penguin completes its breeding cycle in Antarctica (covered by a permanent ice-sheet since 6–15mya; Zachos et al. 2001), at temperatures below –45°C with winds of up to 50m a second, whilst the King penguin breeds in the ice-free and more temperate climates of the sub-Antarctic islands, where temperatures rarely drop below –5°C (Kooyman 1993; Le Maho 1977; Prévost 1961). Fossil evidence, paleoclimatic reconstructions, and ancestral niche construction suggest that whilst all penguins are birds of cold-temperate to subtropical environments, the Emperor penguin represents the furthest evolutionary state of cold adaptation (Pirri et al. 2022; Simpson 1946, 1976; Stonehouse 1969; Vianna et al. 2020). Consequently, as the sister species, the King penguin more closely resembles the ancestral state of *Aptenodytes* penguins (Pirri et al. 2022).

To withstand the harsh environmental conditions of Antarctica, the Emperor penguin hosts a suite of remarkable physiological and behavioural adaptations. Being the largest extant penguin species, with a body mass ranging from 20 to 40kg (the King penguin weighs approximately half this: 10–20kg), its extensive subdermal body fat provides both insulation and energy reserves (Groscolas and Robin 2001; Le Maho, Delclitte, and
Chatonnet 1976). Heat conservation is also highly optimised
via advanced feather and plumule structure and distribution
(Williams, Hagelin, and Kooyman 2015), and counter cur-
rent arterio-venous heat exchangers (Thomas, Ksepka, and
Fordyce 2011). Behavioural adaptations, such as huddling (Ancel
et al. 1997; Gilbert et al. 2006), are also crucial for Emperor pen-
guins to protect themselves against the extreme cold and to lower
energy expenditure. Moreover, Emperor penguin reproduction
is associated with a period of prolonged fasting: up to 45 days for
adult females and 4 months for adult males (Ancel et al. 1992).
These fasting periods are associated with metabolic adaptations
to provide an optimal management of energy resources, includ-
ing the maintenance of stable plasma glucose concentrations,
moderate plasma-free fatty acids, and modest ketosis (Groscolas
and Leloup 1986; Groscolas and Robin 2001). In addition, the re-
productive, moulting, and fasting cycles of the Emperor penguin
have to be completed during the dramatic photoperiodic cycles
of Antarctica (Miché et al. 1991).

Here, we investigate the role of gene expression in adaptation
to the extreme Antarctic environment in the Emperor penguin
by evaluating gene expression patterns across five tissues: brain,
kidney, liver, muscle, and skin, collected from 10 wild chicks of
both the Emperor penguin and the King penguin. Our sampling
concerned tissues that we hypothesised would exhibit the major
gene expression differences for cold adaptation: brain (cold tol-
erance, circadian rhythm, huddling behaviour), kidney (osmo-
regulation, excretion of metabolic waste products), liver (lipid
and fatty-acid metabolism), muscle (thermogenesis), and skin
.feather development and thermal insulation). After generating
the first transcriptomes of both species, we asked: (1) What is
the role of species-specific tissue-enhanced transcripts in ad-
aptation? (2) What is the extent of divergence in tissue-specific
networks of genes between ecologically divergent species? (3)
What are the biological functions of tissue-specific differentially
expressed genes in the Emperor penguin?

2 | Materials & Methods

2.1 | Sample Collection

Emperor and King penguin individuals were sampled from natu-
ral populations between June and September 2016. Emperor pen-
guins were sampled from a colony of Pointe Géologie close to the
Dumont d'Urville Station in Adélie Land, Antarctica. King pen-
guins were sampled from the colony La Baie Du Marin, close to
the Alfred Faure Station in Possession Island, Crozet Archipelago
(Figure S1). All procedures employed during fieldwork were ap-
proved by the French Ethical Committee and the French Polar
Environmental Committee, and authorisations to enter the
breeding site and to collect samples from dead birds were granted
by the "Terres Australes et Antactiques Françaises" (permit
nos 2015-52 and 2015-105). Ten chicks (King penguins aged
between 3 and 8 months, Emperor penguins between 1.5 and
2.5 months) were sampled for each species. Due to differences in
developmental time between the two species (i.e., King penguins
take ca. 11 months from hatching to fledging, whilst Emperor
penguins take ca. 5 months; Ancel, Beaulieu, and Gilbert 2013;
Borboroglu 2015), a relationship between age and developmental
stage is very challenging to estimate. To correct for this variation,

we used a high sample size: for each individual, five tissues
(brain, kidney, liver, skin, and muscle) were sampled for a total
of 100 samples (50 samples per species). All tissue samples were
collected from freshly predated chick corpses, which were col-
lected either immediately after observed predation, or less than
1 h after death (Table S1). Although individual variation in time
since death could affect gene expression, major expression pro-
files should have been preserved by the timely collection of tis-
sues (Shao et al. 2023). Tissues were dissected and chopped into
pieces ~4mm, which were directly fixed in RNAlater (Applied
Biosystems, UK) and stored at -80°C.

## 2.2 | RNA Extraction, Library Construction, and Sequencing

For brain, kidney and liver tissues, total RNA was isolated using a
standard TRIzol-chloroform protocol. RNA extraction from skin
and muscle tissue using this method resulted in poor RNA yields.
Total RNA for these tissues was therefore extracted using the
RNeasy Fibrous Tissue Mini Kit (Qiagen, Germany) according to
the manufacturer's instructions. RNA purity and concentration
were assessed using a Nanodrop 2000 (Thermo Fisher Scientific,
CA) and a Qubit 4.0 fluorometer (Thermo Fisher Scientific, CA),
respectively. RNA integrity (RIN) was evaluated by UV transil-
luminator and Agilent 2100 Bioanalyzer (Agilent technologies,
CA), accepting a RIN score > 6 for library preparation.

For the construction of sequencing libraries aimed at de novo
transcriptome assembly, the total RNA samples extracted from
five tissues (brain, kidney, liver, muscle and skin) in three in-
dividuals of each species, were pooled in equimolar concen-
trations. This resulted in a total of three RNA pools for each
species. Library preparation and sequencing were carried out by
BMR Genomics Service (Padova, Italy). Libraries were synthe-
sised using the Illumina TruSeq Stranded mRNA Sample Prep
kit. Poly-A mRNA was fragmented for 3 min at 94°C, and each
purification step was carried out with 1×Agencourt AMPure
XP beads. Paired-end sequencing (100 bp × 2) was performed on
an Illumina NovaSeq 6000 at a targeted sequencing depth of 20
million reads per library.

For the analysis of gene expression, library preparation was per-
formed using the QuantSeq 3'mRNA-Seq Library Prep Kit FWD
V1 (Lexogen, Vienna, Austria), using the total RNA extracted
from the five tissues. QuantSeq provides an accurate method for
gene expression (even at low read depths) by generating only one
fragment per transcript, making the number of reads mapping
to a gene proportional to its expression (Corley et al. 2019; Moll
et al. 2014). Individual libraries were barcoded and pooled for
sequencing on an Illumina NextSeq500 (single-end, 75 bp) at a
target sequencing depth of 5 million reads per library.

2.3 | De Novo Transcriptome Assembly
and Annotation

Read quality was examined using FastQC v0.11.9 (Andrews 2010).
Trimming was performed using fastp v0.20.1 (Chen et al. 2018),
enabling the elimination of homopolymer tails. Sequences
<71 bp were discarded. Transcriptome assembly for each species
was performed using the Oyster River Protocol (ORP) v2.2.5 (MacManes 2018). By taking advantage of the strengths of different assembly tools, the ORP method recovers transcripts that may be missed by individual assemblers and minimises redundancy, improving overall assembly quality. Within ORP, trimmed reads were error-corrected using Rcorrector v1.0.4 (Song and Florea 2015). Reads were then assembled using three different *de novo* assemblers and different k-mer sizes: Trinity v2.11.0 with 25-mers (Grabherr et al. 2011; Haas et al. 2013), rnaSPAdes v3.14.1 with both 55-mers and 75-mers (Bushmanova et al. 2019), and Shannon v0.0.2 with 75-mers (Kannan et al. 2016). Contigs that were expressed at <1 transcript per million (TPM) reads were removed. This resulted in four distinct assemblies that were merged and clustered into isoform groups using OrthoFuse within ORP. The resulting assembly was used as input to TransRate v1.0.3 (Smith-Unna et al. 2016), and the best transcript (i.e., highest contig score) from each group was placed in a new assembly file to represent the entire group. Only transcripts >250 bp were retained. We removed contigs derived from ribosomal and mtDNA from the assembly (Text S1).

The overall quality of the assemblies was evaluated with BUSCO v5.2.2 (Simão et al. 2015), using the Aves (n = 8782) OrthoDB v10 database (Kriventseva et al. 2019). We also used TransRate v1.0.3 (Smith-Unna et al. 2016) to further assess completeness. The Ex90N50 statistic was quantified as a measure of contiguity, where the N50 value is computed based on the top 90% of transcripts based on expression levels (Haas et al. 2013). Expression levels were estimated by mapping the trimmed reads used for the assembly back to the assembled transcripts using Salmon v1.10.2 (Patro et al. 2017) and statistics were calculated using the TrinityStats.pl script (Haas et al. 2013).

The Trinotate software (github.com/Trinotate/Trinotate/wiki, v4.0.1) was used to annotate the transcriptomes and predict coding regions. First, coding regions were predicted within using TransDecoder (github.com/TransDecoder/TransDecoder/wiki, v5.7.1). The transcript sequences and predicted peptides were used as input to Trinotate for automated functional annotation including identification of protein domains (HMMER v3.1; Finn, Clements, and Eddy 2011) search of Pfam v35 (Finn et al. 2014), top BLASTP and BLASTX match to SwissProt (Release 2023_03) searched using DIAMOND; Buchfink, Reuter, and Drost 2021), signal peptides via SignalP v6.0 (Teufel et al. 2022), transmembrane domains via TMHMM v2.0c (Krogh et al. 2001), non-coding RNA genes using Infernal v1.1.4 (Nawrocki and Eddy 2013), and additional functional annotations captured from eggNOG-mapper v2.1.9 (Cantalapiedra et al. 2021). Long non-coding RNAs (lncRNAs) are unlikely to produce proteins but they may still play a key role in eukaryotic gene regulation (Mattick et al. 2023). We therefore also generated catalogues of lncRNAs for both species (see details in Text S2, available on GitHub—github.com/josieparis/EP-KP-transcriptomics).

### 2.3.1 | Transcript Mapping to de Novo Transcriptomes

Raw QuantSeq reads were cleaned using bbduk in BBmap v39.01 (Bushnell 2014) with parameters specific to 3' QuantSeq data

(Text S3). Quality was assessed using FastQC (Andrews 2010). Trimmed reads from each species were mapped to the respective transcriptome using Salmon v1.10.2 (Patro et al. 2017) to perform transcript-level abundance estimation, setting the mapping parameters for 3' QuantSeq data (Text S3). Raw counts were filtered to remove transcripts with low counts (<10 reads in at least 9 or 10 samples in the Emperor penguin and King penguin groups, respectively—one Emperor penguin sample was removed; see Section 3). This filter reduced our raw table of transcript counts to 21,977 in the Emperor penguin and to 22,793 in the King penguin.

### 2.3.2 | Analysis of Tissue-Enhanced Transcripts (TEGs)

Analysis of tissue-enhanced transcripts was performed using DESeq2 in R (Love, Huber, and Anders 2014), specifying “Tissue” as the negative binomial in the GLM. To quantify tissue-enhancement, contrasts were performed on the counts of each target tissue versus the mean counts of all other tissues. An FDR threshold <0.01 and a log2foldchange ≥5 were used to determine significant tissue-enhanced transcripts (sensu Sjöstedt et al. 2020; proteinatlas.org/humanproteome/tissue/tissue+specific) and ≤-5 as tissue-inhibited transcripts (significantly lower in a particular tissue compared to the average level in all other tissues). To explore the patterns of tissue-enhanced and tissue-inhibited transcripts in each species, we visualised patterns of enhancement or inhibition using the variance stabilising transformation (VST) normalised counts. To identify species-specific genes, transcripts were intersected with the transcriptome annotation to identify putative transcript function (and are therefore referred to as genes). To explore the biological processes of the unique tissue-enhanced genes (TEGs), we performed Gene Ontology (GO) enrichment analysis on the Biological Process (BP) categories using the enricher function in clusterProfiler v4.8.2 (Yu et al. 2012), using a minimum gene size of 10 and a maximum gene size of 500 and an FDR ≤0.05 to test for significance. For both species, genes annotated with Trinotate with associated GO IDs (17,874 in the Emperor penguin and 17,098 in the King penguin) were used as the background.

### 2.3.3 | Transcript Mapping to the Emperor Penguin Genome

For analysis of regulatory networks and differential gene expression, we aligned the Quantseq data of both species to the Emperor penguin genome (GCF_000699145.1; Li et al. 2014). Genome-aligned data were used so that a direct comparison of the same gene sets between the two species could be performed. Note that attempts at reciprocal transcript mapping between the Emperor and King penguin transcriptomes suffered from incomplete isoform resolution. The Emperor penguin genome was chosen because it has improved assembly statistics and a higher number of annotated protein-coding genes in comparison to the King penguin genome (Table S3). Trimmed sequencing reads were aligned using STAR v2.7.9a and counted with htseq-count v0.11.3 with parameters specific to QuantSeq data (Text S3). The raw table of gene counts including samples from both species was filtered to remove
genes with < 10 reads in at least 9 samples of each tissue (re-
ducing the table of gene counts from 15,333 to 10,797 for the
King penguin and Emperor penguin, respectively). Clustering
of the biological replicates of each of the five tissues was as-
sessed using the VST-normalised count data via principal
component analysis (PCA).

2.3.4 | Analysis of Co-Expression Modules

Construction of co-expression modules was performed using
weighted gene co-expression network analysis (WGCNA;
Langfelder and Horvath 2008). Lowly expressed genes are
not informative in WGCNA, and can even generate noise
(Langfelder 2018). It is therefore recommended to reduce the
number of input variables (i.e., number of expressed genes)
prior to analysis. We retained the upper 75th quantile of the nor-
malised count data, resulting in a dataset of 2699 genes. All sam-
ples and genes passed the missing data filter implemented using
the default parameters in the *goodSamplesGenesMS* function.

Gene co-expression networks were constructed for both spe-
cies using the *blockwiseModules* function with default options,
except for the network type, which was set to be signed. In
weighted co-expression networks, genes (network nodes) are
connected by the pairwise correlation of their expression (net-
work edges) in an adjacency matrix. Signed networks exclu-
sively detect positive correlations amongst genes, avoiding the
confounding influence of negative correlations, which were
not of interest in this study. In WGCNA, the adjacency (cor-
relation) matrix is raised to a power $\beta \ge 1$ (soft thresholding),
so that high correlations are accentuated over low correlations
(Zhang and Horvath 2005). We designated a soft threshold
value $\beta$ of 12 for both the Emperor penguin and King penguin
adjacency matrices, corresponding to the value for which the
scale-free topology model fit ($R^2$) stabilised with increasing
power (Figure S2). Co-expression modules for each species
were defined through the topological overlap measure (TOM)
of their respective adjacency matrices, with dissimilarity based
on hierarchical clustering. The resulting gene dendrogram was
trimmed using the dynamic tree-cutting method (Langfelder
and Horvath 2008) and a mergeCutHeight of 0.25, generating
modules with at least 30 genes.

To associate co-expression modules to the five different tissues,
we defined the first principal component (PC1) of each module's
gene expression profile, also known as the module eigengene
(ME) (Langfelder and Horvath 2008). The correlation between
the MEs and the tissue samples was tested using a Pearson
correlation (cor function), with p-values determined using a
Student's asymptotic test (corPvalueStudent function). We ran
five correlation analyses, one per tissue, by setting the samples
from the tissue of interest to one and all other samples to zero.
We considered a module to be positively correlated with a tissue
if cor > 0.5.

Co-expression modules from a reference species that are not pre-
served in a test species can be indicative of regulatory novelties
in the reference group (Oldham, Horvath, and Geschwind 2006).
We therefore tested if co-expression modules defined in the net-
work of the Emperor penguin (reference network) were also

preserved in the King penguin (test network). Gene density and
connectivity statistics were calculated using the *modulePreser-
vation* function, outputting a series of preservation statistics (see
Langfelder et al. 2011 for detailed information). The significance
of each preservation statistic was validated through random per-
mutation tests (Langfelder et al. 2011).

To interpret the potential biological function of the genes
within our chosen modules (i.e., the Emperor penguin-specific
modules, and thus the least preserved modules in the King
penguin network), we performed a GO enrichment analysis.
We also reasoned that modules that are conserved between the
two species should be less related to Emperor penguin adapta-
tion. To this end, we also performed GO term enrichment on
the other modules that we identified as correlated to particular
tissues: blue (liver), black (liver), brown (muscle), red (brain),
green yellow (skin), magenta (skin). For all analyses, we used
the chicken (*Gallus gallus*) gene sets and related GO terms
as a starting point, which were downloaded from Ensembl
using biomart v2.56 in R. Intersecting this database with the
Emperor penguin genome annotation (14,572 protein-coding
genes), we harnessed 66,367 Biological Process (BP) GO terms
related to a total of 9132 genes (i.e., 63% of the genes had asso-
ciated GO terms), which was used as the background for GO
enrichment. Analysis was performed using the enricher func-
tion in clusterProfiler v4.8.2 (Yu et al. 2012), using a minimum
gene size of 10 and a maximum gene size of 500, and an FDR
value ≤ 0.05.

Finally, we identified hub genes for all modules, which corre-
spond to highly connected nodes that can be representative of
that module's main function (Langfelder and Horvath 2008).
We first calculated intramodular connectivity (kIM) for genes
within each module in the Emperor penguin network with the
*intramodularConnectivity,fromExpr* function, using the biweight
midcorrelation (*bicor* function) (Langfelder and Horvath 2012).
We then determined hub genes by selecting nodes with kIMs
≥ 90th quantile of the module's kIM distribution.

2.3.5 | Analysis of Differentially Expressed Genes
(DEGs)

To identify the patterns and biological functions of differen-
tially expressed genes (DEGs) across the five tissues, we per-
formed a standard analysis of differential gene expression
using DESeq2 (Love, Huber, and Anders 2014). We specified
“Group” as the negative binomial in the GLM, where “Group”
was defined as “Tissue+Species.” Contrasts were performed
between the tissue of each species using the Emperor pen-
guin as the baseline. An FDR threshold < 0.01 and a log2fold-
change > 2 or < 2 were used to determine significant DEGs,
where Emperor penguin under-expressed genes (and thus
King penguin over-expressed genes) are indicated by a positive
log2foldchange and Emperor penguin over-expressed genes
(thus King penguin under-expressed genes) are indicated by
a negative log2foldchange. Volcano plots were plotted using
EnhancedVolcano (Blighe, Rana, and Lewis 2018). To limit
our search space for exploring the most likely gene candidates
involved in Emperor penguin adaptation, we specifically ex-
plored the gene annotations of the top ten under-expressed or
over-expressed DEGs (ordered by adjusted *p*-value) for each tissue in the Emperor penguin.

3 | Results

**3.1 | Highly Similar Gene Expression Patterns**
**Between the Emperor Penguin and the King**
**Penguin**

We sampled ten biological replicates from five tissues of
Emperor and King penguin chicks and used 3' QuantSeq to de-
scribe the patterns of gene expression for each tissue. Across
the independent samples, we obtained an average of 6,181,711
(SEM ± 184,196) clean reads for the Emperor penguin and an
average of 8,439,985 (SEM ± 315,902) clean reads for the King
penguin (Table S3). We mapped the multi-tissue samples of
each species to their respective *de novo* transcriptomes (full
details in Text S4, Table S4, Figure S3). Assessment of the
counts showed good consistency across replicates and tissues

(Figure S4). One biological replicate of the brain from the
Emperor penguin (EB04) was removed due to a low number of
recovered reads (Table S3).

Evaluation of the ordination of the samples via PCA revealed high tissue-specific clustering of the biological replicates in both the genome-mapped data (Figure S5) and transcriptome-mapped data (Figure 1A,B). In the transcriptome-mapped data, the samples clustered by tissue in both species, with both PC1 (Emperor penguin: 38%; King penguin: 37%) and PC2 (Emperor penguin: 22%; King penguin: 27%) explaining a similar proportion of the variance. PC1 primarily separated all tissues, except for skin and brain, which were separated by PC2. In the Emperor penguin, the muscle was more different on both axes. Separation of the skin tissue was less apparent in the Emperor penguin data, due to two replicates which deviated from the tissue's centroid. Further inspection of these outlier samples identified no artefactual reasons for the removal of these replicates; therefore, they were retained in downstream analyses.

<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | Tissue-specific patterns of gene expression in the Emperor and the King penguin. (A and B) Principal component analysis (PCA) of the normalised expression of tissues in the (A) Emperor penguin or (B) King penguin. In both PCAs, samples are coloured by tissue as shown in the legend: Brain (green); kidney (purple); liver (pink); muscle (orange); skin (black). (C and D) Variance stabilising transformation (VST) normalised expression counts of the tissue-enhanced genes in the (C) Emperor penguin (brain: 1453, kidney: 501, liver: 646, muscle: 455, skin: 426 tissue-enhanced genes; Table S5) or (D) King penguin (brain: 1372, kidney: 521, liver: 646, muscle: 549, skin: 536 tissue-enhanced genes; Table S5).</figcaption></figure>
Overall, assessment of the tissue-enhanced and tissue-
inhibited genes showed very similar patterns between the
two species (Figure S6, Table S5), with brain, muscle, and
skin tissues showing high tissue specificity in both species
(Figure 1C,D). Kidney-enhanced transcripts were shared
mostly with the liver, and similarly, liver-enhanced tran-
scripts were shared mostly with the kidney. Muscle-enhanced
transcripts showed some shared expression with the skin and
the brain, whereas skin-enhanced transcripts were not shared
in the muscle, but were slightly shared in the kidney. Despite
the overall highly shared landscape of tissue expression in
both species, the numbers only reflect general patterns of tis-
sue regulation.

**3.2 | Identifying Unique Tissue-Specific Genes and Biological Processes in the Emperor Penguin**

Given the high similarity of the patterns of tissue enhance-
ment and tissue inhibition between the species, we next aimed
at identifying the tissue-specific genes that were unique to the
Emperor penguin and the King penguin. For comparative pur-
poses, we used only those transcripts which could be annotated
(referred to as “genes”; see Table S6 for full details). Looking at
all enhanced and inhibited genes (Venn diagrams; Figure 2),
the brain showed the highest percentage of shared genes (45%;
Figure 2A), whilst the liver (31%, Figure 2C) and skin (31%,
Figure 2E) showed the lowest percentage of shared genes.
Compared to the King penguin, the Emperor penguin had more
unique genes in the brain (389) and liver (290).

Inspecting only the unique tissue-enhanced genes (TEGs), the
brain showed a similar proportion between the two species
(Figure 2A), but the King penguin TEGs were associated with a
higher number of enriched GO terms (32, compared to four in
the Emperor penguin; see Table S7 for a full list of significant
GO terms for each tissue). The kidney showed a (slightly) higher
proportion of TEGs in the Emperor penguin, but fewer enriched
GO terms (Figure 2B), and the functions appeared to be unre-
lated to adaptation. In the liver, although the King penguin had
a higher proportion of unique TEGs, these were associated with
only four GO terms. In contrast, Emperor penguin liver TEGs
were associated with 67 GO terms (Figure 2C) related to me-
tabolism—in particular lipid and glucose metabolism, energy
production, and food intake—including digestion, regulation
of feeding behaviour, and response to starvation (Table S7).
The muscle showed the lowest number of total genes unique to
the Emperor penguin (20%), but a slightly higher proportion of
unique TEGs when compared to the King penguin (Figure 2D).
GO analysis of these TEGs identified four enriched terms in
the Emperor penguin, potentially related with adaptation, in-
cluding muscle cell fate commitment, and positive regulation
of skeletal muscle fibre development, whereas no terms were
associated with the King penguin TEG set (Figure 2D). Finally,
although the skin showed a lower proportion of unique TEGs
in the Emperor penguin, these were associated with seven
GO terms (Figure 2E), with potential relevance to adaptation,
including epidermis development, canonical Wnt signalling
pathway involved in positive regulation of epithelial to mesen-
chymal transition, and positive regulation of heparan sulphate
proteoglycan biosynthetic process.

**3.3 | Identifying Regulatory Changes in the Emperor Penguin**

We next investigated differences in the regulatory landscape of gene expression by constructing weighted co-expression networks in WGCNA. The Emperor penguin network was characterised by 1943 genes assigned to 12 modules (i.e., gene clusters, defined by colour names) and 756 genes that were not assigned to any module (“grey module”), whilst the King penguin network contained 2042 genes assigned to 15 modules and 657 genes not assigned to any module (Figure 3A, Table S8). For further analyses, we focus on the relationship between Emperor penguin modules with the different tissues and their preservation in the King penguin network.

In accordance with the strong tissue-specificity enhancement
observed at the transcript-level (Figure 1), we also observed
an association between the module eigengenes (ME) and the
tissues in both species (Figure 3B). More specifically, three
Emperor penguin modules were correlated with the brain (tur-
quoise, red, and brown), four modules were correlated with the
liver (green, black, blue, and pink), five modules with the muscle
(purple, yellow, tan, pink, and brown), and two modules with
the skin (magenta and green yellow) (Figure S7). No modules
were significantly correlated with the kidney.

Computing preservation statistics on the 12 Emperor penguin
modules, we identified six modules as weakly preserved in the
King penguin: green, pink, tan, yellow, purple, and turquoise.
Overall, the purple module (associated with the muscle) was
the least preserved (Table S9, Figure S8). In particular, the
statistics that explained this module's low preservation were
the mean Maximum Adjacency Ratio (meanMAR = 0.149) and
the Proportion of Variance Explained (propVarExpl = 0.146)
by the ME (Figure 3C, Table S10). Low meanMAR values
indicate that genes from an Emperor penguin module have
weak correlations in the King penguin network, even if those
genes have a high number of connections (i.e., high density).
propVarExpl is a density statistic that measures how tightly
correlated are genes from a module in the test dataset based
on the module's eigengene (i.e., expression variance in that
module) (Langfelder et al. 2011). Two other muscle modules,
yellow and tan, the muscle-liver pink module, and the liver
green module also showed low values of meanMAR statistics
(Figure 3C, Table S10). In accordance with low meanMAR val-
ues, we also detected a tendency towards lower intramodular
connectivity (kIM) values (i.e., weaker connections) amongst
genes of the three muscle modules, purple, yellow, and tan,
in the King penguin network in comparison to the Emperor
penguin network (kIM distribution is shifted to the left in the
King penguin; Figure 3D; further details on kIM in Langfelder
and Horvath 2008). Another module that showed signals
of weak preservation in the King penguin dataset was the
brain-related turquoise module, with a propVarExpl = 0.200
(Figure 3C, Table S10).

To investigate the biological meaning of the six candidate
Emperor penguin modules: purple (muscle); yellow (muscle);
tan (muscle); pink (muscle & liver); green (liver); turquoise
(brain), we performed GO term enrichment analysis on all genes
identified within each module (see Table S11 for full results).
Emperor

A

Brain

proportion of unique
tissue-enhanced
genes (TEGs)

Number of enriched
GO Terms

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

King

B

Kidney

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

389
(30%)

587

(45%)

317
(25%)

proportion of unique
tissue-enhanced
genes (TEGs)

Number of enriched
GO Terms

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

C

Liver

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

160
(30%)

178
(33%)

200
(37%)

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

290
(36%)

253

(31%)

268
(33%)

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

D

Muscle

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

161
(20%)

253

(32%)

384
(48%)

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

E

Skin

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

177
(31%)

173

(31%)

217
(38%)

<figure><img src="image_13.png" /><figcaption>FIGURE 2 | Legend on next page</figcaption></figure>
**FIGURE 2** | Quantification of unique tissue-specific genes in the Emperor penguin and King penguin. For each tissue (A: brain; B: kidney; C: liver; D: muscle; E: skin), the Venn diagrams represent the total number of unique genes (both tissue-enhanced and tissue-inhibited) in the Emperor penguin (in blue, left) and the King penguin (in yellow, right). Intersect reflects the number of shared tissue-enhanced genes. Bar plots show the proportion of unique tissue-enhanced genes (TEGs) for each species, and the number of enriched GO BP Terms related to each species' set of TEGs.

The purple module was significantly enriched for 11 GO terms related to muscle development and organisation, as well as heart activity regulation. The two other muscle modules, tan and yellow, were significantly enriched for seven and 11 GO terms, respectively. GO terms of both modules related to embryonic cell organisation and differentiation, whilst the yellow module was also enriched for two GO terms related to stress response and hypoxia. The green liver module had one enhanced GO term related to blood coagulation. The pink (muscle and liver) module had one enhanced GO term related to fatty acid beta-oxidation. The turquoise brain module had six enriched GO terms related to neural cell formation and function.

An additional method to assess the biological meaning of module function is through the detection of hub genes, which are highly connected genes within a module (i.e., nodes with high kIM values), usually representative of the module's main regulatory role (Langfelder and Horvath 2008). By selecting genes using the upper 90th quantile of the module's kIM distribution for the Emperor penguin (Figure 3D), we detected nine hub genes in the purple muscle module, 20 genes in the yellow module, and five genes in the tan module. The pink and green modules contained 11 and 18 hub genes, respectively, and the turquoise brain module had 30 hub genes (Table S12).

Finally, as a null hypothesis, we reasoned that conserved modules should be less related to the biological processes governing Emperor penguin adaptation. To explore this, we also performed GO term enrichment on conserved modules (Table S13). We found one shared GO term (axon guidance) between the turquoise module (brain, not conserved) and the red module (brain, conserved). We also identified one GO term (fatty acid beta-oxidation) shared between the pink module (muscle-liver, not conserved) and the blue module (liver, conserved). Hub genes for all modules can be found in Table S12.

## 3.4 | Differentially Expressed Genes Between the Emperor Penguin and the King Penguin

Finally, using the transcript data mapped to the Emperor penguin genome, we sought to identify the patterns and potential function of differentially expressed genes (DEGs) between the species across the five tissues (Figure 4). We found that the muscle showed more than double the number of DEGs ($n = 786$), where the majority ($n = 598, 78\%$) were over-expressed in the Emperor penguin (Table S14). This was followed by the skin ($n = 300$; 63% over-expressed, 37% under-expressed), liver ($n = 287$; 46% over-expressed, 54% under-expressed), kidney ($n = 183$; 74% over-expressed, 26% under-expressed), and finally the brain ($n = 77, 35\%$ over-expressed, 65% under-expressed). We specifically explored the top ten over-expressed or under-expressed genes in each tissue (Table S15), and performed an extensive literature search on these candidates, considering genes related more intimately to the biology of the Emperor penguin in the discussion.

Assessment of the top candidate genes revealed that some DEGs were identified as consistently over-expressed (seven genes) or under-expressed (3 genes) in more than one tissue (Table 1).

## 4 | Discussion

By performing a comparative multi-tissue profiling of the Emperor penguin and the King penguin, we evaluated the role of gene expression in adaptation to the extreme conditions of Antarctica experienced by the Emperor penguin. We hypothesised that if gene regulation is important, the signal of gene expression divergence should be more related to the specific phenotype of the Emperor penguin, rather than to broader biological processes important for both species. We emphasise that due to sampling restrictions in the wild, our analyses were performed on tissues sampled from chicks, but that our interpretation concerns gene expression related to hardcoded adaptive differences between the species. We generated the first high-quality reference transcriptomes for both species. Enrichment analysis on the annotated tissue-enhanced transcripts unique to the Emperor penguin identified a suite of relevant biological processes related to the metabolic demands imposed by the cold and fasting, particularly in the liver and muscle. We next assessed the regulatory landscapes of both species, classifying six Emperor penguin co-expressed modules associated with the muscle, liver and brain which were weakly preserved in the King penguin. Through analysis of differential gene expression between the two species, we found that the muscle showed more than double the number of differentially expressed genes compared to all other tissues. By exploring the genes and processes identified across the analyses, we are able to pinpoint the tissues and potential candidate genes responsible for gene expression-mediated adaptive processes in the Emperor penguin. More broadly, we demonstrate that quantifying patterns of gene expression regulation across multiple tissues and at multiple levels (i.e., individual transcripts, individual genes, co-expression networks) highlights the complex role of adaptive gene expression in ecological divergence and adaptation to extreme environments.

## 4.1 | Brain

As a highly specialised and complex tissue, the brain often shows unique patterns of gene expression in comparative studies (Bentz et al. 2019; GTEx Consortium 2015; Qi et al. 2019) and our analysis of brain-enhanced transcripts also highlighted this tissue as particularly unique (Figure 1C,D). We also observed the lowest annotation rate for the unique brain-enhanced transcripts (Table S6), which is expected given the high specificity of brain genes and an obvious bias in annotation towards model species (Sjöstedt et al. 2020). When comparing gene expression patterns between the species, the brain had the lowest signature of differential gene expression. Yet at the regulatory level, a lower degree of conservation was observed, as the largest co-expression module in the
<figure><img src="image_14.png" /><figcaption>FIGURE 3 | Legend on next page</figcaption></figure>
**FIGURE 3** | Dendrogram and modules of interest in the Emperor penguin co-expression network. (A) Dendrogram of genes clustered based on topological overlap between nodes (i.e., genes), and their assigned modules in the Emperor penguin network (top colours). King penguin modules defined in the King penguin network are represented in the bottom colours to illustrate module assignment differences between the two networks, whilst all other plots (B–D) are based on the Emperor penguin network modules. (B) Module eigengenes (ME, corresponding to the PC1 of the module's expression profile) of weakly preserved Emperor penguin modules in the King penguin network according to the Z permutation and medianRank statistics (Figure S8, Table S9), where B stands for brain, K for kidney, L for liver, M for muscle, and S for skin samples. Filled and empty dots represent Emperor penguin and King penguin samples, respectively. (C) Main summary statistics that indicate the weak preservation of purple, yellow, tan, pink, turquoise, and green Emperor penguin modules in the King penguin network. In the left panel, the mean maximum adjacency ratio (*meanMAR*) statistics, in which low *meanMAR* indicates modules composed of nodes with weak connections (i.e., co-expression) with neighbour nodes in the King penguin dataset. In the right panel, the Proportion of Variance Explained by the module eigengenes (*propVarExpl*), in which low *propVarExpl* indicates low module density in the King penguin dataset, as correlation between gene expression patterns is lower. (D) kIM densities of genes of the least preserved modules (green, pink, tan, yellow, purple, and turquoise) in Emperor penguin (top) and King penguin (bottom) networks. Dashed vertical lines in the Emperor penguin plots indicate the 90th quantile of kIM distribution threshold used to determine hub genes (Table S12).

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

**FIGURE 4** | Volcano plots of the 10,797 genes included in the differential expression analyses between Emperor penguin and King penguin per tissue. Coloured dots correspond to DEGs with an FDR > 0.01 and log₂foldchanges ≤ 2 (blue) and > 2 (yellow), meaning their over-expression in the Emperor penguin and King penguin, respectively. The gene IDs of highly differentiated genes are included.

<table><caption>**TABLE 1** | Differentially expressed genes (DEGs) identified as over- or under-expressed in more than one tissue in the Emperor penguin.</caption><thead><tr><td>Gene symbol</td><td>Tissues</td><td>Expression in Emperor penguin</td></tr></thead><tbody><tr><td>BOK</td><td>Brain, skin, kidney</td><td>Over-expressed</td></tr><tr><td>LOC103895152</td><td>Brain, kidney, liver</td><td>Over-expressed</td></tr><tr><td>LARP1</td><td>Brain, liver</td><td>Over-expressed</td></tr><tr><td>ACP6</td><td>Brain, kidney</td><td>Over-expressed</td></tr><tr><td>TEX33</td><td>Brain, kidney, liver</td><td>Over-expressed</td></tr><tr><td>MLST8</td><td>Kidney, liver</td><td>Over-expressed</td></tr><tr><td>STC2</td><td>Liver, muscle</td><td>Over-expressed</td></tr><tr><td>ARL1</td><td>Brain, kidney, muscle, skin</td><td>Under-expressed</td></tr><tr><td>PPP1R37</td><td>Brain, kidney, skin</td><td>Under-expressed</td></tr><tr><td>GOLPH3L</td><td>Brain, kidney</td><td>Under-expressed</td></tr></tbody></table>

Emperor penguin (turquoise, *n* = 294) was related to the brain and was weakly preserved in the King penguin. This pattern indicates that, despite the high correspondence of the over-expressed and under-expressed sets of genes in both species, the regulatory pathways through which such genes interact in the brain are likely different. This is supported by the fact that our GO analysis of the

red brain module (preserved between the species), contained only five significant terms, which are broad in function (Table S13). However, interpreting differences in regulatory differentiation and their biological relevance to the Emperor penguin ecological (cold, absence of daylight) and physiological (metabolic balance during long fast periods) adaptation is challenging, as the
majority of enriched GO categories in our analyses appeared to be primarily functional (i.e., involved in cell-cell signalling, cell adhesion, and synapse, axon, and neuron processes).

Still, amongst the turquoise module hub genes and the most significant over-expressed genes, we identified a set of potentially relevant candidates (Tables S12 and S15). Related to cold adaptation we identified two genes: SLC6A1, which encodes the GAT1-GABA transporter and has been previously implicated in adaptive cold-stress responses in *Anolis* lizards (Campbell-Staton et al. 2017); and ITPR1, a neuronal intracellular calcium release regulator critical to nociception and pain behaviour (Petrenko et al. 2003), which in particular, has been shown to be regulated under repeated cold stress in rats (Kozaki et al. 2015; Wakatsuki et al. 2021). We also identified genes related to nutrient regulation and metabolism. The most significantly over-expressed gene in the brain, BOK, was also the top significant gene in the skin, and was also highly over-expressed in the kidney (Table 1). BOK is one of the main mitochondrial apoptosis regulators, which is activated in response to stress stimuli, such as DNA damage, oxidative stress, aberrant calcium fluxes, and nutrient deprivation (Lopez and Tait 2015). The known role of this gene in regulating fasting in humans and chickens is suggestive of a similar mechanism in the Emperor penguin (Giménez-Cassina and Danial 2015). LARP1, a key regulator of the mTOR pathway, was the second most over-expressed gene in the brain, and the fifth most over-expressed gene in the liver. The mTOR pathway is regulated in response to growth signals and nutrient availability (Kim et al. 2002) and is involved in metabolic stress adaptation (Wu and Storey 2021). Interestingly, under nutritional limitations, mTOR acts as a metabolic switch to lower global metabolism (Fuentes et al. 2021), suggesting this gene may be important in the regulation of metabolism under fasting conditions. ACP6 (also highly over-expressed in the kidney) regulates mitochondrial lipid biosynthesis by balancing lipid composition within the cell (Li et al. 2013). We also identified MTNR1A as differentially expressed (although not within the top candidates), which mediates melatonin effects on circadian rhythm and reproductive alterations affected by day length (Fishman and Tauber 2023), which could be relevant as Emperor penguin reproduction occurs during the almost constant darkness of austral winter (Miché et al. 1991).

## 4.2 | Kidney

Due to the kidney's role in osmoregulation and renal function, we hypothesised that this tissue might reflect processes related to the high metabolic demands of the Emperor penguin. Indeed, we observed that kidney-enhanced transcripts shared some expression patterns with the liver, potentially reflecting a common landscape of metabolism (especially glucose), which has been previously observed in birds (Braun 2015; Sparr, Braun, and Sweazea 2010; Tinker, Brosnan, and Herzberg 1984). However, aside from the identification of an enriched pathway involved in the regulation of insulin-like growth factor receptor signalling (Table S7), which is directly implicated in renal function under fasting conditions (Feld and Hirschberg 1996; Hirschberg and Kopple 1989), we observed few other biological processes related specifically to Emperor penguin adaptation. Additionally, no co-expressed modules were significantly correlated with this tissue,

which also showed the second lowest number of DEGs. However, amongst the over-expressed genes, we did identify BOK and APC6 (discussed above). RBP4, the third most over-expressed gene, encodes for a specific retinol carrier in the blood, transporting retinoids (vitamin A and its derivatives) to target tissues, such as adipose tissue (Steinhoff, Lass, and Schupp 2022). Previous reports in mice revealed that retinol transport proteins are potent regulators of cold-induced thermogenic responses and cold adaptation by promoting expression of marker genes in the adipocytes (Fenzl et al. 2020; Ribot et al. 2004). Finally, no gene candidates overlapped amongst our analyses for the kidney (Table S16), and the overall concerted signal of the role of this tissue in Emperor penguin adaptation was low.

## 4.3 | Skin

We reasoned that the skin would show a gene expression profile related to thermal insulation and the development of feathers, as the Emperor penguin's thick skin and morphologically specialised plumage is one of the keys to survival in Antarctica (Williams, Hagelin, and Kooyman 2015). These specialised dermic structures can provide 80%-90% of insulation requirements (Le Maho 1977; Le Maho, Delclitte, and Chatonnet 1976), enabling them to maintain a core body temperature of 38°C (Ponganis et al. 2003). In particular, such adaptations should be observed in the early life stages, as Emperor penguin chicks are covered in a unique down that confers a thermoregulatory advantage (Taylor 1986). Aside from the yellow muscle module (which also showed a low degree of correlation with the skin; Figure S7), we did not find any co-expression modules uniquely preserved in the Emperor penguin and overlap amongst analyses was low (Table S16). However, in the analysis of Emperor penguin unique skin-enhanced TEGs, we identified three enriched GO term processes which could be involved in the development of feathers, including epidermis development (Yu et al. 2004), canonical Wnt signalling pathway (involved in positive regulation of epithelial-to-mesenchymal transition; Heller and Fuchs 2015), and heparan sulphate proteoglycan biosynthetic process (related to hair follicle and sebaceous gland morphogenesis; Coulson-Thomas et al. 2014). Amongst the most highly over-expressed genes, we identified BOK (discussed above), ANO1, which acts as an excitatory heat-activated channel in mammalian thermosensation (Cho et al. 2012; Vriens, Nilius, and Voets 2014), a gene (LOC103893312) encoding an adipocyte fatty acid-binding protein, LPL, a key triglyceride metabolism enzyme, involved in cold-stress responses in birds (Bénistant et al. 1998; Radomski and Orme 1971), and PLIN1, a lipid droplet coat protein, which affects body weight and fat deposition (Bickel, Tansey, and Welte 2009). ALPK3 was identified in all analyses (Table S16), and although primarily implicated cardiomyopathy, this gene has been identified as a plasticity candidate in a gene expression analysis of temperature responses in stickleback (Metzger and Schulte 2018). The role of these pathways and genes in feather development, thermosensation, and energy storage are distinctly relatable to the Emperor penguin's ecology. However, overall, the skin had a lower signal compared to the muscle and liver, which may be due to its more highly shared thermoregulatory function in both Emperor and King penguins (Lewden et al. 2017), especially in chicks (Duchamp, Rouanet, and Barré 2002).
Since muscle contraction is coupled to heat production, the muscle is the primary thermogenic organ in most vertebrates (Periasamy, Herrera, and Reis 2017). We found that the unique expression profile in this tissue was particularly conspicuous, encompassing a high proportion of tissue-enhanced transcripts related to enriched GO processes relevant to cold adaptation (Figure 2), four of the six significantly divergent modules in the co-expression analysis (Figure 3), and the largest fraction of differentially expressed genes (Figure 4; 76% of which were over-expressed in the Emperor penguin). The muscle was also the only tissue showing a high amount of overlapping candidate genes across our analyses (Table S16). Although the muscle is a highly specialised tissue, the signatures of gene expression were not limited to a narrow set of genes or processes involved only in the structural organisation of muscle fibres. Instead, our findings are consistent with the crucial role played by the muscle in temperature homeostasis and whole-body energy metabolism (Swanson, Zhang, and Jimenez 2022).

Skeletal muscle produces heat primarily through shivering thermogenesis (Periasamy, Herrera, and Reis 2017), which involves rapid, repeated muscle contractions (Nakamura and Morrison 2011). Accordingly, we identified muscle contraction and striated muscle contraction as enriched GO processes in the least preserved co-expression module (purple) and in uniquely enhanced genes in the Emperor penguin. Although these could also be related to locomotion, locomotor activity in chicks is very limited as they are brooded on the adults' feet, and activity only increases prior to fledging (Corbel et al. 2009). Moreover, relative to total body mass, the skeletal muscle of birds is commonly much larger than that of reptiles and mammals, and it has been hypothesised that this enlargement allows for more effective heat generation, especially in breast and thigh muscles (Newman 2011; Rowland, Bal, and Periasamy 2015). We detected positive regulation of skeletal muscle fibre development as an enriched GO process in the analysis of unique genes and skeletal muscle tissue development was identified as the most significantly enhanced process in the purple module. Moreover, the most central gene in the purple module was VWDE (also present as an Emperor penguin unique TEG), which in humans, is associated with breast hypertrophy (Rappaport et al. 2013). Amongst the most highly over-expressed genes, we also identified FN1, a glycoprotein that links fibrillar matrix tissue, providing strength and elasticity to muscle fibre (Smith, Meyer, and Lieber 2013), and METRNL, which is induced in skeletal muscle upon exercise and during exposure to cold (Bae et al. 2018; Rao et al. 2014). In mouse models, METRNL increased in vivo in exercised mice, suggesting that METRNL is secreted during repeated muscle contraction, and injected METRNL improved glucose tolerance in mice with high-fat-diet-induced obesity or diabetes (Lee et al. 2020). We also identified the gene EYA1 as overlapping in all analyses, which has been shown to be involved in the metabolic and contractile specialisation of mature muscle fibres (Grifone et al. 2007). Moreover, EYA1 proteins accumulate preferentially in the nuclei of fast-twitch glycolytic muscles (Grifone et al. 2004), which have been shown to play a prominent role under acute cold exposure in rats (5°C; Wang et al. 2003).

Whilst shivering is considered the dominant mechanism of heat production in birds (Bicudo, Vianna, and Chaui-Berlinck 2001),

recent empirical evidence also argues for the role of non-shivering thermogenesis (NST; Teulier et al. 2010). In mammals, NST is regulated by the activity of uncoupling protein 1 (UCP1) in brown adipose tissue (Klingenspor 2003). Although birds lack brown adipose tissue and mammalian UCP1 (Mezentseva, Kumaratilake, and Newman 2008), evidence of a homologous functional avian UCP accompanied by empirical evidence of whole-body metabolic heat production unrelated to shivering suggests NST represents a plausible additional mechanism of thermogenesis in birds (Hohtola 2002; Raimbault et al. 2001; Teulier et al. 2014). Although Emperor penguins lack empirical evidence of NST, data from the King penguin suggests NST occurs via concerted mechanisms in the skeletal muscle (Duchamp et al. 1991; Duchamp and Barré 1993; Talbot et al. 2004; Teulier et al. 2010). The first mechanism includes increased activity of avian UCP, which is dependent on superoxide, a reactive oxygen species (ROS). Compared to other Antarctic birds, Emperor penguins exhibit higher ROS scavenging capacity, suggesting they are naturally exposed to higher basal prooxidant pressures (Corsolini et al. 2001). Related to this, we identified two enriched GO processes in the yellow (muscle) module related to the cellular response to ROS. Moreover, the negative regulation of S-nitrosylation, which is a selective protein modification mediated by ROS (Stomberski, Hess, and Stamler 2019) was identified in the purple module, and has been previously implicated in mammalian NST (Sebag et al. 2021). Amongst the top over-expressed genes in the Emperor penguin, we also identified IGF2BP1, an insulin-like growth factor which stabilises mRNA molecules that are recruited to stress granules under oxidative stress and heat shock (Bley et al. 2015), which is specifically induced in chickens under high-fat diets (Chen et al. 2019). The second potential mechanism of NST is increased activity of the adenine nucleotide translocase, which is dependent on the oxidation of fatty acids. The only enriched GO process in the pink module (related to both the muscle and the liver) was fatty acid beta-oxidation, but we note that this term was also identified in the conserved blue module (liver). Whilst these pathways are indicative of a potential role of NST in Emperor penguin gene expression adaptation to the cold, such responses should be confirmed by physiology studies.

Adipose tissue is also a fundamental regulator of energy homeostasis (Zhuo et al. 2015) and has an important role in metabolic energy storage (Choe et al. 2016). In Emperor penguins, fat stored in adipose tissue is by far the major energy reservoir (Groscolas and Robin 2001), with body fat at the beginning of the breeding fast representing approximately 30% of the body mass and 80% of body energy (Groscolas 1990). Our decision to not sample this tissue was due to evidence of the higher importance of skeletal muscle for NST in birds, in addition to the low nucleic acid content of adipocytes (Cirera 2013; Sinitsky et al. 2018). However, we still identified recurrent signals of genes related to adipogenesis in our analyses, including several significant genes identified as adipocyte fatty acid-binding proteins and two genes (MEOX1 and POSTN), which were uncovered as overlapping in all three analyses of the muscle. The homeobox protein MEOX1 is a mesodermal transcription factor involved in somitogenesis and the commitment of cells to the skeletal muscle lineage (Petropoulos et al. 2004; Wu et al. 2018). MEOX1 expressing progenitors in the somite revealed that both brown and white adipocytes arise from the somitic mesoderm
(Shamsi, Wang, and Tseng 2021). POSTN (also identified as a hub gene in the yellow module) is a secreted extracellular matrix protein that functions in tissue development and regeneration and is produced and secreted by adipose tissue. A transcriptionic analysis of porcine adipocytes showed that POSTN was a crucial candidate gene associated with adipogenesis during the differentiation of intramuscular adipocytes (Mo et al. 2017). In mice, the adaptation of adipose tissue to high-fat diet feeding is impaired in animals with systemic ablation of POSTN, suggesting that a loss of periostin attenuates lipid metabolism in adipose tissue (Graja et al. 2018).

## 4.5 | Liver

With the liver's key role in digestion and metabolism, specifically in regulating the production, storage, and release of lipids, carbohydrates and proteins (Zaefarian et al. 2019), we hypothesised gene expression differences in this tissue to be related to the strong metabolic demands typified by the life cycle of the Emperor penguin. In line with this, the liver showed marked differences between the two species, associating to more than 60 enriched GO terms in the analysis of unique TEGs, showing an association to two co-expression modules (green and pink—also correlated to muscle tissue), which were weakly conserved in the King penguin, and a high number of differentially expressed genes (n=287). Fasting in the Emperor penguin is characterised in three phases. In particular, Phase II is defined by a long period of protein sparing and preferential mobilisation of fat stores, where proteins account for 4% of energy reserves, and lipids for the remaining 96%, followed by Phase III which is defined as a period of increased net protein catabolism (Groscolas and Robin 2001; Robin et al. 1988). Such fasting periods are carefully controlled by metabolic and endocrine shifts, which we can assume are hardcoded early in life. We observed a striking signal of processes directly relatable to this biology in the analysis of unique liver TEGs, including GO terms enriched for digestion, regulation of feeding behaviour, adult feeding behaviour, intestinal absorption, regulation of peptide hormone, cellular response to parathyroid hormone stimulus, and response to starvation. In addition, we identified ten processes related to lipid metabolism, 12 terms related to glucose, glycolysis, and gluconeogenesis, and four terms related to energy production, including NADH regeneration and response to cAMP. We also uncovered the enrichment of fatty acid beta-oxidation in the pink co-expression module (also correlated with the muscle). Moreover, in fasting birds, including Emperor penguins, plasma uric acid concentrations are known to reflect the intensity of protein catabolism (Cherel, Robin, and Maho 1988; Robin et al. 1998), and we identified four enriched GO terms related to uric acid production, including the urea cycle and urate biosynthesis. Additionally, processes involved in vitamin D and retinoic acid were identified, which are also involved in mTOR signalling (Lisse and Hewison 2011). Overall, this strong signature of highly relevant biological processes enriched uniquely in the liver of the Emperor penguin highlights the important role of this organ for crucial metabolic processes later in life. Importantly, we detected these pathways as already under potentially adaptive gene regulation in chicks.

Interestingly, in the analysis of the liver-specific co-expression module (green) the only enriched GO term was blood

coagulation, and this was also identified in analysis of unique tissue-enhanced genes, providing additional evidence on the importance of this process. The liver plays a key role in blood coagulation, where it synthesises the majority of proteins involved in fibrinolysis (Mitra and Metcalf 2012). In humans, metabolic fatty liver disease is associated with heightened blood coagulation (Tripodi et al. 2022; Valenti et al. 2022), indicating that the high rate of lipid storage and subsequent metabolism likely affects the balance of blood coagulation in Emperor penguins. Previous studies have identified blood coagulation genes as under positive selection in *Aptenodytes* penguins (Cole et al. 2022), but the finding of enrichment of this process in Emperor penguins in this study, suggests disparate mechanisms of gene regulation between the two species. Finally, in the analysis of the differentially expressed genes, we identified a suite of candidates involved in lipid metabolism, and glucose and insulin regulation. RET was the most highly over-expressed gene in the liver, which is a well-known oncogene in humans (Jhiang 2000), and has been linked to lipid accumulation in humans (Söhle et al. 2012) and chickens (Sun et al. 2013). LARP1 (also identified in the brain) was also highly over-expressed, and as mentioned above, is a key regulator of the mTOR pathway. Related to this, we also identified MLST8 (also highly over-expressed in the kidney), which is another key regulator of the mTOR pathway, essential for growth signalling and metabolism (Wullschleger, Loewith, and Hall 2006). STC2, a secreted glycoprotein hormone, was identified in all three analyses (Table S16) and is involved in the regulation of insulin-like growth factors and glucose (Sarapio et al. 2019). In obese mice, overexpression of STC2 significantly attenuated fatty liver and hypertriglyceridemia (Zhao et al. 2018). Overall, the strong signature of the liver, together with the muscle, suggests that these two tissues are the most important in gene expression adaptation in Emperor penguins.

## 4.6 | The Role of Gene Expression in Adaptive Evolution

Whilst many of our results are specific to our particular study organisms, their implications for the role of gene expression in accelerating ecological divergence are more general. When a population colonises a new environment, gene expression becomes crucial in guaranteeing persistence, mediating phenotypic plasticity, and contributing to the early phases of adaptive divergence (Pavey et al. 2010; Schlichting and Smith 2002). Genetic divergence in adaptive traits over time eventually leads to reproductive isolation between populations (Fay and Wittkopp 2008; Pavey et al. 2010). Unlike changes in coding genes, changes in gene expression can speed up adaptive evolution, especially under strong selection pressures over short timescales (Brauer, Unmack, and Beheregaray 2017; Jones et al. 2012; Uusi-Heikkilä et al. 2017). More explicitly, changes in gene expression are expected to respond quicker than mutations in coding sequences when a population colonises a novel environment (Josephs 2021). In fact, only gene expression effects can be moulded to enhance or inhibit the role of genes in a particular tissue(s). On the other hand, given recent evidence that a large fraction of genes take part in multiple, if not all, cellular pathways (i.e., the omnigenic model; Boyle, Li, and Pritchard 2017), novel protein mutations which can confer an advantage under specific environmental
conditions whilst not disrupting any of the other protein func-
tions are expected to be rare (i.e., slowly appearing). Under this
perspective, regulatory changes in gene expression could match,
or likely surpass, standing coding sequence diversity as the first
adaptive tool to quickly respond to fast environmental change
(Burton, Ratikainen, and Einum 2022).

Tissue specificity is thought to be hardwired into gene regula-
tory networks, which activate cohorts of genes in particular tis-
sues at particular times (Davidson 2010). However, the evolution
of gene regulatory networks sometimes occurs by mechanisms
that sacrifice specificity, e.g., via network co-option, offering an
efficient mechanism for generating novel phenotypes (McQueen
and Rebeiz 2020). With growing evidence that expression
regulation is heritable (Schadt et al. 2003; Whitehead and
Crawford 2006b), initially plastic gene expression networks may
become hard-wired as ecological divergence proceeds (i.e., ge-
netic assimilation; López-Maury, Marguerat, and Bähler 2008).
Assessing gene networks which are not conserved between spe-
cies offers a valuable approach to measure the process of gene
expression hard-wiring. For example, aside from the plethora of
over-expressed genes related to cold adaptation detected in the
different tissues of the Emperor penguin, we stress the evolu-
tionary relevance of gene network rewiring in the divergence of
the *Aptenodytes* species. Conserved gene co-expression modules
are representative of core physiological functions at an ancestral
level and are evolutionarily maintained through different strong
selective regimes (Barua and Mikheyev 2021; Miller, Horvath,
and Geschwind 2010; Stuart et al. 2003). Alternatively, a lack
of conservation in co-expression networks reflects a rewiring of
gene interactions and can be related to adaptations to different
environmental contexts (Filteau et al. 2013; Oldham, Horvath,
and Geschwind 2006).

Likewise, further investigation of genes from a network per-
spective can give important hints on how regulatory elements
respond to novel selective regimes affecting complex polygenic
traits (Fagny and Austerlitz 2021). For example, the weak pres-
ervation of the Emperor penguin's muscle, liver and brain mod-
ules in the sub-Antarctic sister species suggests gene expression
rewiring related to adaptation to extreme Antarctic conditions.
We can posit that rewiring of the muscle and liver could be
tightly linked to morphological adaptations in the Emperor pen-
guin, such as larger body size and adipose mass (Le Maho 1977).
As previously stressed, rewiring of brain expression is more dif-
ficult to interpret, particularly when using a wild experimental
design, and especially without clear information about brain
morphology in the studied species. However, our study gives a
first hint that a key factor for adaptation to extreme cold may also
lie in the rewiring of tissue expression, like the brain in this case.
Brain rewiring leading to significant morphological differences
over short evolutionary timespans has been previously detected
in the expansion of the human cortex when compared to the
chimpanzee cortex (Oldham, Horvath, and Geschwind 2006).
Thus, even though adaptive differences between the Emperor
penguin and the King penguin are not as clearly connected with
brain evolution (as in the human-chimpanzee comparison),
the Emperor penguin's brain rewiring could be related to be-
havioural adaptations allowing the species to thrive under the
extreme conditions of the Antarctic winter, such as huddling

behaviour and breeding during complete absence of sunlight
(Borboroglu 2015).

5 | Conclusions

Through complementary analyses of a multi-tissue dataset,
we were able to determine which tissues show the strongest
signal of gene expression shifts in the Emperor penguin. We
identified a prominent signature of the muscle and liver and
an apparent rewiring in the brain expression network. The
lack of a gene expression signature in the skin suggests we
should explore the role of coding sequence changes in this
tissue. However, quantifying gene expression differences and
relating these to divergent ecological adaptations is a complex
task and involves the interrogation of only a single aspect of the
extremely complex genotype-to-phenotype association. This
study therefore represents a first contribution to understanding
the role of gene expression in adaptation to the most extreme
environmental conditions experienced on Earth for an endo-
thermic vertebrate. Future research on the functions and roles
of genes showing signals of divergence, either in expression or
in coding sequences, is needed to better understand the initial
phases of ecological divergence and rapid adaptation to novel
selective constraints.

Author Contributions

J.R.P. performed research, analysed data, wrote the paper. F.A.N.F. per-
formed research, analysed data, wrote the paper. F.P. designed research,
performed research, analysed data. S.G. performed research, analysed
data. M.G. performed research, analysed data. A.P. performed research.
M.B. performed research. C.C. performed research. L.Z. performed
research. B.H. analysed data. C.L.B. designed research, performed re-
search, wrote the paper. E.T. designed research, performed research,
wrote the paper. All authors commented and contributed to the final
version of the manuscript.

Acknowledgements

This study was supported by the Institut Polaire Français Paul-Emile
Victor (IPEV) within the framework of the Project 137-ANTAVIA
(PI: CLB), by the Centre Scientifique de Monaco with additional
support from the LIA-647 and RTPI-NUTRESS (CSM/CNRS-
UNISTRA), by the Centre National de la Recherche Scientifique
(CNRS) through the long-term Studies in Ecology and Evolution
(SEE-Life) and Zone Atelier Antarctique et Subantarctique (ZATA)
programs. The study was approved by the French ethics commit-
tee (last: APAFIS#29338-2020070210516365) and the French Polar
Environmental Committee, and permits to handle animals and
access breeding sites were delivered by the "Terres Australes et
Antarctiques Françaises" (TAAF). JRP is supported by funding from
the European Union's Horizon Europe research and innovation pro-
gramme under the Marie Skłodowska-Curie grant agreement No.
101068395-"Poly2Adapt." SG was supported by PNRA16_00099 and
PNRA16_00234. ET was supported by PNRA_16 00164 ("Programma
Nazionale di Ricerca in Antartide." Bando PNRA 5 aprile 2016, n. 651.
- Linea B "Genomica degli adattamenti estremi alla vita in Antartide").
Bioinformatic analyses were performed on the HPC clusters at the
Department of Life and Environmental Sciences ("HappyComputing@
DiSVA"), Marche Polytechnic University, and the Department of Life
Sciences and Biotechnology, University of Ferrara. Open access pub-
lishing facilitated by Universita Politecnica delle Marche, as part of the
Wiley - CRUI-CARE agreement.
Conflicts of Interest

The authors declare no conflicts of interest.

Data Availability Statement

All data pertaining to this study are accessible at the European Nucleotide Archive (ENA) under the study accession PRJEB64484. The Emperor penguin transcriptome can be found under ERS16107788 (HCAY01000001-HCAY01106060). The King penguin transcriptome can be found under the ERS16107787 (HCAX01000001-HCAX01080605). Clean reads from the 3'-QuantSeq mRNA data for each tissue and each species can also be found on the ENA (Study: PRJEB64484) under the sample accessions: ERS16201405—ERS16201453 (Emperor penguin), and ERS16201454—ERS16201503 (King penguin). Additional Data in- cluding the transcriptome annotations, catalogues of long non-coding RNAs, tissue-enhanced transcripts, and Differentially Expressed Genes (DEGs) per tissue can be found on GitHub at https://github.com/ josieparis/EP-KP-transcriptomics, and archived with Zenodo: https:// doi.org/10.5281/zenodo.10218729.

References

Ancel, A., M. Beaulieu, and C. Gilbert. 2013. "The Different Breeding Strategies of Penguins: A Review." *Comptes Rendus Biologies* 336, no. 1:1-12.

Ancel, A., G. L. Kooyman, P. J. Ponganis, et al. 1992. "Foraging Behaviour of Emperor Penguins as a Resource Detector in Winter and Summer." *Nature* 360, no. 6402: 336-339.

Ancel, A., H. Visser, Y. Handrich, D. Masman, and Y. Le Maho. 1997. "Energy Saving in Huddling Penguins." *Nature* 385, no. 6614: 304-305.

Andrews, S. 2010. "FastQC: A Quality Control Tool for High Throughput Sequence Data." https://www.bioinformatics.babraham.ac.uk/proje cts/fastqc/.

Bae, J. Y., J. Woo, S. Kang, and K. O. Shin. 2018. "Effects of Detraining and Retraining on Muscle Energy-Sensing Network and Meteorin-Like Levels in Obese Mice." *Lipids in Health and Disease* 17, no. 1:97.

Barua, A., and A. S. Mikheyev. 2021. "An Ancient, Conserved Gene Regulatory Network Led to the Rise of Oral Venom Systems." *Proceedings of the National Academy of Sciences of the United States of America* 118, no. 14: e2021311118. https://doi.org/10.1073/pnas.20213 11118.

Bell, E. M., and T. V. Callaghan. 2012. "What Are Extreme Environments and What Lives in Them?" In *Life at Extremes: Environments, Organisms and Strategies for Survival*, 1-12. Oxfordshire, UK: CABI.

Bénistant, C., C. Duchamp, F. Cohen-Adad, J. L. Rouanet, and H. Barré. 1998. "Increased In Vitro Fatty Acid Supply and Cellular Transport Capacities in Cold-Acclimated Ducklings (Cairina Moschata)." *American Journal of Physiology* 275, no. 3: R683-R690.

Bentz, A. B., G. W. C. Thomas, D. B. Rusch, and K. A. Rosvall. 2019. "Tissue-Specific Expression Profiles and Positive Selection Analysis in the Tree Swallow (Tachycineta Bicolor) Using a de Novo Transcriptome Assembly." *Scientific Reports* 9, no. 1:1-12.

Bickel, P. E., J. T. Tansey, and M. A. Welte. 2009. "PAT Proteins, an Ancient Family of Lipid Droplet Proteins That Regulate Cellular Lipid Stores." *Biochimica et Biophysica Acta* 1791, no. 6: 419-440.

Bicudo, J. E., C. R. Vianna, and J. G. Chaui-Berlinck. 2001. "Thermogenesis in Birds." *Bioscience Reports* 21, no. 2: 181-188.

Bird, J. P., R. Martin, H. R. Akçakaya, et al. 2020. "Generation Lengths of the World's Birds and Their Implications for Extinction Risk." *Conservation Biology: The Journal of the Society for Conservation Biology* 34, no. 5: 1252-1261.

Bley, N., M. Lederer, B. Pfalz, et al. 2015. "Stress Granules Are Dispensable for mRNA Stabilization During Cellular Stress." *Nucleic Acids Research* 43, no. 4: e26.

Blighe, K., S. Rana, and M. Lewis. 2018. "EnhancedVolcano: Publication-ready volcano plots with enhanced colouring and labeling." https://github.com/kevinblighe/EnhancedVolcano.

Borboroglu, P. G., and P. Dee Boersma. 2015. "Penguins: Natural History and Conservation." University of Washington Press.

Botero, C. A., F. J. Weissing, J. Wright, and D. R. Rubenstein. 2015. "Evolutionary Tipping Points in the Capacity to Adapt to Environmental Change." *Proceedings of the National Academy of Sciences of the United States of America* 112, no. 1: 184-189.

Boyle, E. A., Y. I. Li, and J. K. Pritchard. 2017. "An Expanded View of Complex Traits: From Polygenic to Omnipenic." *Cell* 169, no. 7: 1177-1186.

Brauer, C. J., P. J. Unmack, and L. B. Beheregaray. 2017. "Comparative Ecological Transcriptomics and the Contribution of Gene Expression to the Evolutionary Potential of a Threatened Fish." *Molecular Ecology* 26, no. 24: 6841-6856.

Braun, E. J. 2015. "Chapter 12 Osmoregulatory Systems of Birds." In *Sturkie's Avian Physiology*, edited by C. G. Scanes, 6th ed., 285-300. Amsterdam, The Netherlands: Academic Press.

Buchfink, B., K. Reuter, and H.-G. Drost. 2021. "Sensitive Protein Alignments at Tree-of-Life Scale Using DIAMOND." *Nature Methods* 18, no. 4: 366-368.

Burton, T., I. I. Ratikainen, and S. Einum. 2022. "Environmental Change and the Rate of Phenotypic Plasticity." *Global Change Biology* 28, no. 18: 5337-5345.

Bushmanova, E., D. Antipov, A. Lapidus, and A. D. Prjibelski. 2019. "rnaSPAdes: A de Novo Transcriptome Assembler and Its Application to RNA-Seq Data." *GigaScience* 8, no. 9: giz100.

Bushnell, B. 2014. "BBMap: A Fast, Accurate, Splice-Aware Aligner." https://sourceforge.net/projects/bbmap/.

Campbell-Staton, S. C., Z. A. Cheviron, N. Rochette, J. Catchen, J. B. Losos, and S. V. Edwards. 2017. "Winter Storms Drive Rapid Phenotypic, Regulatory, and Genomic Shifts in the Green Anole Lizard." *Science* 357, no. 6350: 495-498.

Cantalapiedra, C. P., A. Hernández-Plaza, I. Letunic, P. Bork, and J. Huerta-Cepas. 2021. "eggNOG-Mapper v2: Functional Annotation, Orthology Assignments, and Domain Prediction at the Metagenomic Scale." *Molecular Biology and Evolution* 38, no. 12: 5825-5829.

Chapin, F. S., K. Autumn, and F. Pugnaire. 1993. "Evolution of Suites of Traits in Response to Environmental Stress." *American Naturalist* 142: S78-S92.

Chen, J., X. Ren, L. Li, et al. 2019. "Integrative Analyses of mRNA Expression Profile Reveal the Involvement of IGF2BP1 in Chicken Adipogenesis." *International Journal of Molecular Sciences* 20, no. 12: 2923. https://doi.org/10.3390/ijms20122923.

Chen, S., Y. Zhou, Y. Chen, and J. Gu. 2018. "Fastp: An Ultra-Fast All-In-One FASTQ Preprocessor." *Bioinformatics* 34, no. 17: i884-i890.

Cherel, Y., J.-P. Robin, and Y. L. Maho. 1988. "Physiology and Biochemistry of Long-Term Fasting in Birds." *Canadian Journal of Zoology* 66, no. 1: 159-166.

Chesson, P., and N. Huntly. 1997. "The Roles of Harsh and Fluctuating Conditions in the Dynamics of Ecological Communities." *American Naturalist* 150, no. 5: 519-553.

Chevin, L.-M., and A. A. Hoffmann. 2017. "Evolution of Phenotypic Plasticity in Extreme Environments." *Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences* 372, no. 1723: 20160138. https://doi.org/10.1098/rstb.2016.0138.
Chevin, L.-M., R. Lande, and G. M. Mace. 2010. "Adaptation, Plasticity, and Extinction in a Changing Environment: Towards a Predictive Theory." *PLoS Biology* 8, no. 4: e1000357.

Cho, H., Y. D. Yang, J. Lee, et al. 2012. "The Calcium-Activated Chloride Channel Anoctamin 1 Acts as a Heat Sensor in Nociceptive Neurons." *Nature Neuroscience* 15, no. 7: 1015-1021.

Choe, S. S., J. Y. Huh, I. J. Hwang, J. I. Kim, and J. B. Kim. 2016. "Adipose Tissue Remodeling: Its Role in Energy Metabolism and Metabolic Disorders." *Frontiers in Endocrinology* 7: 30.

Cirera, S. 2013. "Highly Efficient Method for Isolation of Total RNA From Adipose Tissue." *BMC Research Notes* 6: 472.

Cole, T. L., C. Zhou, M. Fang, et al. 2022. "Genomic Insights Into the Secondary Aquatic Transition of Penguins." *Nature Communications* 13, no. 1: 3912.

Corbel, H., F. Morlon, S. Geiger, and R. Groscolas. 2009. "State-Dependent Decisions During the Fledging Process of King Penguin Chicks." *Animal Behaviour* 78, no. 4: 829-838.

Corley, S. M., N. M. Troy, A. Bosco, and M. R. Wilkins. 2019. "QuantSeq. 3' Sequencing Combined With Salmon Provides a Fast, Reliable Approach for High Throughput RNA Expression Analysis." *Scientific Reports* 9, no. 1: 1-15.

Corsolini, S., M. Nigro, S. Olmastroni, S. Focardi, and F. Regoli. 2001. "Susceptibility to Oxidative Stress in Adélie and Emperor Penguin." *Polar Biology* 24, no. 5: 365-368.

Coulson-Thomas, V. J., T. F. Gesteira, J. Esko, and W. Kao. 2014. "Heparan Sulfate Regulates Hair Follicle and Sebaceous Gland Morphogenesis and Homeostasis." *Journal of Biological Chemistry* 289, no. 36: 25211-25226.

Davidson, E. H. 2010. "Emerging Properties of Animal Gene Regulatory Networks." *Nature* 468, no. 7326: 911-920.

Davidson, E. H., and D. H. Erwin. 2006. "Gene Regulatory Networks and the Evolution of Animal Body Plans." *Science* 311, no. 5762: 796-800.

Duchamp, C., and H. Barré. 1993. "Skeletal Muscle as the Major Site of Nonshivering Thermogenesis in Cold-Acclimated Ducklings." *American Journal of Physiology* 265, no. 5 Pt 2: R1076-R1083.

Duchamp, C., H. Barré, J. L. Rouanet, et al. 1991. "Nonshivering Thermogenesis in King Penguin Chicks. I. Role of Skeletal Muscle." *American Journal of Physiology* 261: R1438-R1445.

Duchamp, C., J. L. Rouanet, and H. Barré. 2002. "Ontogeny of Thermoregulatory Mechanisms in King Penguin Chicks (Aptenodytes Patagonicus)." *Comparative Biochemistry and Physiology. Part A, Molecular & Integrative Physiology* 131, no. 4: 765-773.

Fagny, M., and F. Austerlitz. 2021. "Polygenic Adaptation: Integrating Population Genetics and Gene Regulatory Networks." *Trends in Genetics: TIG* 37, no. 7: 631-638.

Fay, J. C., and P. J. Wittkopp. 2008. "Evaluating the Role of Natural Selection in the Evolution of Gene Regulation." *Heredity* 100, no. 2: 191-199.

Feld, S., and R. Hirschberg. 1996. "Growth Hormone, the Insulin-Like Growth Factor System, and the Kidney." *Endocrine Reviews* 17, no. 5: 423-480.

Fenzl, A., O. C. Kulterer, K. Spirk, et al. 2020. "Intact Vitamin A Transport Is Critical for Cold-Mediated Adipose Tissue Browning and Thermogenesis." *Molecular Metabolism* 42: 101088.

Filteau, M., S. A. Pavey, J. St-Cyr, and L. Bernatchez. 2013. "Gene Coexpression Networks Reveal Key Drivers of Phenotypic Divergence in Lake Whitefish." *Molecular Biology and Evolution* 30, no. 6: 1384-1396.

Finn, R. D., A. Bateman, J. Clements, et al. 2014. "Pfam: The Protein Families Database." *Nucleic Acids Research* 42: D222-D230.

Finn, R. D., J. Clements, and S. R. Eddy. 2011. "HMMER Web Server: Interactive Sequence Similarity Searching." *Nucleic Acids Research* 39: W29-W37.

Fishman, B., and E. Tauber. 2023. "Epigenetics and Seasonal Timing in Animals: A Concise Review." *Journal of Comparative Physiology A, Neuroethology, Sensory, Neural, and Behavioral Physiology* 210: 565-574. https://doi.org/10.1007/s00359-023-01673-3.

Freddolino, P. L., J. Yang, A. Momen-Roknabadi, and S. Tavazoie. 2018. "Stochastic Tuning of Gene Expression Enables Cellular Adaptation in the Absence of Pre-Existing Regulatory Circuitry." *eLife* 7: e31867. https://doi.org/10.7554/eLife.31867.

Fuentes, P., J. Pelletier, C. Martinez-Herráez, et al. 2021. "The 40S-LARP1 Complex Reprograms the Cellular Translatome Upon mTOR Inhibition to Preserve the Protein Synthetic Capacity." *Science Advances* 7, no. 48: eabg9275.

Gavryushkina, A., T. A. Heath, D. T. Ksepka, T. Stadler, D. Welch, and A. J. Drummond. 2017. "Bayesian Total-Evidence Dating Reveals the Recent Crown Radiation of Penguins." *Systematic Biology* 66, no. 1: 57-73.

Gilbert, C., G. Robertson, Y. Le Maho, Y. Naito, and A. Ancel. 2006. "Huddling Behavior in Emperor Penguins: Dynamics of Huddling." *Physiology & Behavior* 88, no. 4-5: 479-488.

Giménez-Cassina, A., and N. N. Danial. 2015. "Regulation of Mitochondrial Nutrient and Energy Metabolism by BCL-2 Family Proteins." *Trends in Endocrinology and Metabolism: TEM* 26, no. 4: 165-175.

Grabherr, M. G., B. J. Haas, M. Yassour, et al. 2011. "Full-Length Transcriptome Assembly From RNA-Seq Data Without a Reference Genome." *Nature Biotechnology* 29, no. 7: 644-652.

Graja, A., F. Garcia-Carrizo, A.-M. Jank, et al. 2018. "Loss of Periostin Occurs in Aging Adipose Tissue of Mice and Its Genetic Ablation Impairs Adipose Tissue Lipid Metabolism." *Aging Cell* 17, no. 5: e12810.

Grant, P. R., B. R. Grant, R. B. Huey, M. T. J. Johnson, A. H. Knoll, and J. Schmitt. 2017. "Evolution Caused by Extreme Events." *Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences* 372, no. 1723: 20160146. https://doi.org/10.1098/rstb.2016.0146.

Grifone, R., J. Demignon, J. Giordani, et al. 2007. "Eya1 and Eya2 Proteins Are Required for Hypaxial Somitic Myogenesis in the Mouse Embryo." *Developmental Biology* 302, no. 2: 602-616.

Grifone, R., C. Laclef, F. Spitz, et al. 2004. "Six1 and Eya1 Expression Can Reprogram Adult Muscle From the Slow-Twitch Phenotype Into the Fast-Twitch Phenotype." *Molecular and Cellular Biology* 24, no. 14: 6253-6267.

Groscolas, R. 1990. "Metabolic Adaptations to Fasting in Emperor and King Penguins." In *Penguin Biology*, edited by L. S. Davis and J. T. Darby, 269-296. London, UK: Academic Press Limited.

Groscolas, R., and J. Leloup. 1986. "The Endocrine Control of Reproduction and Molt in Male and Female Emperor (*Aptenodytes forsteri*) and Adelie (*Pygoscelis adeliae*) Penguins. II. Annual Changes in Plasma Levels of Thyroxine and Triiodothyronine." *General and Comparative Endocrinology* 63, no. 2: 264-274.

Groscolas, R., and J. P. Robin. 2001. "Long-Term Fasting and Re-Feeding in Penguins." *Comparative Biochemistry and Physiology. Part A, Molecular & Integrative Physiology* 128, no. 3: 645-655.

GTEx Consortium. 2015. "Human Genomics. The Genotype-Tissue Expression (GTEx) Pilot Analysis: Multitissue Gene Regulation in Humans." *Science* 348, no. 6235: 648-660.

Haas, B. J., A. Papanicolaou, M. Yassour, et al. 2013. "De Novo Transcript Sequence Reconstruction From RNA-Seq Using the Trinity Platform for Reference Generation and Analysis." *Nature Protocols* 8, no. 8: 1494-1512.
Heller, E., and E. Fuchs. 2015. "Tissue Patterning and Cellular Mechanics." *Journal of Cell Biology* 211, no. 2: 219-231.

Hirschberg, R., and J. D. Kopple. 1989. "Evidence That Insulin-Like Growth Factor I Increases Renal Plasma Flow and Glomerular Filtration Rate in Fasted Rats." *Journal of Clinical Investigation* 83, no. 1: 326-330.

Hoffmann, A. A., and M. J. Hercus. 2000. "Environmental Stress as an Evolutionary Force." *Bioscience* 50, no. 3: 217-226.

Hoffmann, A. A., and P. A. Parsons. 1997. *Extreme Environmental Change and Evolution*. Cambridge, UK: Cambridge University Press.

Hoffmann, A. A., J. G. Sørensen, and V. Loeschcke. 2003. "Adaptation of Drosophila to Temperature Extremes: Bringing Together Quantitative and Molecular Approaches." *Journal of Thermal Biology* 28, no. 3: 175-216.

Hohtola, E. 2002. "Facultative and Obligatory Thermogenesis in Young Birds: A Cautionary Note." *Comparative Biochemistry and Physiology. Part A, Molecular & Integrative Physiology* 131, no. 4: 733-739.

Jackson, S. T., J. L. Betancourt, R. K. Booth, and S. T. Gray. 2009. "Ecology and the Ratchet of Events: Climate Variability, Niche Dimensions, and Species Distributions." *Proceedings of the National Academy of Sciences of the United States of America* 2, no. Suppl 2: 19685-19692.

Jacob, F., and J. Monod. 1961. "Genetic Regulatory Mechanisms in the Synthesis of Proteins." *Journal of Molecular Biology* 3: 318-356.

Jhiang, S. M. 2000. "The RET Proto-Oncogene in Human Cancers." *Oncogene* 19: 5590-5597.

Jones, F. C., M. G. Grabherr, Y. F. Chan, et al. 2012. "The Genomic Basis of Adaptive Evolution in Threespine Sticklebacks." *Nature* 484, no. 7392: 55-61.

Josephs, E. B. 2021. "Gene Expression Links Genotype and Phenotype During Rapid Adaptation." *Molecular Ecology* 30, no. 1: 30-32.

Kannan, S., J. Hui, K. Mazooji, L. Pachter, and D. Tse. 2016. "Shannon: An Information-Optimal de Novo RNA-Seq Assembler." *bioRxiv:* 039230. https://doi.org/10.1101/039230.

Kim, D.-H., D. D. Sarbassov, S. M. Ali, et al. 2002. "mTOR Interacts With Raptor to Form a Nutrient-Sensitive Complex That Signals to the Cell Growth Machinery." *Cell* 110, no. 2: 163-175.

Klingenspor, M. 2003. "Cold-Induced Recruitment of Brown Adipose Tissue Thermogenesis." *Experimental Physiology* 88, no. 1: 141-148.

Kooyman, G. L. 1993. "Breeding Habitats of Emperor Penguins in the Western Ross Sea." *Antarctic Science/Blackwell Scientific Publications* 5, no. 2: 143-148.

Kozaki, Y., R. Umetsu, Y. Mizukami, et al. 2015. "Peripheral Gene Expression Profile of Mechanical Hyperalgesia Induced by Repeated Cold Stress in SHRSP5/Dmr Rats." *Journal of Physiological Sciences: JPS* 65, no. 5: 417-425.

Kriventseva, E. V., D. Kuznetsov, F. Tegenfeldt, et al. 2019. "OrthoDB v10: Sampling the Diversity of Animal, Plant, Fungal, Protist, Bacterial and Viral Genomes for Evolutionary and Functional Annotations of Orthologs." *Nucleic Acids Research* 47, no. D1: D807-D811.

Krogh, A., B. Larsson, G. von Heijne, and E. L. Sonnhammer. 2001. "Predicting Transmembrane Protein Topology With a Hidden Markov Model: Application to Complete Genomes." *Journal of Molecular Biology* 305, no. 3: 567-580.

Langfelder, P. 2018. "Filtering and collapsing data. Insights from a Billion Correlations." https://peterlangfelder.com/2018/11/25/filtering-and-collapsing-data/.

Langfelder, P., and S. Horvath. 2008. "WGCNA: An R Package for Weighted Correlation Network Analysis." *BMC Bioinformatics* 9: 559.

Langfelder, P., and S. Horvath. 2012. "Fast R Functions for Robust Correlations and Hierarchical Clustering." *Journal of Statistical Software* 46, no. 11: 1-17. https://doi.org/10.18637/jss.v046.i11.

Langfelder, P., R. Luo, M. C. Oldham, and S. Horvath. 2011. "Is My Network Module Preserved and Reproducible?" *PLoS Computational Biology* 7, no. 1: e1001057.

Lekevičius, E. 2009. "Vacant Niches in Nature, Ecology, and Evolutionary Theory: A Mini-Review." *Ekologija* 55, no. 3-4: 165-174.

Le Maho, Y. 1977. "The Emperor Penguin: A Strategy to Live and Breed in the Cold: Morphology, Physiology, Ecology, and Behavior Distinguish the Polar Emperor Penguin From Other Penguin Species, Particularly From Its Close Relative, the King Penguin." *American Scientist* 65, no. 6: 680-693.

Le Maho, Y., P. Delclitte, and J. Chatonnet. 1976. "Thermoregulation in Fasting Emperor Penguins Under Natural Conditions." *American Journal of Physiology* 231, no. 3: 913-922.

Lee, J. O., W. S. Byun, M. J. Kang, et al. 2020. "The Myokine Meteorin-Like (Metrnl) Improves Glucose Tolerance in Both Skeletal Muscle Cells and Mice by Targeting AMPKα2." *FEBS Journal* 287, no. 10: 2087-2104.

Lewden, A., M. R. Enstipp, B. Picard, T. van Walsum, and Y. Handrich. 2017. "High Peripheral Temperatures in King Penguins While Resting at Sea: Thermoregulation Versus Fat Deposition." *Journal of Experimental Biology* 220, no. Pt 17: 3084-3094.

Li, C., Y. Zhang, J. Li, et al. 2014. "Two Antarctic Penguin Genomes Reveal Insights Into Their Evolutionary History and Molecular Changes Related to the Antarctic Environment." *GigaScience* 3, no. 1: 27.

Li, J., Y. Dong, X. Lü, et al. 2013. "Crystal Structures and Biochemical Studies of Human Lysophosphatidic Acid Phosphatase Type 6." *Protein & Cell* 4, no. 7: 548-561.

Lisse, T. S., and M. Hewison. 2011. "Vitamin D: A New Player in the World of mTOR Signaling." *Cell Cycle* 10, no. 12: 1888-1889.

Liu, S., E. D. Lorenzen, M. Fumagalli, et al. 2014. "Population Genomics Reveal Recent Speciation and Rapid Evolutionary Adaptation in Polar Bears." *Cell* 157, no. 4: 785-794.

Lopez, J., and S. W. G. Tait. 2015. "Mitochondrial Apoptosis: Killing Cancer Using the Enemy Within." *British Journal of Cancer* 112, no. 6: 957-962.

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, no. 8: 583-593.

Love, M. I., W. Huber, and S. Anders. 2014. "Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data With DESeq2." *Genome Biology* 15, no. 12: 550.

MacManes, M. D. 2018. "The Oyster River Protocol: A Multi-Assembler and Kmer Approach for de Novo Transcriptome Assembly." *PeerJ* 6: e5428.

Mattick, J. S., P. P. Amaral, P. Carninci, et al. 2023. "Long Non-coding RNAs: Definitions, Functions, Challenges and Recommendations." *Nature Reviews. Molecular Cell Biology* 24: 1-17.

McQueen, E., and M. Rebeiz. 2020. "On the Specificity of Gene Regulatory Networks: How Does Network Co-Option Affect Subsequent Evolution?" *Current Topics in Developmental Biology* 139: 375-405.

Metzger, D. C. H., and P. M. Schulte. 2018. "Similarities in Temperature-Dependent Gene Expression Plasticity Across Timescales in Threespine Stickleback (*Gasterosteus aculeatus*)." *Molecular Ecology* 27, no. 10: 2381-2396.

Mezentseva, N. V., J. S. Kumaratilake, and S. A. Newman. 2008. "The Brown Adipocyte Differentiation Pathway in Birds: An Evolutionary Road Not Taken." *BMC Biology* 6: 17.
Miché, F., B. Vivien-Roels, P. Pévet, C. Spehner, J. P. Robin, and Y. Le Maho. 1991. "Daily Pattern of Melatonin Secretion in an Antarctic Bird, the Emperor Penguin, Aptenodytes Forsteri: Seasonal Variations, Effect of Constant Illumination and of Administration of Isoproterenol or Propranolol." *General and Comparative Endocrinology* 84, no. 2: 249–263.

Miller, J. A., S. Horvath, and D. H. Geschwind. 2010. "Divergence of Human and Mouse Brain Transcriptome Highlights Alzheimer Disease Pathways." *Proceedings of the National Academy of Sciences of the United States of America* 107, no. 28: 12698–12703.

Mitra, V., and J. Metcalf. 2012. "Metabolic Functions of the Liver." *Anaesthesia & Intensive Care Medicine* 13, no. 2: 54–55.

Mo, D., K. Yu, H. Chen, et al. 2017. "Transcriptome Landscape of Porcine Intramuscular Adipocytes During Differentiation." *Journal of Agricultural and Food Chemistry* 65, no. 30: 6317–6328.

Moll, P., M. Ante, A. Seitz, and T. Reda. 2014. "QuantSeq 3' mRNA sequencing for RNA quantification." *Nature Methods* 11, no. 12: i–iii.

Nakamura, K., and S. F. Morrison. 2011. "Central Efferent Pathways for Cold-Defensive and Febrile Shivering." *Journal of Physiology* 589, no. Pt 14: 3641–3658.

Nawrocki, E. P., and S. R. Eddy. 2013. "Infernal 1.1: 100-Fold Faster RNA Homology Searches." *Bioinformatics* 29, no. 22: 2933–2935.

Newman, S. A. 2011. "Thermogenesis, Muscle Hyperplasia, and the Origin of Birds." *BioEssays: News and Reviews in Molecular, Cellular and Developmental Biology* 33, no. 9: 653–656.

Oldham, M. C., S. Horvath, and D. H. Geschwind. 2006. "Conservation and Evolution of Gene Coexpression Networks in Human and Chimpanzee Brains." *Proceedings of the National Academy of Sciences of the United States of America* 103, no. 47: 17973–17978.

Patro, R., G. Duggal, M. I. Love, R. A. Irizarry, and C. Kingsford. 2017. "Salmon Provides Fast and Bias-Aware Quantification of Transcript Expression." *Nature Methods* 14, no. 4: 417–419.

Pavey, S. A., H. Collin, P. Nosil, and S. M. Rogers. 2010. "The Role of Gene Expression in Ecological Speciation." *Annals of the New York Academy of Sciences* 1206, no. 1: 110–129.

Periasamy, M., J. L. Herrera, and F. C. G. Reis. 2017. "Skeletal Muscle Thermogenesis and Its Role in Whole Body Energy Metabolism." *Diabetes & Metabolism Journal* 41, no. 5: 327–336.

Petrenko, A. B., T. Yamakura, H. Baba, and K. Shimoji. 2003. "The Role of N-Methyl-d-Aspartate (NMDA) Receptors in Pain: A Review." *Anesthesia & Analgesia* 97, no. 4: 1108–1116.

Petropoulos, H., P. J. Gianakopoulos, A. G. Ridgeway, and I. S. Skerjanc. 2004. "Disruption of Meox or Gli Activity Ablates Skeletal Myogenesis in P19 Cells." *Journal of Biological Chemistry* 279, no. 23: 23874–23881.

Pirri, F., L. Ometto, S. Fuselli, et al. 2022. "Selection-Driven Adaptation to the Extreme Antarctic Environment in the Emperor Penguin." *Heredity* 129: 317–326. https://doi.org/10.1038/s41437-022-00564-8

Ponganis, P. J., R. P. Van Dam, D. H. Levenson, T. Knower, K. V. Ponganis, and G. Marshall. 2003. "Regional Heterothermy and Conservation of Core Temperature in Emperor Penguins Diving Under Sea Ice." *Comparative Biochemistry and Physiology. Part A, Molecular & Integrative Physiology* 135, no. 3: 477–487.

Prévost, J. 1961. *L'écologie du Manchot empereur Aptenodytes forsteri* (Gray). Paris, France: Éditions Scientifiques Hermann.

Qi, X., Q. Zhang, Y. He, et al. 2019. "The Transcriptomic Landscape of Yaks Reveals Molecular Pathways for High Altitude Adaptation." *Genome Biology and Evolution* 11, no. 1: 72–85.

Radomski, M. W., and T. Orme. 1971. "Response of Lipoprotein Lipase in Various Tissues to Cold Exposure." *American Journal of Physiology* 220, no. 6: 1852–1856.

Raimbault, S., S. Dridi, F. Denjean, et al. 2001. "An Uncoupling Protein Homologue Putatively Involved in Facultative Muscle Thermogenesis in Birds." *Biochemical Journal* 353, no. Pt 3: 441–444.

Rao, R. R., J. Z. Long, J. P. White, et al. 2014. "Meteorin-Like Is a Hormone That Regulates Immune-Adipose Interactions to Increase Beige Fat Thermogenesis." *Cell* 157, no. 6: 1279–1291.

Rappaport, N., N. Nativ, G. Stelzer, et al. 2013. "MalaCards: An Integrated Compendium for Diseases and Their Annotation." *Database: The Journal of Biological Databases and Curation* 2013: bat018.

Ribot, J., F. Felipe, M. L. Bonet, and A. Palou. 2004. "Retinoic Acid Administration and Vitamin A Status Modulate Retinoid X Receptor Alpha and Retinoic Acid Receptor Alpha Levels in Mouse Brown Adipose Tissue." *Molecular and Cellular Biochemistry* 266, no. 1–2: 25–30.

Robin, J.-P., L. Boucontet, P. Chillet, and R. Groscolas. 1998. "Behavioral Changes in Fasting Emperor Penguins: Evidence for a "Refeeding Signal" Linked to a Metabolic Shift." *American Journal of Physiology. Regulatory, Integrative and Comparative Physiology* 274, no. 3: R746–R753.

Robin, J.-P., M. Frain, C. Sardet, R. Groscolas, and Y. Le Maho. 1988. "Protein and Lipid Utilization During Long-Term Fasting in Emperor Penguins." *American Journal of Physiology. Regulatory, Integrative and Comparative Physiology* 254, no. 1: R61–R68.

Romero, I. G., I. Ruvinsky, and Y. Gilad. 2012. "Comparative Studies of Gene Expression and the Evolution of Gene Regulation." *Nature Reviews. Genetics* 13, no. 7: 505–516.

Rothschild, L. J., and R. L. Mancinelli. 2001. "Life in Extreme Environments." *Nature* 409, no. 6823: 1092–1101.

Rowland, L. A., N. C. Bal, and M. Periasamy. 2015. "The Role of Skeletal-Muscle-Based Thermogenic Mechanisms in Vertebrate Endothermy." *Biological Reviews of the Cambridge Philosophical Society* 90, no. 4: 1279–1297.

Sarapio, E., S. K. De Souza, J. F. A. Model, M. Trapp, and R. S. M. Da Silva. 2019. "Stanniocalcin-1 and -2 Effects on Glucose and Lipid Metabolism in White Adipose Tissue From Fed and Fasted Rats." *Canadian Journal of Physiology and Pharmacology* 97, no. 10: 916–923.

Schadt, E. E., S. A. Monks, T. A. Drake, et al. 2003. "Genetics of Gene Expression Surveyed in Maize, Mouse and Man." *Nature* 422, no. 6929: 297–302.

Schlichting, C. D., and H. Smith. 2002. "Phenotypic Plasticity: Linking Molecular Mechanisms With Evolutionary Outcomes." *Evolutionary Ecology* 16, no. 3: 189–211.

Sebag, S. C., Z. Zhang, Q. Qian, et al. 2021. "ADH5-Mediated NO Bioactivity Maintains Metabolic Homeostasis in Brown Adipose Tissue." *Cell Reports* 37, no. 7: 110003.

Shamsi, F., C.-H. Wang, and Y.-H. Tseng. 2021. "The Evolving View of Thermogenic Adipocyte-Ontogeny, Niche and Function." *Nature Reviews. Endocrinology* 17, no. 12: 726–744.

Shao, Z., X. Hou, Y. Chen, Z. Lin, and X. Zhan. 2023. "RNA Sampling Time on Postmortem Avian Carcasses in the Wild." *Wildlife Biology* 2023, no. 6: e01157. https://doi.org/10.1002/wlb3.01157

Simão, F. A., R. M. Waterhouse, P. Ioannidis, E. V. Kriventseva, and E. M. Zdobnov. 2015. "BUSCO: Assessing Genome Assembly and Annotation Completeness With Single-Copy Orthologs." *Bioinformatics* 31, no. 19: 3210–3212.

Simpson, G. G. 1946. "Fossil Penguins." *Bulletin of the American Museum of Natural History* 87, no. 7: 1–100.

Simpson, G. G. 1976. *Penguins: Past and Present, Here and There*. New Haven, Connecticut, USA: Yale University Press.
Sinitsky, M. Y., V. G. Matveeva, M. A. Asanov, and A. V. Ponasenko. 2018. "Modifications in Routine Protocol of RNA Isolation Can Improve Quality of RNA Purified From Adipocytes." *Analytical Biochemistry* 543: 128-131.

Sjöstedt, E., W. Zhong, L. Fagerberg, et al. 2020. "An Atlas of the Protein-Coding Genes in the Human, Pig, and Mouse Brain." *Science* 367, no. 6482: eaay5947. https://doi.org/10.1126/science.aay5947.

Smith, L. R., G. Meyer, and R. L. Lieber. 2013. "Systems Analysis of Biological Networks in Skeletal Muscle Function." *Wiley Interdisciplinary Reviews. Systems Biology and Medicine* 5, no. 1: 55-71.

Smith-Unna, R., C. Boursnell, R. Patro, J. M. Hibberd, and S. Kelly. 2016. "TransRate: Reference-Free Quality Assessment of de Novo Transcriptome Assemblies." *Genome Research* 26, no. 8: 1134-1144.

Söhle, J., N. Machuy, E. Smailbegovic, et al. 2012. "Identification of New Genes Involved in Human Adipogenesis and Fat Storage." *PLoS One* 7, no. 2: e31193.

Song, L., and L. Florea. 2015. "Rcorrector: Efficient and Accurate Error Correction for Illumina RNA-Seq Reads." *GigaScience* 4: 48.

Sparri, R., E. J. Braun, and K. L. Sweazea. 2010. "Putative Role for Insulin in Glucose Regulation by Avian Kidneys." *FASEB Journal: Official Publication of the Federation of American Societies for Experimental Biology* 24, no. S1: lb614-lb614.

Steinhoff, J. S., A. Lass, and M. Schupp. 2022. "Retinoid Homeostasis and Beyond: How Retinol Binding Protein 4 Contributes to Health and Disease." *Nutrients* 14, no. 6: 1236.

Stomberski, C. T., D. T. Hess, and J. S. Stamler. 2019. "Protein S-Nitrosylation: Determinants of Specificity and Enzymatic Regulation of S-Nitrosothiol-Based Signaling." *Antioxidants & Redox Signaling* 30, no. 10: 1331-1351.

Stonehouse, B. 1969. "Environmental Temperatures of Tertiary Penguins." *Science* 163, no. 3868: 673-675.

Stuart, J. M., E. Segal, D. Koller, and S. K. Kim. 2003. "A Gene-Coexpression Network for Global Discovery of Conserved Genetic Modules." *Science* 302, no. 5643: 249-255.

Sun, Y., G. Zhao, R. Liu, et al. 2013. "The Identification of 14 New Genes for Meat Quality Traits in Chicken Using a Genome-Wide Association Study." *BMC Genomics* 14: 458.

Swanson, D. L., Y. Zhang, and A. G. Jimenez. 2022. "Skeletal Muscle and Metabolic Flexibility in Response to Changing Energy Demands in Wild Birds." *Frontiers in Physiology* 13: 961392.

Talbot, D. A., C. Duchamp, B. Rey, et al. 2004. "Uncoupling Protein and ATP/ADP Carrier Increase Mitochondrial Proton Conductance After Cold Adaptation of King Penguins." *Journal of Physiology* 558, no. Pt 1: 123-135.

Taylor, J. R. E. 1986. "Thermal Insulation of the Down and Feathers of Pygoscelid Penguin Chicks and the Unique Properties of Penguin Feathers." *Auk* 103, no. 1: 160-168.

Teufel, F., J. J. Almagro Armenteros, A. R. Johansen, et al. 2022. "SignalP 6.0 Predicts all Five Types of Signal Peptides Using Protein Language Models." *Nature Biotechnology* 40, no. 7: 1023-1025.

Teulier, L., J.-L. Rouanet, D. Letexier, et al. 2010. "Cold-Acclimation-Induced Non-shivering Thermogenesis in Birds Is Associated With Upregulation of Avian UCP but Not With Innate Uncoupling or Altered ATP Efficiency." *Journal of Experimental Biology* 213, no. 4: 2476-2482.

Teulier, L., J.-L. Rouanet, B. Rey, and D. Roussel. 2014. "Ontogeny of Non-shivering Thermogenesis in Muscovy Ducklings (*Cairina moschata*)." *Comparative Biochemistry and Physiology Part A, Molecular & Integrative Physiology* 175: 82-89.

Thomas, D. B., D. T. Ksepka, and R. E. Fordyce. 2011. "Penguin Heat-Retention Structures Evolved in a Greenhouse Earth." *Biology Letters* 7, no. 3: 461-464.

Tinker, D. A., J. T. Brosnan, and G. R. Herzberg. 1984. "Metabolite Levels in Fed and Fasted Chicken Liver and Kidney While Employing a Technique for Artificial Respiration in Birds." *Comparative Biochemistry and Physiology Part B: Comparative Biochemistry* 77, no. 1: 119-124.

Tirosh, I., A. Weinberger, M. Carmi, and N. Barkai. 2006. "A Genetic Signature of Interspecies Variations in Gene Expression." *Nature Genetics* 38, no. 7: 830-834.

Tobler, M., J. L. Kelley, M. Plath, and R. Riesch. 2018. "Extreme Environments and the Origins of Biodiversity: Adaptation and Speciation in Sulphide Spring Fishes." *Molecular Ecology* 27, no. 4: 843-859.

Tripodi, A., R. Lombardi, M. Primignani, V. La Mura, F. Peyvandi, and A. L. Fracanzani. 2022. "Hypercoagulability in Patients With Non-Alcoholic Fatty Liver Disease (NAFLD): Causes and Consequences." *Biomedicine* 10, no. 2: 249.

Uusi-Heikkilä, S., T. Sävilammi, E. Leder, R. Arlinghaus, and C. R. Primmer. 2017. "Rapid, Broad-Scale Gene Expression Evolution in Experimentally Harvested Fish Populations." *Molecular Ecology* 26, no. 15: 3954-3967.

Valenti, L., A. Tripodi, V. La Mura, et al. 2022. "Clinical and Genetic Determinants of the Fatty Liver-Coagulation Balance Interplay in Individuals With Metabolic Dysfunction." *JHEP Reports: Innovation in Hepatology* 4, no. 12: 100598.

Vianna, J. A., F. A. N. Fernandes, M. J. Frugone, et al. 2020. "Genome-Wide Analyses Reveal Drivers of Penguin Diversification." *Proceedings of the National Academy of Sciences of the United States of America* 117, no. 36: 22303-22310.

Vriens, J., B. Nilius, and T. Voets. 2014. "Peripheral Thermosensation in Mammals." *Nature Reviews. Neuroscience* 15, no. 9: 573-589.

Wakatsuki, K., Y. T-Uchimura, T. Matsubara, T. Nasu, K. Mizumura, and T. Taguchi. 2021. "Peripheral Nociceptive Mechanisms in an Experimental Rat Model of Fibromyalgia Induced by Repeated Cold Stress." *Neuroscience Research* 162: 22-30.

Wang, Z., Y. Kontani, Y. Sato, T. Mizuno, N. Mori, and H. Yamashita. 2003. "Muscle Type Difference in the Regulation of UCP3 Under Cold Conditions." *Biochemical and Biophysical Research Communications* 305, no. 2: 244-249.

Wani, A. K., N. Akhtar, F. Sher, A. A. Navarrete, and J. H. P. Américo-Pinheiro. 2022. "Microbial Adaptation to Different Environmental Conditions: Molecular Perspective of Evolved Genetic and Cellular Systems." *Archives of Microbiology* 204, no. 2: 144.

Whitehead, A., and D. L. Crawford. 2006a. "Neutral and Adaptive Variation in Gene Expression." *Proceedings of the National Academy of Sciences of the United States of America* 103, no. 14: 5425-5430.

Whitehead, A., and D. L. Crawford. 2006b. "Variation Within and Among Species in Gene Expression: Raw Material for Evolution." *Molecular Ecology* 15, no. 5: 1197-1211.

Williams, C. L., J. C. Hagelin, and G. L. Kooyman. 2015. "Hidden Keys to Survival: The Type, Density, Pattern and Functional Role of Emperor Penguin Body Feathers." *Proceedings of the Royal Society B: Biological Sciences* 282, no. 1817: 20152033.

Wingfield, J. C., J. P. Kelley, and F. Angelier. 2011. "What Are Extreme Environmental Conditions and How Do Organisms Cope With Them?" *Current Zoology* 57, no. 3: 363-374.

Wittkopp, P. J. 2007. "Variable Gene Expression in Eukaryotes: A Network Perspective." *Journal of Experimental Biology* 210, no. Pt 9: 1567-1575.

Wu, B., L. Zhang, Y.-H. Zhu, et al. 2018. "Mesoderm/Mesenchyme Homeobox Gene 1 Promotes Vascular Smooth Muscle Cell Phenotypic Modulation and Vascular Remodeling." *International Journal of Cardiology* 251: 82-89.
Wu, C.-W., and K. B. Storey. 2021. "mTOR Signaling in Metabolic Stress Adaptation." *Biomolecules* 11, no. 5: 681.

Wullschleger, S., R. Loewith, and M. N. Hall. 2006. "TOR Signaling in Growth and Metabolism." *Cell* 124, no. 3: 471-484.

Xiao, L., Z. Zhao, F. He, and Z. Du. 2019. "Multivariable Regulation of Gene Expression Plasticity in Metazoans." *Open Biology* 9, no. 12: 190150.

Yang, J., W.-R. Li, F.-H. Lv, et al. 2016. "Whole-Genome Sequencing of Native Sheep Provides Insights Into Rapid Adaptations to Extreme Environments." *Molecular Biology and Evolution* 33, no. 10: 2576-2592.

Yu, G., L.-G. Wang, Y. Han, and Q.-Y. He. 2012. "clusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters." *OMICS: A Journal of Integrative Biology* 16, no. 5: 284-287.

Yu, M., Z. Yue, P. Wu, et al. 2004. "The Biology of Feather Follicles." *International Journal of Developmental Biology* 48, no. 2-3: 181-191.

Zachos, J., M. Pagani, L. Sloan, E. Thomas, and K. Billups. 2001. "Trends, Rhythms, and Aberrations in Global Climate 65 ma to Present." *Science* 292, no. 5517: 686-693.

Zaefarian, F., M. R. Abdollahi, A. Cowieson, and V. Ravindran. 2019. "Avian Liver: The Forgotten Organ." *Animals: An Open Access Journal From MDPI* 9, no. 2: 63.

Zhang, B., and S. Horvath. 2005. "A General Framework for Weighted Gene Co-Expression Network Analysis." *Statistical Applications in Genetics and Molecular Biology* 4: Article17.

Zhao, J., Y. Jiao, Y. Song, et al. 2018. "Stanniocalcin 2 Ameliorates Hepatosteatosis Through Activation of STAT3 Signaling." *Frontiers in Physiology* 9: 873.

Zhuo, Z., S. J. Lamont, W. R. Lee, and B. Abasht. 2015. "RNA-Seq Analysis of Abdominal Fat Reveals Differences Between Modern Commercial Broiler Chickens With High and Low Feed Efficiencies." *PLoS One* 10, no. 8: e0135810.

## Supporting Information

Additional supporting information can be found online in the Supporting Information section.