Accumulation and functional architecture of deleterious genetic variants during the extinction of Wrangel Island Mammoths

Woolly mammoths were among the most abundant cold adapted species during the Pleistocene. Their once large populations went extinct in two waves, an end-Pleistocene extinction of continental populations1-4 followed by the mid-Holocene extinction of relict populations on St. Paul Island ∼5,600 years ago5 and Wrangel Island ∼4,000 years ago1-4. Population and conservation genetics theory predicts that deleterious alleles will accumulate as populations decline, leading to downward spiral of declining population size and fitness ending in extinction (mutational meltdown). While the extinction of Wrangel Island mammoths was preceded by prolonged demographic decline, reduced population size and genetic diversity, recurrent inbreeding, and the fixation of deleterious alleles3,4,6-8, the functional consequences of these processes are unclear. Here we show that the extinction of Wrangel Island mammoths was accompanied by an accumulation of deleterious mutations that are predicted to cause diverse behavioral and developmental defects. Resurrecting and functional characterization of Wrangel Island mammoth genes with putative deleterious substitutions identified both loss and gain of function mutations associated with ciliopathies, oligozoospermia and reduced fertility, and neonatal diabetes. These results indicate that last mammoths likely suffered from genetic disease that reduced fitness and directly contributed to their extinction.


Island mammoth genes with putative deleterious substitutions identified both loss and
gain of function mutations associated with ciliopathies, oligozoospermia and reduced fertility, and neonatal diabetes. These results indicate that last mammoths likely suffered from genetic disease that reduced fitness and directly contributed to their extinction.
Woolly mammoths (Mammuthus primigenius) were among the most abundant cold adapted megafaunal species during the Middle to Late Pleistocene (ca 116-12 Kyr BP), inhabiting a large swath of steppe-tundra that extended from Western Europe, through Asia and Beringia, into North America. Paleontological and genetic data indicate mammoths went extinct in at least two waves, an end-Pleistocene decline and extinction of continental populations [1][2][3][4] followed by the mid-Holocene extinction of relict populations on St. Paul Island ~5,600 years ago 5 and Wrangel Island ~4,000 years ago [1][2][3][4] (Fig. 1). Extinction of Wrangel Island mammoths was preceded by prolonged demographic decline, resulting in a small population, reduced genetic diversity, recurrent breeding among distant relatives, and the fixation of deleterious alleles 3,4,6-8 . The functional significance of putatively deleterious amino acid substitutions in the Wrangel Island mammoth, however, is unknown. Thus it is unclear if an accumulation of deleterious variants contributed to the extinction of Wrangel Island mammoths.
To infer if the extinction of Woolly mammoths was associated with an accumulation of putatively deleterious amino acid variants, we identified homozygous nonsynonymous substitutions unique (private) in three extant Asian elephants (Elephas maximus) and three woolly mammoths, the ~44,800 year old Oimyakon mammoth 1 , the ~20,000 year old M4 mammoth [9][10][11][12][13] , and the ~4,300 years old Wrangel Island mammoth 1 . These mammoths span the age from when mammoth populations were large and wide-spread (Oimyakon), to near the beginning of their final decline (M4), and their last known population (Wrangel Island); thus these individuals allow us to infer if deleterious mutations accumulated in the mammoth genome during their decline and eventual extinction. Next we used PolyPhen-2 14,15 to computationally predict the functional impact of each private homozygous amino acid substitution ( Table 1). 213300]), a milder disorder characterized by defects in the cerebellum and brain stem leading to impaired balance and coordination (Oka et al., 2016).
The Wrangel Island mammoth P119L mutation occurs in a highly conserved region of the protein, which is invariant for proline or serine across vertebrates and is therefore potentially deleterious (Fig. 3a). To infer if this mutation had functional consequences, we used a Xenopus model of ciliogenesis 16  Another gene with a 'probably damaging' amino acid substitution in the Wrangel Island mammoth associated with male infertility is Naked cuticle 1 (NKD1), which encodes passive antagonist of the Wnt/TCF-LEF signaling pathway [18][19][20] . The Wrangel Island mammoth A88V substitution occurred at a site that is nearly invariant for alanine or serine in diverse vertebrates ( Fig. 4a), suggesting it may have had functional consequences. To infer if the A88V substitution had functional affects, we resurrected the Wrangel Island/M4 ancestral (AncYakut, Fig. 2a) and Wrangel Island NKD1 genes and tested their ability to antagonize Wnt-signaling in elephant dermal fibroblasts transiently transfected with a luciferase reporter vector containing a minimal promoter and eight copies of a TCF-LEF response element (pGL4.49[luc2P/TCF-LEF/Hygro]) and treated with a small molecule agonist of the Wnt-signaling pathway (CHIR99021). The AncYakut NKD1 reduced luminescence to background levels in response to CHIR99021 treatment (Wilcox test, P=2.71×10 -6 ), in stark contrast however, the Wrangel Island NKD1 did not affect luciferase expression (Wilcox P=0.98), indicating the A88V substitution is a loss of function mutation (Fig. 4b). Transgenic mice with loss of function mutations in NKD1 have dysregulated Wnt/beta-catenin signaling in the testis leading to abnormal seminiferous tubule morphology, small seminiferous tubules, small testis, oligozoospermia, and reduced fertility 21,22 , suggesting the A88V substitution may have affected fertility in Wrangel Island mammoths.
Although Wrangel Island mammoths were generally smaller than mainland mammoth populations, their size varied and large and small individuals coexisted on the Island. Four genes with 'probably damaging' amino acid substitutions in the Wrangel Island genome caused 'decreased body size' (MP:0001265) and five genes with 'probably damaging' amino acid substitutions caused 'decreased body weight' (MP:0001262) in mouse knockouts. In stark contrast, genes with 'probably damaging' amino acid substitutions were not associated with mouse knockout phenotypes such as 'increased body size' (MP:0001264) and only a single gene was associated with 'increased body weight' (MP:0001262). The proportion of genes with 'probably damaging' substitutions associated with 'decreased body size' and 'decreased body weight' phenotypes (7/38 = 0.184) was significantly greater (binomial test N≥7, P=4.48×10 -5 ) than the proportion of genes associated with 'increased body size' and 'increased body weight' phenotypes (1/38 = 0.026), suggesting these variants may have contributed to variation in size among Wrangel Island mammoths.
Among the genes with 'probably damaging' mutations in the Wrangel Island mammoth associated with decreased body size was NEUROGENIN 3 (NEUROG3), which encodes a basic helix-loop-helix transcription factor that is required for endocrine cell development in the pancreas and intestine. The 'probably damaging' G195E substitution in the Wrangel Island mammoth NEUROG3 protein occurred at a site that is nearly invariant for glycine across mammals ( Fig. 5a/b) within an LXXLL motif in the C-terminal transcriptional activation domain 23 . To infer if the G195E substitution had functional affects, we resurrected the AncYakut and Wrangel Island NEUROG3 genes and tested their ability to transactivate luciferase expression from a reporter vector containing a minimal promoter and six repeats of the PAX4 E-box (pGL3[luc/6x-PAX4E/minP]) in transiently transfected elephant dermal fibroblasts (Fig. 5c) Previous studies of Wrangel Island mammoths have found signatures of reduced genetic diversity, reduced population size, and recurrent breeding among distant relatives 1 , and an increased burden of deletions, retrogenes, genes with deleted exons, and genes with premature stop codons 6 . The functional and phenotypic consequences of these findings, however, are unclear. We found that the Wrangel Island mammoth genome had significantly more putatively deleterious mutations than older continental mammoths, and that these mutations were associated abnormal phenotypes in knockout mice. These data suggest that Wrangel Island mammoths suffered from genetic disorders resulting from the accumulation of deleterious mutations in an inbred and rapidly declining population. Indeed, our functional assays suggest that the Wrangel Island mammoth suffered from ciliopathies, reduced male fertility, and pathological insulin signaling. Thus an accumulating burden of deleterious mutations in Wrangel Island mammoths may have contributed to their extinction.

