We generated Illumina sequence data for two woolly mammoths that died ∼20,000 and ∼60,000 years ago (), including individuals from the two major lineages of woolly mammoths, clade I (individual M4) and clade II (M25), which are estimated to have diverged ∼1.5 Ma (), and three extant Asian elephants (Elephas maximus). We aligned sequencing reads to the genome assembly for the African Savannah elephant (Loxodonta africana), resulting in non-redundant average sequence coverage of ∼20-fold for each mammoth and ∼30-fold for each Asian elephant ( Figure S1 ). We identified ∼33 million putative single-nucleotide variants (SNVs) among the three elephantid species (see Experimental Procedures for details), including ∼1.4 million nucleotide variants fixed for the derived allele in the two mammoths, but for the ancestral allele in the African and Asian elephants. Among the variants were 2,020 fixed, mammoth-derived amino acid substitutions in 1,642 protein-coding genes and 26 protein-coding genes with premature stop codons (putative LOF substitutions).

We also inferred the functional significance of fixed, derived LOF substitutions in woolly mammoth genes. We identified a single KEGG term enriched among the genes with LOF substitutions, fat digestion and absorption (E = 127.64, hypergeometric p = 1.0 × 10, FDR q = 1.0 × 10), and 48 KO terms enriched among these genes at an FDR ≤ 0.10 ( Figure 3 A). Enriched KO terms were almost exclusively related to cholesterol, sterol, triglyceride, and lipid homeostasis and metabolism ( Figure 3 B), such as decreased circulating cholesterol level (E = 33.15, hypergeometric p = 5.7 × 10, FDR q = 4.3 × 10), decreased sterol level (E = 30.15, hypergeometric p = 7.6 × 10, FDR q = 4.3 × 10), and abnormal circulating lipid level (E = 30.15, hypergeometric p = 7.6 × 10, FDR q = 4.3 × 10).

(B) Word cloud of 26 selected mouse KO phenotypes enriched among the protein-coding genes with fixed, derived premature stop codons in woolly mammoths. Phenotype terms are scaled to the log2 enrichment of that phenotype and color coded by −log10 p value of phenotype enrichment (hypergeometric test).

(A) Manhattan plot of mouse KO phenotypes enriched among genes with fixed, derived loss premature stop codons in woolly mammoths. The –log10(hypergeometric p values) are shown for each phenotype; phenotypes are grouped by anatomical system effected. Vertical red line, FDR = 0.1.

We found that genes with fixed, derived woolly mammoth substitutions were enriched for 40 KEGG pathways and 859 mouse KO phenotypes, at a false discovery rate (FDR) ≤ 0.10 ( Figure 2 A). Significantly enriched KEGG pathways included circadian rhythm – mammal (enrichment [E] = 6.71, hypergeometric p = 2.7 × 10, FDR q = 0.02), fat digestion and absorption (E = 4.01, hypergeometric p = 7.9 × 10, FDR q = 0.05), complement and coagulation cascades (E = 4.28, hypergeometric p = 5.0 × 10, FDR q = 6.7 × 10), and metabolic pathways (E = 8.39, hypergeometric p = 2.2 × 10, FDR q = 1.6 × 10) ( Table S1 ). Enriched KO phenotypes included decreased core body temperature (E = 4.15, hypergeometric p = 8.0 × 10, FDR q = 7.2 × 10), abnormal brown adipose tissue morphology (E = 2.99, hypergeometric p = 1.4 × 10, FDR q = 4.0 × 10), abnormal thermal nociception (E = 2.25, hypergeometric p = 5.4 × 10, FDR q = 0.05), abnormal glucose homeostasis (E = 1.46, hypergeometric p = 2.6 × 10, FDR q = 3.2 × 10), and many body mass-/weight-related phenotypes ( Figure 2 B).

(B) Word cloud of 40 selected mouse KO phenotypes enriched among the protein-coding genes with fixed, derived mammoth amino acid changes. Phenotype terms are scaled to the log2 enrichment of that phenotype and color coded by −log10 p value of phenotype enrichment (hypergeometric test).

(A) Manhattan plot of mouse KO phenotypes enriched among genes with fixed, derived mammoth amino acid changes. The –log10(hypergeometric p values) are shown for each phenotype; phenotypes are grouped by anatomical system effected. The 551 genes with fixed, derived mammoth amino acid changes have mouse KO data. Horizontal red line, FDR = 0.1.

