Excess of genomic defects in a woolly mammoth on Wrangel island.

Woolly mammoths (Mammuthus primigenius) populated Siberia, Beringia, and North America during the Pleistocene and early Holocene. Recent breakthroughs in ancient DNA sequencing have allowed for complete genome sequencing for two specimens of woolly mammoths (Palkopoulou et al. 2015). One mammoth specimen is from a mainland population 45,000 years ago when mammoths were plentiful. The second, a 4300 yr old specimen, is derived from an isolated population on Wrangel island where mammoths subsisted with small effective population size more than 43-fold lower than previous populations. These extreme differences in effective population size offer a rare opportunity to test nearly neutral models of genome architecture evolution within a single species. Using these previously published mammoth sequences, we identify deletions, retrogenes, and non-functionalizing point mutations. In the Wrangel island mammoth, we identify a greater number of deletions, a larger proportion of deletions affecting gene sequences, a greater number of candidate retrogenes, and an increased number of premature stop codons. This accumulation of detrimental mutations is consistent with genomic meltdown in response to low effective population sizes in the dwindling mammoth population on Wrangel island. In addition, we observe high rates of loss of olfactory receptors and urinary proteins, either because these loci are non-essential or because they were favored by divergent selective pressures in island environments. Finally, at the locus of FOXQ1 we observe two independent loss-of-function mutations, which would confer a satin coat phenotype in this island woolly mammoth.

Introduction

Woolly mammoths (Mammuthus primigenius) were among the most populous large herbivores in North America, Siberia, and Beringia during the Pleistocene and early Holocene [1]. However warming climates and human predation led to extinction on the mainland roughly 10,000 years ago [2]. Lone isolated island populations persisted out of human reach until roughly 3,700 years ago when the species finally went extinct [3]. Recently, two complete high-quality high-coverage genomes were produced for two woolly mammoths [4]. One specimen is derived from the Siberian mainland at Oimyakon, dated to 45,000 years ago [4]. This sample comes from a time when mammoth populations were plentiful, with estimated effective population size of N = 13,000 individuals [4]. The second specimen is from Wrangel Island off the north Siberian coast [4]. This sample from 4,300 years ago represents one of the last known mammoth specimens. This individual comes from a small population estimated to contain roughly 300 individuals [4]. These two specimens offer the rare chance to explore the ways the genome responds to pre-extinction population dynamics.

Nearly neutral theories of genome evolution predict that small population sizes will lead to an accumulation of detrimental variation in the genome [5]. Such explanations have previously been invoked to explain genome content and genome size differences across multiple species [6]. Yet, within-species comparisons of how genomes are changed by small effective population sizes remain necessarily rare. These mammoth specimens offer the unique opportunity for within-species comparative genomics under a 43-fold reduction in population size. This comparison offers a major advantage as it will be free from confounding biological variables that are present in cross species comparisons. If nearly neutral dynamics lead to an excess of detrimental variation, we should observe an excess of harmful mutations in pre-extinction mammoths from Wrangel Island.

We use these two ancient DNA sequences to identify retrogenes, deletions, premature stop codons, and point mutations found in the Wrangel Island and Oimyakon mammoths. We identify an excess of putatively detrimental mutations, with an excess of stop codons, an excess of deletions, an increase in the proportion of deletions affecting gene sequences, an increase in non-synonymous substitutions relative to synonymous substitutions, and an excess of retrogenes, reflecting increased transposable element activity. These data bear the signature of genomic meltdown in small populations, consistent with nearly-neutral genome evolution. They furthermore suggest large numbers of detrimental variants collecting in pre-extinction genomes, a warning for continued efforts to protect current endangered species with small population sizes.

Results

Excess of amino acid substitutions and stop codons