Genome assembly
Details of the sequencing protocol for the Oimyakon and Wrangel Island mammoths can be found in Palkopoulou et al. (2015) and for the Asian elephants, M25, and M4 in Lynch et al. The sequences from the mammoths were treated separately to account for DNA damage in the sequences. Putative adapter sequences were removed and we merged overlapping paired-end reads using available scripts (Kircher, 2012). We required an overlap of at least 11 nucleotides between the mates, and only pairs that could be merged were retained for subsequent analyses. The merged reads were aligned to the genome from the African elephant (loxAfr3) using BWA with default parameters, and only the mapped reads that were longer than 20 bps were retained for the subsequent SNP calls. The reads were realigned using the GATK IndelRealigner and putative PCR duplicates were flagged using MarkDuplicates, similar to the process described for the modern genomes. We also limited the incorporation of damaged sites into the variant-calling pipeline by hard-masking all sites that would be potentially affected by the characteristic ancient DNA patterns of cytosine deamination in single stranded overhangs. This mask was applied to 10 nucleotides on both ends of the merged reads from the ancient samples.  1.19), which was applied with "-C50" to adjust the mapping quality of the reads with multiple mismatches. We did not call differences in regions where the reference base was unknown, and the calls were limited to regions that were covered at least 4 times, and at most 350 times by the sequences in these samples. We then identified homozygous nonsynonymous substitutions unique (private) to each elephant and mammoth genome; these homozygous nonsynonymous substitutions were used in downstream analyses.
The M25 genome is likely a composite of multiple mammoth individuals 6 , therefore, we did not include M25 in downstream analyses. However, because we are primarily interested in private nonsynonymous variants within each elephant and mammoth genome we were able advantage of the composite nature of the M25 genome by excluding homozygous nonsynonymous substitutions identified in the three Asian elephants and the three mammoths that were also observed in M25. The rationale for this filtering process is that any homozygous nonsynonymous substitution observed in either the Asian elephants or the three mammoths that is also observed in M25 is not truly a private variant. We thus identified 106 fixed amino acid substitutions in 99 genes in the Oimyakon mammoth, 162 fixed amino acid substitutions in 143 genes in the M4 mammoth, and 594 fixed amino acid substitutions in 525 genes in the Wrangel Island mammoth.