We used several complementary approaches to infer the putative functional consequences of mammoth-specific amino acid substitutions, including classifying substitutions based on their BLOSUM80 exchangeabilities (), predicted functional consequences based on PolyPhen-2 (), and the inter-species conservation of sites at which substitutions occurred, as well as identifying Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways () and mouse knockout (KO) phenotypes () enriched among protein-coding genes with fixed, derived amino acid substitutions in the wooly mammoth. Finally, we manually selected gene-pathway and gene-phenotype associations for further study according to the following two criteria: (1) the richness of literature supporting the role of each gene in specific pathways and phenotypes; and (2) the exchangeability, PolyPhen-2 score, and strength of sequence conservation at sites with mammoth-specific substitutions.

Mouse Genome Database Group The Mouse Genome Database: integration of and access to knowledge about the laboratory mouse.

Previous studies have found hypomorphic polymorphisms in the woolly mammoth MC1R associated with reddish fur color (), but these variants may have been relatively rare in mammoth populations (). Thus, variants in other genes also may have contributed to coat color variability in mammoths, which varied from blonde to orange to nearly black (). We identified 38 genes with mammoth-specific amino acid changes associated with abnormal coat/hair morphology in KO mice, including derived substitutions in eight genes specifically associated with diluted coat color. We also found that the expression of genes with fixed, derived woolly mammoth substitutions were enriched in hair root sheath (hypergeometric p = 0.006), coat hair follicle (hypergeometric p = 0.013), hair follicle (hypergeometric p = 0.016), skin (hypergeometric p = 0.018), and hair outer root sheath (hypergeometric p = 0.018).

Woolly mammoths evolved a suite of morphological adaptations to life in the extreme cold, including small ears and tails, a long thick coat, and, unlike other elephants, numerous large sebaceous glands, which are thought to have helped repel water and improve insulation (). Woolly mammoths also evolved a characteristic set of skeletal traits, including a high, domed skull with dorsally expanded parietals; an anterio-posteriorly compressed skull; and a sloping back. Consistent with mammoth-specific amino acid changes contributing to these traits, we found that genes with mammoth-specific substitutions were enriched in KO phenotypes such as abnormal tail morphology (E = 1.71, hypergeometric p = 2.0 × 10, FDR q = 2.2 × 10), abnormal tail bud morphology (E = 5.06, hypergeometric p = 3.0 × 10, FDR q = 3.2 × 10), small tail bud (E = 18.00, hypergeometric p = 4.1 × 10, FDR q = 1.6 × 10), abnormal ear morphology (E = 1.60, hypergeometric p = 9.0 × 10, FDR q = 6.4 × 10), cup-shaped ears (E = 5.06, hypergeometric p = 3.0 × 10, FDR q = 1.4 × 10), and abnormal sebaceous gland morphology (E = 2.33, hypergeometric p = 8.0 × 10, FDR q = 6.3 × 10), including substitutions in three genes leading to enlarged sebaceous glands. We also found numerous enriched KO phenotypes that affected the skeleton, including domed cranium (E = 2.26, hypergeometric p = 1.3 × 10, FDR q = 8.3 × 10), abnormal parietal bone morphology (E = 2.66, hypergeometric p = 2.6 × 10, FDR q = 1.9 × 10), and short snout (E = 2.39, hypergeometric p = 2.0 × 10, FDR q = 2.3 × 10).

Fixed, derived mammoth-specific amino acid substitutions occurred in eight genes associated with circadian biology, including those that play central roles in maintaining normal circadian rhythms and entraining the circadian clock to external stimuli such as temperature. HRH3 and PER2 KO mice, for example, have abnormal circadian temperature homeostasis (). PER2 directly mediates the early adaptive response to shifted temperature cycles and coordinates adaptive thermogenesis by synchronizing UCP1 expression and activation in brown adipose tissue (). Similarly, neuronal histamine receptors regulate circadian energy homeostasis through UCP1 expression in brown adipose tissue. HRH1 KO mice, for example, have abnormal circadian rhythms and abnormal circadian feeding behaviors, including a shift in food consumption from day to night (). These observations suggest that the circadian system in woolly mammoths may have adapted to the extreme seasonal light-dark oscillations of the high arctic.