We identified all SNPs in each mammoth genome as well as one Indian elephant specimen, Maya, using GATK [7]. We identified all non-synonymous and synonymous changes relative to the L. africana reference genome (https://www.broadinstitute.org/scientific-community/science/projects/mammals-models/elephant/elephant-genome-project) using r3.7 annotations lifted over to L. africana 4.0 genome sequences. We observe a significant increase in the number of heterozygous non-synonymous changes relative to synonymous changes in the Wrangel island genome compared with Oimyakon (χ2 = 68.799, df = 1, P < 2.2 × 10−16; S1 Table). There is also a significant increase in the number of homozygous mutations at non-synonymous sites relative to synonymous sites (χ2 = 9.96, df = 1, P < 0.0016; S1 Table). We further observe an excess of premature stop codons in the genome of the Wrangel Island mammoth, with 1.8X as many genes affected. There are 503 premature stop codons in the Oimyakon genome (adjusting for a 30% false negative rate at heterozygous sites) compared with 819 in the Wrangel island genome (Fig 1, Table 1). There are 318 genes that have premature stop codons that are shared across the two mammoths, and 357 genes that are truncated in both mammoths, including mutations that form at independent sites. A total of 120 of these genes have stop codons in the two mammoths as well as in Maya the Indian elephant, suggesting read through in the L. africana reference. Among truncated genes, there is a significant excess of olfactory genes and oderant binding receptors that appear to be pseudogenized with an EASE enrichment score of 9.1 (S2 Table) [8, 9]. We observe 85 truncated olfactory receptors and 3 vomeronasal receptors as well as multiple signal transduction peptides compared with 44 olfactory receptors and 2 vomeronasal receptors pseudogenized in the mainland mammoth.

Fig 1

Excess of putatively detrimental mutations in the Wrangel island genome.

A) Deletions B)Genes deleted C) Retrogenes D) Premature stop codons. Numbers shown are corrected for false negative rates of 30% for heterozygous SNPs and 0.5% for deletions in the lower coverage Oimyakon mammoth.

Table 1

Mutations identified in mammoth genomes.

Mutation	Oimyakon	Wrangel	Maya
Deletions	213461	27228	28095
Retrogenes	2130	2853	1575
Genes with exons deleted	11151	1628	3427
Stop Codons	5032	819	450
Stop Codons, excluding damage	377	645	401

1 Corrected for false negative rate of 0.5% in heterozygotes

2 Corrected for false negative rate of 30% at heterozygous sites established by Palkopoulou et al 2015.

It is possible that DNA damage in the archaic specimens could contribute to a portion of the observed stop codons. When we exclude A/G and C/T mutations, there is still a gross excess of premature stop codons, with 645 genes truncated in the Wrangel Island mammoth compared with 377 in the Oimyakon mammoth. Hence, the patterns are not explained solely by differential DNA damage in the two mammoths. Maya, the Indian Elephant specimen shows 450 premature stop codons, but 401 when A/G and T/C mutations are excluded. When putative damage to ancient DNA is excluded, Maya appears to house an intermediate number of premature stop codons, with a 6% increase compared to the Oimyakon mammoth.

Deletions

We identify 27228 deletions over 1 kb long in the Wrangel island genome, and 21346 (correcting for a 0.5% false negative rate at heterozygous sites) in the Oimyakon genome (Table 1). There are 6147 deletions (23%) identified in the Wrangel Island mammoth that are homozygous (≤ 10% coverage) compared with 5035 (24%) in the Oimyakon mammoth. (S3 Table). A total of 13,459 deletions are identified in both mammoth genomes (S4 Table). Some 4813 deletions in the Wrangel Island mammoth and 4598 in the Oimyakon mammoth appear hemizygous but have stretches of zero coverage for at least 50% of their length. These sites may represent multiple independent homozygous deletions that cannot be differentiated via change point statistics. Alternatively, they might indicate smaller secondary deletions that appear on hemizygous haplotypes. Such secondary deletions are common when large loop mismatch repair attacks unpaired, hemizygous stretches of DNA [10, 11]. The Wrangel Island Mammoth has sharply increased heterozygosity for deletions in comparison with the Oimyakon mammoth (S3 Table). Some portion of the inflated heterozygosity for deletions in the Wrangel Island mammoth could be due to this difficulty in inferring genotypes in a high throughput setting. Alternatively, the effective mutation rate may have increased as fewer deletions were removed from the population via purifying selection, inflating θ. It is also possible that there was an increase in the rate of deletions in the Wrangel Island lineage due to defective DNA repair mechanisms. An increase in non-homologous end joining after DNA breaks rather than double stranded break repair could putatively induce such a change in the deletion rate.