Organisms living at high latitudes in the arctic experience long periods of darkness during the winter and near constant light in the summer, which prevents polar-adapted species from utilizing daily light-dark cycles to entrain their circadian clocks. Svalbard reindeer (Rangifer tarandus platyrhynchus), for example, have lost functioning circadian clocks and circadian rhythmicity in PER2 and BMAL1 (ARNTL) expression (). Moreover, several other arctic species also are known to have derived circadian systems (), and we observed that several enriched KO and KEGG terms were related to circadian biology, motivating us to explore circadian genes in greater detail.

Animal activity around the clock with no overt circadian rhythms: patterns, mechanisms and adaptive value.

The enrichment of genes with derived amino acid ( Figure 2 ) and LOF substitutions ( Figure 3 ) in woolly mammoths that function in lipid metabolism, adipose development, and physiology suggests modifications of these processes may have played an important role in the evolution of woolly mammoths and adaptation to arctic life. We identified 54 genes with fixed, derived amino acid substitutions and KO phenotypes that affect adipose tissue, including phenotypes that alter both the location and abundance of white and brown fat deposits throughout the body. Among the genes with woolly mammoth-specific substitutions were the leptin receptor (LEPR); DLK1 (also known as preadipocyte factor 1), an epidermal growth factor repeat-containing transmembrane protein that regulates adipocyte differentiation; the growth hormone receptor (GHR); and corticotropin-releasing hormone (CRH). We also identified 39 genes with KO phenotypes that affect insulin signaling, and found that genes with mammoth-specific amino acid substitutions were enriched in several KO phenotypes related to insulin signaling, including abnormal circulating insulin level (E = 1.82, hypergeometric p = 1.0 × 10, FDR q = 1.5 × 10), insulin resistance (E = 2.23, hypergeometric p = 3.0 × 10, FDR q = 3.5 × 10), and impaired glucose tolerance (E = 1.91, hypergeometric p = 4.0 × 10, FDR q = 3.7 × 10).

To infer the putative consequences of woolly mammoth-specific amino acid substitutions in thermoTRPs, we generated structural models of the ancestral Asian elephant/mammoth (AncGajah; Figure 1 A) and ancestral mammoth (AncMammoth; Figure 1 A) TRPA1, which mediates nocifensive () and vascular responses to noxious cold () as well as generally potentiating responses to noxious stimuli (), and TRPV4, which mediates autonomic and behavioral responses to cold (). We found that the elephantid TRPA1 and TRPV4 proteins were predicted to adopt the common TRP channel structure, which is composed of a series of amino terminal ankyrin repeats (ARD), separated by a membrane proximal domain (MPD) from the six transmembrane helices (S1–S6) that form the ion-permeable pore in tetrameric channels ( Figure 4 B). We found that the TRPA1 R1031T substitution occurred in an unstructured loop between the TRP-like domain and the C-terminal coiled-coil domain ( Figures 4 C and 4D), which is predicted to alter the electrostatic surface by reducing the local positive charge ( Figure 4 E). The mammoth-specific TRPV4 V658I substitution occurred at the first site in the S6 helix ( Figure 5 A), which is part of the outer pore region important for activation of the channel in response to heat ( Figure 5 B). Indeed, we found that site 658 is located within a cluster of sites that mediate heat activation in the related TRPV3 channel (), and it is homologous to a site in TRPV3 that adopts temperature-dependent conformations ( Figures 5 C and 5D;). Site 658 also mediates the interaction between TRPV channels and the agonist vanillotoxin DkTx ( Figure 5 E;). These data suggest the mammoth-specific R1031T and V658I substitutions may have affected the gating dynamics in TRPA1 and TRPV4, respectively.

(E) Site 658 also mediates the interaction between TRPV channels and the spider vanillotoxin DkTx (pink). I658 is shown as a magenta sphere in the woolly mammoth TRPV4 model (blue), and the experimentally determined structure of mouse TRPV1 (gray) complexed with DkTx is superimposed onto the mammoth structural model.

(D) Close-up of the region boxed in (B). I658 is shown as a magenta-colored sphere, homologous residues in mouse TRPV3 that mediate temperature sensitivity are shown as red spheres, and residues with temperature-dependent conformations in TRPV3 are shown as yellow spheres.

(C) Conservation of the pore helix-pore loop-S6 region between TRPV3 (top) and TRPV4 (bottom). Residues in TRPV3 that have been experimentally shown to mediate temperature sensing and that have temperature-dependent conformations are shown in red and yellow, respectively. Homologous residues in TRPV4 are shown in dark gray, and site 658 is shown in magenta.