Maya the Indian elephant shows a larger number of deletions than the Oimyakon mammoth, but with different character from the Wrangel Island mammoth. The bulk of these are derived from 22,954 hemizygous deletions (S3 Table). Maya houses only 5141 homozygous deletions, similar to the mainland mammoth (S3 Table). There is an increase in the number of hemizygous deletions that affect gene sequences, but only a modest increase in the number of homozygous deletions that affect gene sequences (S3 Table). Competing pressures of higher N, longer time frames to accumulate mutations toward equilibrium frequencies, differences in mutation rates between the mammoths and elephants, differences in selective pressures, differences in the distribution of selective coefficients for deletions, different effective mutation rates due to different selective constraints, or differences in dominance coefficients might all contribute to differences in the number of deletions observed in elephants and mammoths. Additional samples would be necessary to determine the extent to which genetic declines may be influencing the diversity of deletions in modern Indian elephants. We currently have no basis for conclusions given this single sample, with no prior comparison.

There is a significant difference in the size distribution of deletions identified in the two mammoth samples, with a mean of 1707 bp in Oimyakon and 1606 bp in the Wrangel mammoth (Wilcox W = 304430000, P < 2.2e − 16; Fig 2). This difference could reflect either differences in DNA replication or repair mechanisms in the two mammoths, or altered selective constraints for different types of duplications. No significant difference is observed between the Wrangel island mammoth down sampled sequence data (W = 2004400, P = 0.3917) suggesting that the observed decrease in size is not due to differences in coverage. Some 1628 genes have deleted exons in the Wrangel Island mammoth compared to 1110 in Oimyakon (Table 1), a significant excess of genes deleted compared to expectations based on the number of deletions (χ2 = 12.717, df = 1,P = 0.0003623). Among these deleted genes, 112 in the mainland mammoth are homozygous compared to 133 homozygous exon deletions in the Wrangel Island Mammoth. Gene functions for affected genes in the Oimyakon mammoth include synapse functions, PHD domains, zinc fingers, aldo-keto metabolism, calcium dependent membrane targeting, DNA repair, transcription regulation, and development (S5 Table). Gene functions overrepresented among deletions in the Wrangel Island mammoth include major urinary proteins, lipocalins, and pheromones, pleckstrins, transcription regulation, cell transport, DNA repair, chromatin regulation, hox domains, and development (S5 Table).

Fig 2

eCDF for the size distribution of deletions in the Oimyakon and Wrangel island genomes.

There is a significant reduction in the size of deletions identified in the Wrangel Island Genome.

Among the genes deleted in the Wrangel Island mammoth, several have phenotypes of interest in other organisms. We observe a hemizygous deletion in riboflavin kinase RFK in the Wrangel Island mammoth, but normal coverage in the Oimyakon mainland mammoth (S1 Fig). Homozygous knockouts of riboflavin kinase, essential for B2 utilization/FAD synthesis, are embryonic lethal in mice [12]. Finally, we identify a hemizygous deletion in the Wrangel island mammoth that would remove the entire gene sequence at the FOXQ1 locus (S2 Fig). The alternative haplotype carries a frameshift mutation that disrupts the FOXQ1 functional domain. FOXQ1 knock-outs in mice are associated with the satin coat phenotype, which results in translucent fur but normal pigmentation due to abnormal development of the inner medulla of hairs [13], with two independent mutations producing this phenotype [13]. FOXQ1 also regulates mucin secretion in the GI tract, a case of pleiotropic functions from a single gene [14]. If the phenotype in elephantids matches the phenotype exhibited in mice, this mammoth would have translucent hairs and a shiny satin coat, caused by two independently formed knock-out alleles at the same locus. These genes each have functions that are conserved across mammals, though there is no guarantee that they would produce identical phenotypes in other species.

Retrogene formation