(B) Cartoon representation of the pore domain of the woolly mammoth TRPV4 homology model. The location of the V658I substitution is shown as a magenta-colored sphere.

(A) Diagram of the major structural domains of TRPV4. Gray regions were not included in the TRPV4 structural model. The location of the mammoth-specific V658I substitution is shown in magenta.

TRPA1 contributes to cold, mechanical, and chemical nociception but is not essential for hair-cell transduction.

The most intriguing mouse KO phenotype enriched among genes with woolly mammoth-specific amino acid changes was abnormal thermal nociception (13 genes). For example, we identified woolly mammoth-specific amino acid changes in five temperature-sensitive transient receptor potential (thermoTRP) channels ( Figure 4 A) that sense noxious cold (TRPM8) (), innocuous warmth (TRPV3 and TRPV4) (), or noxious cold or heat depending on species (TRPA1) () or that are heat sensitive but not known to be involved in temperature sensation (TRPM4) (). We also identified a mammoth-specific amino acid change in PIRT, a small phosphoinositide-binding protein that functions as a regulatory subunit of TRPM8 and the noxious heat sensor TRPV1 ().

(E) Close-up of the region boxed in (D) in the TRPA1 homology model of the AncGajah ancestor (left) and woolly mammoth (right). The electrostatic surfaces of the proteins are shown (blue, positive; red, negative). The superimposed square indicates the location of the charge altering R1031T substitution.

(D) Cartoon representation of the pore domain of the TRPA1 homology model. The location of the R1031T substitution is shown as magenta-colored spheres; helix coloring follows (C).

(C) Diagram of the major structural domains of TRPA1. Gray regions were not included in the TRPA1 structural model. The location of the mammoth-specific R1031T substitution is shown.

(B) Structure of a single TRP subunit (left) and the tetrameric channel (right) viewed from the side. The ankyrin-repeat domain (ARD), transmembrane domains (S1–S6), membrane-proximal domain (MPD), C-terminal domain (CTD), TRP box, pore loop, and pore turret are labeled. Amino acids within the ARD, MPD, pore turret, the outer pore region and in the initial part of S6, the TRP box, and the CTD influence temperature sensing in TRPV and TRPA channels.

(A) Temperature ranges over which TRPM8, TRPA1, TRPM4, TRPV4, and TRPV3 are active. Threshold temperatures for each channel are shown.

TRPM8, but not TRPA1, is required for neural and behavioral responses to acute noxious cold temperatures and cold-mimetics in vivo.

TRPM8, but not TRPA1, is required for neural and behavioral responses to acute noxious cold temperatures and cold-mimetics in vivo.

Thermal Tuning of the Woolly Mammoth Temperature Sensor TRPV3

Cao et al., 2013 Cao E.

Liao M.

Cheng Y.

Julius D. TRPV1 structures in distinct conformations reveal activation mechanisms. Liao et al., 2013 Liao M.

Cao E.

Julius D.

Cheng Y. Structure of the TRPV1 ion channel determined by electron cryo-microscopy. We inferred the structural consequences of the N647D substitution by generating homology models of the AncGajah and AncMammoth TRPV3 protein and tetrameric channel (as described above). We found that the TRPV3 models were structurally very similar to the TRPV1 reference structure on which they were based (root-mean-square deviation [RMSD] = 1.93–1.99 Å; Figure S4 ), particularly in the α helices that form the pore of the tetrameric channel, suggesting that these are realistic models. In the TRPV1 structure, hydrogen bonds between residues in the pore loop are thought to maintain the outer pore in a non-conductive conformation in the closed state; conformational changes in the pore helix and pore loops disrupt these local hydrogen bonds to facilitate gating and widening of the selectivity filter upon channel activation (). Our structural model suggests that the carbonyl oxygen of the ancestral N647 residue forms a hydrogen bond with the neighboring side chain of Q645, whereas in the AncMammoth structure these hydrogen bonds are replaced by a pair of hydrogen bonds between D647 and K610, potentially impeding full opening of the channel in mammoths ( Figure 6 C). Consistent with this prediction, in the model of the AncGajah open channel, the distance between diagonally opposed G637 residues is 8.5 Å, whereas this distance is only 6.3 Å in the AncMammoth open channel ( Figure 6 D), which narrows the pore diameter at the selectivity filter by ∼26% and the pore radius at the selectivity filter by ∼60% ( Figure 6 E; Figures S5 A–S5H).