Retrogene formation can serve as a proxy for retrotransposon activity. We identify retrogenes that display exon-exon junction reads in genomic DNA. We observe 1.3X more retrogenes formed in the Wrangel island mammoth. The Wrangel Island mammoth has 2853 candidate retrogenes, in comparison with 2130 in the Oimyakon mammoth and 1575 in Maya (Table 1). There are 436 retrogenes that are shared between the two mammoths, though some of these could arise via independent mutations. This excess of retrogenes is consistent with increased retroelement activity in the Wrangel Island lineage. During retrogene formation, highly expressed genes, especially those expressed in the germline, are expected to contribute to new retrogenes. To determine the types of loci that had been copied by retrotransposons, we performed a gene ontology analysis using DAVID [8, 9]. Functional categories overrepresented among candidate retrogenes include genes involved in transcription, translation, cell division/cytoskeleton, post translational modification, ubiquitination, and chaperones for protein folding (S6 and S7 Tables). All of these are expected to be highly expressed during cell divisions or constitutively expressed, consistent with expectations that highly expressed genes will be overrepresented. Gene ontologies represented are similar for both mammoths (S6 and S7 Tables). Although these retrogenes are unlikely to be detrimental in and of themselves, they may point to a burst of transposable element activity in the lineage that led to the Wrangel island individual. Such a burst of TE activity would be expected to have detrimental consequences, additionally contributing to genomic decline.

Genomic effects of demography

Under nearly-neutral theory of genome evolution, detrimental mutations should accumulate in small populations as selection becomes less efficient [5]. This increase in non-neutral amino acid changes and premature stop codons is consistent with reduced efficacy of selection in small populations. We attempted to determine whether the data is consistent with this nearly-neutral theory at silent and amino acid replacement substitutions whose mutation rates and selection coefficients are well estimated in the literature. Under nearly neutral theory, population level variation for non-synonymous amino acid changes should accelerate toward parity with population level variation at synonymous sites.

Given the decreased population size on Wrangel Island, we expect to observe an accumulation of detrimental changes that would increase heterozygosity at non-synonymous sites (H) relative to synonymous sites (H) in the island mammoth. Heterozygosity depends directly on effective population sizes. We observe H = 0.00130 ± 0.00002 in the Wrangel Island mammoth, which is 80% of H = 0.00161 ± 0.00002 observed in the Oimyakon mammoth (Table 2). The magnitude of the difference between H in these two mammoths is 28 standard deviations apart, suggesting that these two mammoths could not have come from populations with the same effective population sizes. The specimens are well beyond the limits of expected segregating variation for a single population. To determine whether such results are consistent with theory, we fitted a model using PSMC [42] inferred population sizes for the Wrangel island mammoth, based on decay of heterozygosity of (1 − 1/2N) H0. The observed reduction in heterozygosity is directly consistent theoretical expectations that decreased effective population sizes would lower heterozygosity to H = 0.00131.

Table 2

Non-synonymous and synonymous heterozygosity.

	Wrangel	Oimyakon1
HS ± 2σ	0.00130 ± 0.00002	0.00161 ± 0.00002
HN ± 2σ	0.000490 ± 0.000012	0.000506 ± 0.000012
HN/HS	0.370	0.314

1 Oimyakon corrected for false negative rate of 30% established by Palkopoulou et al 2015.

At non-synonymous sites, however, there are no closed-form solutions for how H would decay under reduced population sizes. We observe H = 0.000490 in the Wrangel Island Mammoth, 95% of H = 0.000506 in the Oimyakon mammoth (Table 2). To determine whether such results could be caused by accumulation of nearly-neutral variation, we simulated population trajectories estimated using PSMC [42]. This trajectory shows ancient populations with N ≈ 104, followed by a population decline prior to extinction. These numbers are slightly lower than previous estimates of ancestral N based on mitochondrial DNA [43]. We were able to qualitatively confirm results that population trajectories from PSMC with previously described mutation rates and selection coefficients can lead to an accumulation of detrimental alleles in populations. However, the magnitude of the effects is difficult to fit precisely. The simulations show a mean H = 0.00148 and H = 0.000339 in Oimyakon and H = 0.00126 and H = 0.000295 for the Wrangel Island Mammoth (S3 Fig). In simulations, we estimate H/H = 0.229 both for the Oimyakon mammoth and directly after the bottleneck, but H/H = 0.233 in the Wrangel Island Mammoth at the time of the Wrangel Island mammoth. These numbers are less than empirical observations of H/H = 0.370 (Table 2). Several possibilities might explain the observed disparity between precise estimates from simulations versus the data. The simulations may be particularly sensitive to perturbations from PSMC population levels or time intervals. Similarly, selection coefficients that differ from the gamma distribution previously estimated for humans might lead to greater or lesser changes in small populations. Additionally, an acceleration in generation time on Wrangel Island is conceivable, especially given the reduced size of Wrangel Island mammoths [15]. Finally, positive selection altering nucleotide variation on the island or the mainland could influence diversity levels.

Founder effects during island invasion sometimes alter genetic diversity in populations. However, it is unlikely that a bottleneck alone could cause an increase in H/H. There is no evidence in effective population sizes inferred using PSMC to suggest a strong bottleneck during Island colonization [4]. The power of such genetic analyses may be limited, but these results are in agreement with paleontological evidence showing no phenotypic differentiation from the mainland around 12,000 years ago followed by island dwarfism much later [15]. During glacial maxima, the island was fully connected to the mainland, becoming cut off as ice melted and sea levels rose. The timing of separation between the island and mainland lies between 10,000 years and 14,000 years before present [3, 15–17], but strontium isotope data for mammoth fossils suggests full isolation of island populations was not complete until 10,000-10,500 years ago [18]. Forward simulations suggest that hundreds of generations at small N are required for detrimental mutations to appear and accumulate in the population. These results are consistent with recent theory suggesting extended bottlenecks are required to diminish population fitness [19]. Thus, we suggest that a bottleneck alone could not produce the accumulation of H/H that we observe.

E. maximus indicus specimen, Maya shows an independent population decline in the past 100,000 years, with current estimates of N = 1000 individuals (S4 Fig). This specimen shows a parallel case of declining population sizes in a similar species of elephantid. Maya houses hemizygous deletions in similar numbers with the Wrangel Island Mammoth. However, the number of stop codons and homozygous deletions is intermediate in comparison with the Oimyakon and Wrangel mammoths (Table 1). It is possible that Indian elephants, with their recently reduced population sizes may be subject to similar accumulation of detrimental mutations, a prospect that would need to be more fully addressed in the future using population genomic samples for multiple individuals or timepoints and more thorough analyses.

Discussion

Nearly neutral theories of genome evolution

Nearly-neutral theories of genome evolution have attempted to explain the accumulation of genome architecture changes across taxa [5]. Under such models, mutations with selection coefficients less than the nearly neutral threshold will accumulate in genomes over time. Here, we test this hypothesis using data from a woolly mammoth sample from just prior to extinction. We observe an excess of retrogenes, deletions, amino acid substitutions, and premature stop codons in woolly mammoths on Wrangel Island. Given the long period of isolation and extreme population sizes observed in pre-extinction mammoths on Wrangel Island, it is expected that genomes would deteriorate over time. These results offer genetic support for the nearly-neutral theory of genome evolution, that under small effective population sizes, detrimental mutations can accumulate in genomes. Independent analysis supporting a reduction in nucleotide diversity across multiple individuals at MHC loci suggests a loss of balancing selection, further support the hypothesis that detrimental variants accumulated in small populations [20].

We observe two independent loss-of-function mutations in the Wrangel Island mammoth at the locus of FOXQ1. One mutation removes the entire gene sequence via a deletion, while the other produces a frameshift in the CDS. Based on phenotypes observed in mouse models, these two independent mutations would result in a satin fur coat, as well as gastric irritation [14]. Many phenotypic screens search for homozygous mutations as causative genetic variants that could produce disease. More recently, it has been proposed that the causative genetic variation for disease phenotypes may be heterozygous non-complementing detrimental mutations [21]. These data offer one case study of independent non-functionalizing mutations in a single individual, genetic support for independent non-functionalizing mutations at a single locus. Woolly mammoth outer hairs house multiple medullae, creating a stiff outer coat that may have protected animals from cold climates [22] (though see [23] for alternative interpretations). Putative loss of these medullae through loss of FOXQ1 could compromise this adaptation, leading to lower fitness.

Island specific changes

One of the two specimens comes from Wrangel Island, off the northern coast of Siberia. This mammoth population had been separated from the mainland population for at least 6000 years after all mainland mammoths had died off. Prior to extinction, some level of geographic differentiation combined with differing selective pressures led to phenotypic differentiation on Wrangel island [15]. Island mammoths had diminished size, but not until 12,000 years ago when mainland populations had reduced and ice sheets melted [15]. One possible explanation for the poor fit of simulations is that generation time may have decreased. Previous work suggested a very high mutation rate for woolly mammoths based on comparisons between island and mainland mammoths. It is possible that an acceleration in generation times could cause the accumulation of more mutations over time, and that the real mutation rate is similar to humans (1 − 2 × 10−8 [24] rather than 3.8 × 10−8 [4]). Such changes would be consistent with island dwarfism being correlated with shorter generation times, and would explain the unusually high mutation rate estimate for mammoths based on branch shortening observed in [4].

We observe large numbers of pseudogenized olfactory receptors in the Island mammoth. Olfactory receptors evolve rapidly in many mammals, with high rates of gain and loss [25]. The Wrangel island mammoth has massive excess even compared to the mainland mammoth. Wrangel island had different flora compared to the mainland, with peat and sedges rather than grasslands that characterized the mainland [17]. The island also lacked large predators present on the mainland. It is possible that island habitats created new selective pressures that resulted in selection against some olfactory receptors. Such evolutionary change would echo gain and loss of olfactory receptors in island Drosophila [26]. In parallel, we observe a large number of deletions in major urinary proteins in the island mammoth. In Indian elephants E. maximus indicus, urinary proteins and pheromones ellicit behavioral responses including mate choice and social status [27]. It is possible that coevolution between urinary proteins, olfactory receptors, and vomeronasal receptors led to a feedback loop, allowing for rapid loss in these related genes. It is equally possible that urinary peptides and olfactory receptors are not essential and as such they are more likely to fall within the nearly neutral range [25]. Either of these hypotheses could explain the current data.

Implications for conservation genetics

Many factors contributed to the demise of woolly mammoths in prehistoric times. Climate change led to receding grasslands as forests grew in Beringia and North America and human predation placed a strain on already struggling populations [2]. Unlike many cases of island invasion, Wrangel Island mammoths would not have continuous migration to replenish variation after mainland populations went extinct. Under such circumstances, detrimental variation would quickly accumulate on the island. The putatively detrimental variation observed in these island mammoths, with the excess of deletions, especially recessive lethals may also have limited survival of these struggling pre-extinction populations. Climate change created major limitations for mammoths on other islands [28], and these mammoths may have struggled to overcome similar selective pressures.

Many modern day species, including elephants, are threatened or endangered. Asiatic cheetahs are estimated to have fewer than 100 individuals in the wild [29]. Pandas are estimated to have 1600 individuals living in highly fragmented territories [30]. Mountain Gorilla population census sizes have been estimated as roughly 300 individuals, similar to effective population sizes for pre-extinction mammoths [31]. If nearly neutral dynamics of genome evolution affect contemporary endangered species, detrimental variation would be expected in these genomes. With single nucleotide changes, recovered populations can purge detrimental variation in hundreds to thousands of generations, returning to normal genetic loads [19]. However, with deletions that become fixed in populations, it is difficult to see how genomes could recover quickly. The realm of back mutations to reproduce deleted gene sequences will be limited or impossible. Although compensatory mutations might conceivably correct for some detrimental mutations, with small effective population sizes, adaptation through both new mutation and standing variation may be severely limited [32]. Thus we might expect genomes affected by genomic meltdown to show lasting repercussions that will impede population recovery.

Ethics statement

All sequences are taken from publicly available sequence data in the ENA or SRA. Indian elephant specimens for previously published sequence data were handled by the San Diego Zoo.

Materials and methods

Genome sequence alignments

We used previously aligned bam files from ERR852028 (Oimyakon, 11X) and ERR855944 (Wrangel, 17X) (S8 Table) [4] aligned against the L. africana 4.0 reference genome (available on request from the Broad Institute—vertebrategenomes@broadinstitute.org; https://www.broadinstitute.org/scientific-community/science/projects/mammals-models/elephant/elephant-genome-project). We also aligned 33X coverage of sequencing reads for one modern E. maximus indicus genome Maya (previously described as “Uno”) using bwa 0.7.12-r1044 [33], with parameters set according to [4] bwa aln -l 16500 -o 2 -n 0.01. The E. maximus indicus sample, previously labeled in the SRA as “Uno”, is from Maya, a former resident of the San Diego Zoo wild-born in Assam, India, North American Studbook Number 223, Local ID #141002 (O. Ryder, personal communication). We were not able to use two other mammoth sequences are publicly available, M4 and M25 from Lynch et al. [34]. These sequences display abnormal PSMC results (S4 Fig), high heterozygosity (S5 Fig), and many SNPs with asymmetrical read support (S6 Fig). The unrealistically high heterozygosity as well as abnormal heterozygote calls raise concerns with respect to sequence quality. For further description, please see Supporting Information.

Synonymous and nonsynonymous substitutions

We used the GATK pipleine [7] v3.4-0-g7e26428 to identify SNPs in the aligned sequence files for the Oimyakon and Wrangel Island mammoths. We identified and realigned all indel spanning reads according to the standard GATK pipeline. We then identified all SNPs using the Unified Genotyper, with output mode set to emit all sites. We used all CDS annotations from cDNA annotations from L. africana r3.7 with liftover coordinates provided for L. africana 4.0 to identify SNPs within coding sequences. We identified all stop codons, synonymous substitutions, and non-synonymous substitutions for the Wrangel Island and Oimyakon mammoths at heterozygous and homozygous sites.

Retrogenes

We aligned all reads from the mammoth genome sequencing projects ERR852028 (Oimyakon) and ERR855944 (Wrangel) (S8 Table) against elephant cDNA annotations from L. africana r3.7. Sequences were aligned using bwa 0.7.12-r1044 [33], with parameters set according to [4] bwa aln -l 16500 -o 2 -n 0.01 in order to account for alignments of damaged ancient DNA. We then collected all reads that map to exon-exon boundaries with at least 10 bp of overhang. Reads were then filtered against aligned genomic bam files produced by Palkopoulou et al [4], discarding all exon-exon junction reads that have an alignment with equal or better alignments in the genomic DNA file. We then retained all putative retrogenes that showed signs of loss for two or more introns, using only cases with 3 or more exon-exon junction reads.

We calculated coverage depth using samtools [35] with a quality cutoff of -q 20. We then implemented change point analysis [36] in 20 kb windows. Change point methods have been commonly used to analyze microarray data and single read data for CNVs [37-39] The method seeks compares the difference in the log of sum of the squares of the residuals with one regression line vs. two regression lines [36]. The test statistic follows a chi-squared distribution with a number of degrees of freedom determined by the number of change-points in the data, in this case df = 1. We required significance at a Bonferroni corrected p-value of 0.05 or less. We allowed for a maximum of one CNV tract per window, with minimum of 1 kb and maximum of 10 kb (half the window size) with a 100 bp step size. We did not attempt to identify deletions smaller than 1 kb due to general concerns of ancient DNA sequence quality, limitations to assess small deletions in the face of stochastic coverage variation, and concerns that genotype calls for smaller deletions might not be as robust to differences in coverage between the two mammoths. Sequences with ‘N’s in the reference genome did not contribute to change point detection. We excluded all deletions that were identified as homozygous mutations in both mammoths and in E. maximus indicus specimen Maya, as these suggest insertion in the L. africana reference rather than deletion in other elephantids. To determine the effects that coverage differences would have on deletions, we downsampled the sequence file for the Wrangel Island mammoth using samtools to 11X coverage, using chromosome 1 as a test set. We observe a reduction in the number of deletions for chromosome 1 from 1035 deletions to 999 deletions, resulting in an estimated false negative rate of 0.5% at reduced coverage for deletions greater than 1 kb. Highly diverged haplotypes with greater than 2% divergence might prevent read mapping and mimic effects of deletions, but this would require divergence times within a species that are greater than the divergence between mammoths and L. africana. Mutations were considered homozygous if mean coverage for the region was less than 10% of the background coverage level. Otherwise it was considered to be heterozygous. These methods are high-throughput, and it is possible that multiple small homozygous deletions interspersed with full coverage sequences might mimic heterozygote calls. Whether such mutations might meet the conditions for significant change-point detection would depend on the deletion length, placement, and background coverage level.

Demography

We identified SNPs that differentiate Mammoth genomes from the reference using samtools mpileup (options -C50 -q30 -Q30), and bcftools 1.2 consenus caller (bcftools call -c). The resulting vcf was converted to fastq file using bcftools vcf2fq.pl with a mimimum depth of 3 reads and a maximum depth of twice the mean coverage for each genome. Sequences were then converted to psmc fasta format using fq2psmcfa provided by psmc 0.6.5-r67. We then ran psmc with 25 iterations (-N25), an initial ratio of θ/ρ of 5 (-r5), and parameters 64 atomic time intervals and 28 free parameters (-p “4+25*2+4+6”) as was done in previous analysis of woolly mammoths [4]. Effective population sizes and coalescence times were rescaled using previously estimated mutation rates of 3.8 × 10−8. Using the population size estimates from PSMC, we calculated the expected reduction in heterozygosity at synonymous sites according to for each time period in PSMC output. We compared the number of deletions, number of premature stop codons, proportion affecting gene sequences, and number of putative retrogenes between the two mammoth genomes using chi squared tests.

Simulations

To determine expectations of sequence evolution at non-synonymous sites under population crash, we ran simulations using SLiM v. 2.0 population genetic software [40]. We modeled two classes of sites: neutral and detrimental. For detrimental mutations we used a gamma distributed DFE with a mean of -0.043 and a shape parameter of 0.23 as estimated for humans [41], assuming a dominance coefficient of 0.5 and free recombination across sites. Mutation rates were set as 3.8 × 10−8 based on previously published estimates [4]. The trajectory of population sizes was simulated according to estimates from PSMC, omitting the initial and final time points from PSMC, which are often subject to runaway behavior. We then simulated the accumulation of H/H in the Wrangel Island Mammoths. Simulations were run with a burn-in of 100,000 generations. We simulated 460 replicates of haplotypes with 100 sites for each mutation class.

Gene ontology

To gather a portrait of functional categories captured by deletions, retrogenes, and stop codons, we identified all mouse orthologs based on ENSEMBL annotations for L. africana 3.7 for affected gene sequences. We then used DAVID gene ontology analysis with the clustering threshold set to ‘Low’ (http://david.ncifcrf.gov/; Accessed April 2016) [8, 9]. S2–S7 Tables include all functions overrepresented at an EASE enrichment cutoff of 2.0. Full gene ontology data is included in Supplementary Information.

Supporting text.

(PDF)

Click here for additional data file.

Non synonymous and synonymous changes in Wrangel and Oimyakon mammoths.

(PDF)

Click here for additional data file.

Gene ontology categories for premature stop codons in mammoths.

(PDF)

Click here for additional data file.

Deletions identified in mammoth genomes.

(PDF)

Click here for additional data file.

Shared deletions in mammoth genomes.

(PDF)

Click here for additional data file.

Gene ontology for deleted exons.

(PDF)

Click here for additional data file.

Gene ontology for retrogenes in the Oimyakon mammoth.

(PDF)

Click here for additional data file.

Gene ontology for retrogenes in the Wrangel island mammoth.

(PDF)

Click here for additional data file.

SRA and ENA identifiers for mammoth and Elephant sequence data.

(PDF)

Click here for additional data file.

Heterozygous sites per 10 kb.

(PDF)

Click here for additional data file.

Asymmetrical support.

(PDF)

Click here for additional data file.

Coverage depth at the RFK locus.

(PDF)

Click here for additional data file.

Coverage depth at the FOXQ1 locus.

(PDF)

Click here for additional data file.

Simulations for heterozygosity at synonymous and non-synonymous sites.

(PDF)

Click here for additional data file.

PSMC results.

(PDF)

Click here for additional data file.

Heterozygosity for mammoth and elephant samples.

(PDF)

Click here for additional data file.

Asymmetric SNPs.

(PDF)

Click here for additional data file.

Asymmetric SNPs, excluding damage.

(PDF)

Click here for additional data file.

SuppFiles.zip—Data release archive.

(ZIP)

Click here for additional data file.