Geometric morphometric analyses

Taxon sampling

Taxon sampling was largely based on a previous study6, which itself was an expansion of the 3D cranial dataset assembled by Wroe & Milne5 to examine morphological convergence in carnivorous eutherians (order Carnivora) and marsupials, including the thylacine (Thylacinus cynocephalus). The Carnivora is divided into two suborders, the Feliformia and Caniformia, which together comprise over 280 living species44. Among these are the well-known large-bodied carnivores such as lions, wolves and hyenas, but also omnivores (some bears), insectivores (aardwolf and bat-eared fox) and a diversity of small- to medium-sized mustelids, raccoons, civets, mongooses and others with varying ecologies. Members of the other major clade of carnivorous eutherians, the extinct Creodonta, ranged from cat- to hyena-sized and were the most successful and abundant predators of the early Tertiary period, an ecological niche later dominated by carnivorans45. Among metatherians, carnivory (here defined as the killing and comminuting of vertebrate prey, sensu46) has evolved independently in the extinct stem order Sparassodonta (for example, Thylacosmilus atrox) from South America, the extinct Australian diprotodontian family Thylacoleonidae, and in members of the Australian order Dasyuromorphia.

A previous study5 sampled crania of 28 extant and 2 fossil species from 7 families of Carnivora, covering a wide range of body sizes (3–382 kg) and feeding ecologies (hypercarnivorous, bone-cracking, omnivorous and insectivorous), to which they compared crania of 13 marsupial species from 7 families and 4 orders. To match morphological and ecological variation in carnivorans, they included the insectivorous numbat Myrmecobius fasciatus and other generalists outside of Dasyuromorphia, plus five extinct marsupials: the thylacine and Miocene Nimbacinus dicksoni (Thylacinidae); the highly specialized predator and largest carnivorous Australian mammal Thylacoleo carnifex (comparable to a large felid) and its smaller more primitive relative Wakaleo vanderleuri, both thylacoleonids; and the most morphologically conservative dasyurid known, Barinya wangala 5. A subsequent study6 extended this dataset by an additional 35 species, sampling from members of the stem clades Creodonta and Sparrasodonta, as well as a broad diversity of extinct carnivorans lacking modern counterparts (such as the ‘false’ sabre-toothed cats Nimravidae and bear-dogs of Amphicyonidae), covering nearly the full morphological and ecological range of carnivorous terrestrial mammals.

Here we further expand that 3D cranial dataset6 by an additional 35 species, including members of Canis and Vulpes (the eutherian genera to which the thylacine is most commonly compared46,47,48,49,50,51,52,53), dasyuromorphs and New World didelphids, 8 families of Diprotodontia (the clade containing carnivorous thylacoleonids but also kangaroos, wallabies, koalas, possums, gliders, and wombats), and other insectivorous australidelphids from the orders Microbiotheria, Notoryctemorphia and Peramelemorphia. The majority of newly added taxa are relatively small (0.1–5 kg body mass), representing the non-carnivorous marsupials to which the thylacine is most closely related (such as the southern marsupial mole Notoryctes typhlops and the eastern barred bandicoot Perameles gunnii). Herbivorous diprotodontids were included to increase metatherian morphospace to a level comparable to that of the sampled eutherians, since marsupial crania have been shown to be significantly less disparate in shape than for eutherians54. They were also sampled to improve ancestral state reconstructions of cranial shape for Diprotodontia, which were used in convergence tests between the thylacine and other carnivorous marsupials (such as T. carnifex; see below).

Cranial data acquisition

3D cranial landmark data for 80 extinct and extant mammals were acquired from that previously published dataset6. Additional landmarks were generated from 28 newly sampled museum specimens and 13 cranial surface models made available through DigiMorph.org (Supplementary Table 4). Surface scanning was performed at the School of Engineering, University of Melbourne, using a NextEngine 3D laser scanner (NextEngine, Santa Monica, California). X-ray computed tomography (CT) scans were made at the School of Earth Sciences, University of Melbourne, on a GE Phoenix Nanotom M and reconstructed in datos|x-reconstruction software (GE Sensing and Inspection Technologies GmbH, Wunstorf, Germany). All scanned specimens were adult males unless otherwise noted (Supplementary Table 4). 3D models of the surface scans were produced in Meshlab (Visual Computing Lab, ISTI, CNR) and in VGStudio Max 2.1 (Volume Graphics, Heidelberg, Germany) for CT data. Landmarks were placed on the newly generated models in Landmark Editor (Institute of Data Analysis and Visualization, UC Davis, USA), with landmark locations matching those in the previously published analysis6, so that the datasets could be combined. In total, thirty homologous landmarks were digitized on the dorsal and lateral surfaces of the cranium (Supplementary Fig. 5 and Supplementary Table 5). Because the goal of the original study was to examine cranial convergence related to diet and feeding ecology, landmarks were focused on facial, dental and zygomatic regions, including muscle attachment sites such as the sagittal crest6. However, we note that landmarks on the palatal surface are lacking, which may miss other aspects of cranial morphology related to diet. Preliminary analyses of a subset of our landmarked crania revealed that variation due to measurement error was negligible (<2%; results not shown), indicating that landmarks could be placed confidently on all specimens under study. The final combined data set consisted of 171 individuals from 113 species, including two representatives of the thylacine (Supplementary Table 4).

Geometric information was extracted from the landmark coordinates by a generalized Procrustes fit55. The resulting Procrustes coordinates representing the symmetric component of shape variation after translating, scaling and rotating all individuals to a common average were extracted and used as shape variables in all analyses. The asymmetric component, which is typically the focus of studies on developmental integration and modularity56, was disregarded. Individual size information, preserved in the centroid size, was calculated as the square root of the sum of squared distances of landmarks to the centre of their configuration. Species with multiple specimens in the dataset were represented by their mean Procrustes shape and centroid size (CS) values.

To examine evolutionary trends in cranial morphology, we generated a simplified phylogenetic tree of mammals using a direct supertree approach57. Published phylogenies of eutherian and metatherian species were pruned to match our taxon sampling and combined at the root using the bind.tree function in the R package Ape v3.558. Extinct taxa were placed according to published morphological analyses and, in some cases, phylogenetic analyses including ancient DNA (Supplementary Fig. 10). Although phylogenetic positions of many of the fossils are based on cranial material, and are therefore not strictly independent of the cranial shape data, the majority of coded cranial characters are highly detailed and involve qualitative descriptions of specific cranial elements not likely to be captured in our landmark data. A permutation test using a multivariate generalization of a previously published59 K statistic (K mult )60 revealed significant phylogenetic structure in our data (K mult shape = 0.531, P = 0.0001; K mult log CS = 0.306, P = 0.02); accordingly, we used phylogenetic comparative methods in subsequent tests.

Dominant patterns of cranial shape variation were identified by principal component analysis (PCA). To determine evolutionary patterns of cranial shape change, species PC scores were mapped onto the phylogeny and weighted-square change parsimony was used to reconstruct the morphological state of ancestral nodes61. The mapped tree was then projected back into morphospace to visualize patterns of phenotypic evolution. Inspection of the PCA plots, together with examination of relative changes in landmark positions along each axis aided in interpretations of evolutionary changes in cranial shape.

For statistical analyses of phenotypic variation, we constructed analysis of variance models for individual effects based on the Procrustes distance among species, known as Procrustes ANOVA62. Each species was assigned to a general diet type (carnivore, herbivore, insectivore, omnivore; Supplementary Table 4), and the relative effects of size (represented by centroid size) and diet on cranial shape were compared. Given the extensive range of body sizes in our data set, centroid size was log-transformed prior to analysis63. Due to the non-independence of related species64, tests were also run in a phylogenetic context using our mammalian supertree (Supplementary Fig. 10). A distance-based phylogenetic generalized least squares model (D-PGLS), equivalent to phylogenetically independent contrasts65, was constructed separately for each of the above effects. Shape data were permutated across the tips of the phylogeny, and trait values predicted under a Brownian motion (BM) model of evolution were compared to those observed to assess statistical significance. This model assumes that phenotypic changes are independent from time step to time step, and that phenotypic variation increases proportionally with time66. Although D-PGLS is currently limited to BM, comparisons of phenotypic patterns derived from BM and other processes such as the Ornstein-Uhlenbeck model67 may be more appropriate in the future65. A complete list of statistical results is given in Supplementary Table 6.

Finally, to determine the degree of evolutionary convergence among thylacine and extant canids in cranial morphospace, we calculated three recently developed distance-based measures (C 1 , C 2 and C 3 )23, which quantify the amount of ancestral phenotypic space between putatively convergent taxa that has been ‘closed’ by subsequent evolution. Under this approach, convergence is implied when two or more taxa have evolved to be more similar to one another than their ancestors were to each other23. Unlike other convergence methods like SURFACE68, which identify convergent shifts in selective regimes among lineages, C 1 -C 3 measure the increase in phenotypic similarity between taxa compared to that between the most divergent species in their lineages, without assuming an adaptive process. Therefore, the more dissimilar the ancestor and the more similar the descendants, the greater the strength of the convergence. C 1 23 is estimated as the inverse ratio of the maximum Procrustes distance (D) between lineages from their most recent common ancestor to the distance between their extant tips, measured as C 1 = 1 − D tip / D max . Values range from 0 to 1, with 0 indicating that lineages are as dissimilar as they have ever been, and 1 meaning complete evolutionary convergence, that is, descendants are indistinguishable. While C 1 is a proportion and thus scaled to allow comparisons between taxa, C 2 represents the absolute amount of evolution that has occurred during convergence, indicating the magnitude of change expressed as C 2 = D tip / D max . From C 2 , the proportion of convergence relative to the total evolution (L tot.lineage , the sum of all phenotypic distances from ancestors to descendants) along those lineages since their most recent common ancestor can also be estimated as C 3 = C 2 / L tot.lineage . For all measures, values closer to 1 indicate greater morphological similarity. Species scores from the first 31 PC axes, accounting for 99% of the total morphological variation, were used as phenotypic variables. C 1 –C 3 were estimated for all pairwise comparisons between the thylacine and species of Canis and Vulpes, the two eutherian lineages which the thylacine most superficially resembles46,47,50,51. We also estimated convergence values between the thylacine and its closest carnivorous marsupial relatives, the extinct Nimbacinus dicksoni and Barinya wangala, as well as with the closely related insectivorous taxon, Myrmecobius fasciatus (Fig. 2a), and other carnivorous dasyurids and diprotodontids (Supplementary Table 7). Statistical significance of each measure was determined by simulation of the parameters derived from the observed data on our phylogeny using a BM model of evolution23. As an independent assessment of phenetic relationships among taxa, we additionally conducted an unweighted pair group average (UPGMA) hierarchical cluster analysis on the matrix of Procrustes distances between species using PAST software v3.1669. The resulting UPGMA phenogram (Supplementary Fig. 11) places the thylacine within a carnivorous eutherian, predominantly canid cluster closest to the red fox (Vulpes vulpes).

All other analyses were performed using MorphoJ v1.06d55, except for K mult and ANOVA, which were estimated in the geomorph package of R70, and convergence measures C 1 –C 3 calculated in the R package convevol23. For all tests, statistical significance (α = 0.05) was determined by a random permutation procedure of 10,000 iterations, except in convevol, which was run for 1,000 iterations.

Sequencing and assembly of the thylacine genome

Sample collection and DNA extraction and sequencing

DNA was obtained from a >108-year-old ethanol-preserved pouch young specimen obtained from the Museums Victoria collection (specimen C5757; Melbourne, Victoria, Australia; Fig. 1c). The sampled individual was a female, approximately 1 month old at the time of death. She was part of a litter of four young, three of which, C5755 (male), C5756 (female) and C5757 (female) are preserved intact at Museums Victoria. The fourth pouch young, C5754 (male) was sectioned for microscopy in 199471. All four pouch young, along with their mother were collected dead on 23 June 1909.

Soft tissues were washed several times in 1× PBS before being subjected to a prolonged (overnight) proteinase K digestion. We extracted and purified genomic DNA using a method similar to that described previously72. In short, the digested tissue was mixed with an equal volume of phenol chloroform isoamyl alcohol 25:24:1 and mixed for >1 hour. Samples were then centrifuged and the aqueous phase removed to a new tube; this procedure was repeated until the interphase was clear. The final supernatant was purified and concentrated to a final volume of ~50 µl using the Millipore MRCPRT010 purification column. We obtained >1 µg of total genomic DNA from a single extraction of 100 mg of starting tissue. The isolated DNA was fragmented as expected73, and contained kilobase length fragments as determined by electrophoresis using the BioAnalyzer 2100 high sensitivity DNA assay (Agilent; Supplementary Fig. 1). Two paired-end libraries were prepared using the Illumina TruSeq Nano kits (Revision A, May 2013, initial release, part number 15041110) using the manufacturers protocols for 350 bp (library 1) and 550 bp (library 2) fragments libraries, with the only modification being a reduced number of PCR cycles (6 cycles instead of 8). In short, thylacine genomic DNA was sheared with a Covaris ultrasonicator. 100 ng of input gDNA was used to construct library 1 and 200 ng was used to construct library 2. Fragments were then end repaired to remove 5ʹ and 3ʹ overhangs. End-repaired fragments were then size selected using manufacturer-provided magnetic sample purification beads, 3ʹ ends were acetylated, and adaptors ligated to the 5ʹ and 3ʹ ends. The libraries were then amplified with 6 cycles of PCR (modified from 8 cycles). Finally, the resulting library size distribution was analysed using the BioAnalyzer 2100 HiSensitivity DNA assay (Supplementary Fig. 1). Libraries were sequenced on the Illumina NextSeq 5500 (~1.14 billion reads) at the University of Connecticut and HiSeq 2000 (~192 million reads) by Perkin Elmer, using 2 × 150 bp and 2 × 100 bp chemistry, respectively (Supplementary Table 1 and Supplementary Fig. 1).

Sequence data preprocessing and quality assessment

Data were filtered by PHRED quality score and length using Trimmomatic v0.3274 (parameters: ILLUMINACLIP:2:30:10, LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:50 AVGQUAL:20). Data were depleted of microbial and mitochondrial sequences by mapping to two databases (one containing all Ensembl microbial and fungal genomes and another containing the thylacine and human mitochondrial genomes) using bwa mem v.0.7.12 (parameters −B 8 −O 12,12). Reads that mapped to the contaminant database were removed from the dataset and not subjected to further analysis. To further validate that retained reads were composed largely of thylacine DNA, reads were mapped against the repeat-masked genomes of three extant marsupials: the Tasmanian devil (Sarcophilus harrisii) 10, the tammar wallaby (Macropus eugenii)11 and the gray short-tailed opossum (Monodelphis domestica)12 as well as the complete, unmasked genome of the Tasmanian devil (assembly versions given in Supplementary Table 8). The percentage of reads mapping uniquely to each reference genome was graphed (Supplementary Fig. 2). The vast majority of reads mapped to the unmasked assembly of the Tasmanian devil, the closest relative to the thylacine available, indicating very low levels of microbial contamination. The percentage of mapped reads was observed to decline with increasing phylogenetic distance between the reference genome and the thylacine, indicating that the majority of reads were thylacine in origin (Supplementary Fig. 2).

To assess DNA damage in the thylacine sample, Illumina reads were processed and mapped to the thylacine reference-based contigs (see next section) using the PALEOMIX pipeline75. Briefly, residual adaptor sequences were trimmed and reads shorter than 25 nucleotides were removed, and overlapping paired reads were collapsed using AdapterRemoval v.2.076. Mapping of collapsed and paired reads was performed using the mem algorithm with default parameters in bwa v.0.7.1377. Reads with mapping quality below 25 were removed. Finally, mapDamage 2.013 was used to assess levels of cytosine deamination and depurination on a random subset of 100,000 reads (default parameter) for each of the sequencing datasets. The only damage pattern characteristic of ancient DNA sequence data was a slight increase of depurination (G only) immediately before the reads start (Supplementary Fig. 3)78,79,80,81. Since we did not observe an increase in 5′ C-to-T substitutions and 3′ G-to-A substitutions (Supplementary Fig. 3), we could not apply the mapDamage2.0 model to estimate parameters such as single-stranded overhang length (λ) and deamination rates (δ S and δ D )13. Overall, the lack of characteristic patterns of DNA damage (Supplementary Fig. 3) and the very large DNA fragments (Supplementary Fig. 1) are signs of the exceptional preservation of the DNA.

De novo and reference-based genome assembly

Thylacine de novo contigs were assembled from the reads that survived filtering (see above; Sample collection and DNA extraction and sequencing) using the De Bruijn graph-based assembler, SOAPdenovo282 (parameters −K 73 −M 3 −D 4). The quality of de novo and reference-based contigs were assessed using BUSCO (v1.22)83 with maximum genomic regions set to 3 and the vertebrata profile provided at http://busco.ezlab.org/files/vertebrata_buscos.tar.gz (Supplementary Table 21). The thylacine de novo contigs were used for retrophylogenomic analyses (see below) but not for further comparative genomic analyses.

Thylacine read data were then re-mapped to the repeat-masked Tasmanian devil genome (Ensembl Devil_ref v7.0) using bwa mem v0.7.12 (parameters −B 3 −O 5,5). Alignments were processed and variants called using the GATK Best Practices pipeline (v3.4-46-gbc02625)84,85,86. A reference-based assembly was generated by producing a consensus sequence of thylacine read data using samtools (v0.1.19) mpileup (parameters −I −B −u, without inclusion of the reference genome fasta) and piping to bcftools (v0.1.19) view (parameters −I −cg) and vcfutils vcf2fq87. To preserve the coordinate system of the Tasmanian devil genome, referenced-based scaffolds were padded with ‘N’ characters. For comparative data, reference-based assemblies were also produced for five wild canids (wolf, coyote, golden jackal, red fox and arctic fox). Data for the wild canids was downloaded from the Sequence Read Archive (Supplementary Table 8). The same pipeline was used to produce the wild canid reference-based assemblies as was used for the thylacine, however data was mapped to the repeat-masked dog genome (Ensembl CanFam3.1)88.

Retrophylogenomics

Retroelement identification

WSINE1/1a and WALLSI1A are SINEs (short interspersed elements) mobilized by LINE1 (long interspersed elements) and RTEs (retrotransposon-like elements), respectively. These elements were previously selected to computationally screen the genome sequences of dasyuromorph species for phylogenetic markers, and returned 25 diagnostic markers15. The corresponding loci of these markers were used in the present study as queries for blastn (BLAST version 2.2.28+)89 searches of thylacine contigs. All but 3 of them had orthologous sequences in the thylacine.

Twelve of the investigated genomic loci carried orthologous SINEs (5 WALLSI1A and 7 WSINE1/1a insertions) in all analysed dasyuromorphs, including the thylacine and numbat, providing significant evidence for the placement of the thylacine within the monophyletic order Dasyuromorphia (insectivore–carnivores) (one-directional KKSC insertion significance test16 for 12:0 markers, P < 2.1 × 10−4; Fig. 2a and Supplementary Fig. 4). Nine of the 25 diagnostic WSINE1/1a insertions (plus an additional 3 markers discovered during bioinformatics and experimental screens, markers 41–43; Supplementary Fig. 4) were absent in orthologous loci of the thylacine and numbat but were present in all investigated Dasyuridae, including in dunnart (Sminthopsis spp.), quoll (Dasyurus geoffroii) and the Tasmanian devil, verifying the monophyly of this family to the exclusion of the thylacine and numbat (one-directional KKSC for 12:0 markers, P < 2.1 × 10−4; Fig. 2a and Supplementary Fig. 4). Finally, one of the diagnostic WSINE1s was present in all Dasyuridae and the thylacine but absent in the numbat (Fig. 2a, Supplementary Fig. 4; marker 18).

To search for additional diagnostic markers to further clarify the phylogenetic position of the thylacine within Dasyuromorphia, we also conducted bioinformatics and experimental screens to test all possible tree topologies:

(1) A basal position of thylacine in Dasyuromorphia (2) A sister group relationship of thylacine and dasyurids (3) A close relationship of thylacine and numbat

To test hypothesis no. 1, we extracted 945 copies of WSINE1/1a elements plus 1,000 nucleotides of flanking sequences per side and locus from the Tasmanian devil genome that were absent in thylacine (detected per local blast89) for PCR amplification and sequencing in the numbat genomic DNA. To test hypotheses 2 and 3, we extracted 32,000 copies of WSINE1/1a elements plus their flanks from the thylacine genome. These loci were searched in batches as blastn queries against the Tasmanian devil genome (taxid:9305) at http://blast.ncbi.nlm.nih.gov (nucleotide blast). For hypothesis 2, loci with orthologous WSINE1/1a insertions in the Tasmanian devil (absent in koala) were included in a final dataset for experimental analysis in the numbat. For hypothesis 3, loci with an absence of WSINE1/1a elements in the Tasmanian devil were selected for experimental analysis in the numbat.

We randomly selected 75 of the most conserved (to enable reliable PCR) loci from each of the three datasets. These loci were manually aligned, and sequences from the short-tailed opossum (taxid:9265), the Tammar wallaby (taxid:9315), and the koala (taxid:38626) were aligned along with these where possible (Supplementary Table 23). Conserved, SINE-flanking PCR primers were designed and used for standard PCR amplification of genomic DNA from the numbat and other dasyuromorphs for all loci.

These three hypothesis-testing screens yielded 11 phylogenetically informative markers present in the numbat and dasyurids but absent in the thylacine (Supplementary Fig. 4; markers 30–40); three additional diagnostic WSINE1/1a elements present in thylacine and dasyurids but absent in the numbat and outgroup marsupials (Supplementary Fig. 4; markers 26–28); and one marker present in thylacine and numbat but absent in other dasyuromorphs (Supplementary Fig. 4; marker 29). Multidirectional KKSC analysis thus revealed significant support for the first hypothesis—a basal position of thylacine in Dasyuromorphia (11:4:1 markers; P < 0.005).

Retrophylogenetic tree topology and statistics

All markers and their presence/absence patterns are compiled in Supplementary Fig. 4. Supplementary Table 24a,b presents the data in a format accessible for a Dollo parsimony tree reconstruction (PAUP 4.0b1090; irrev.up character transformation; heuristic search with 1,000 random sequence additions; tree bisection and reconnection branch swapping) and a Bayesian tree reconstruction (MrBayes v3.2.591; standard discrete model; binary; ctype irreversible) to reveal the phylogenetic position of the thylacine shown in Fig. 2a. The branching points were evaluated using the KKSC insertion significance test for presence/absence markers16. Additional PCR primers are presented in Supplementary Table 22. FASTA alignments of phylogenetically informative markers are presented in Supplementary Table 23.

Reconstruction of demography with PSMC′

Retrophylogenomics

We inferred past effective population sizes separately for the thylacine and Tasmanian devil using the PSMCʹ method as implemented in the program MSMC92.

Read mapping and filtering

We downloaded Tasmanian devil reads from NCBI’s short-read archive, with accessions ERR789026, ERR789027, ERR789028, ERR789029, ERR789030, ERR789031 and ERR789032. The Tasmanian devil reads and our thylacine reads were mapped to the Tasmanian devil reference (Ensembl Devil_ref v7.0) with bwa mem77 using default parameters. Following mapping, we excluded scaffolds shorter than 1 Mb in length, and scaffolds that were putatively X-linked.

Determination of X-linked scaffolds

To remove X-linked scaffolds, we started with a list of 2,507 putatively X-linked scaffolds in the tammar wallaby reference (MEUG3.0, unpublished data), identified with homology to the human X (GRCh38) or containing tammar wallaby X-linked loci93. By mapping Tasmanian devil scaffolds to the tammar wallaby reference we could therefore identify X-linked scaffolds for the Tasmanian devil and thylacine. However, we first split the DEVIL7.0 scaffolds into contigs to allow more sensitive mapping, as attempts to map DEVIL7.0 scaffolds directly to the MEUG3.0 reference (A. J. Pask, personal communication) resulted in larger numbers of scaffolds remaining unmatched (25,348 of 35,975). We split individual scaffolds at the occurrence of one or more Ns. Resulting contigs were named for their scaffold of origin with an additional numerical suffix (for example, the 12th contig in scaffold GL834412.1 was named GL834412.1.12). The resulting whole-genome alignment was filtered to retain contigs that mapped to one of the MEUG3.0 X-linked scaffolds (no 0 × 4 flag), were primary mappings (no 0 × 100 or 0 × 800 flag), and had MAPQ ≥ 25. Contig names were truncated to obtain the originating scaffold names and collapsed for uniqueness with ‘sort −u’. This process identified 640 putatively X-linked devil scaffolds. The putatively X-linked Tasmanian devil scaffolds comprised 404 Mb of the total 3027.64 Mb DEVIL7.0 reference. Given the small size of marsupial sex chromosomes94, we expect that this list contains a high number of scaffolds falsely marked as X-linked. The exclusion of these scaffolds is thus considered to represent a conservative autosomal set.

After the removal of short and X-linked scaffolds we were left with 802 scaffolds, which had mean read depths for the Tasmanian devil and the thylacine of 59.85x and 41.79x, respectively. Alignments for each species were converted to the MSMC92 input format with samtools87 mpileup −E −A −q 20 −Q 20 −C 50, and bamCaller.py/generate_multihetsep.py (msmc-tools git~d99320d). Scaled population sizes, and the times to which they apply, were inferred using MSMC v1.0.0 for 50 expectation maximization iterations (−i 50).

Parameter estimation

MSMC produces an output that is scaled by the per-generation mutation rate, requiring estimates of the generation time and the per-generation mutation rate in order to rescale plots into real time. We assume a generation time of three years for both the Tasmanian devil and the thylacine. The Tasmanian devil reaches sexual maturity by two years, suggesting a mode of two years for parental age, with right skew. Thus an estimated generation time of three years is not unreasonable for these species95.

We first attempted to estimate the per-year mutation rate by dividing the number of observed substitutions between Tasmanian devil and thylacine by twice the time to their most recent common ancestor. Because of the large divergence time and low contiguity of the Tasmanian devil reference, only conserved regions of the genomes map to each other and as a consequence this strategy produced an unreasonably low estimate (1.74 × 10−12 mutations per base per year, or for a three-year generation time, 5.22 × 10−12 mutations per base per generation).

As an alternative, we took advantage of the highly conserved 2n = 14 karyotype of dasyurids to estimate a recombination rate, r, that may equally apply to the Tasmanian devil and the thylacine. In mammals, the number of recombinations per-generation is strongly predicted by the number of chromosome arms96. The Tasmanian devil has six metacentric autosomes, leading to the expectation that we should observe 12 recombinations per-generation. To obtain the recombination rate per base, per generation, we divide by the size of the autosome, where the autosome is assumed to be 97% of the genome97. This gives r = 9.147596 × 10−9 recombinations per base per generation.

MSMC provides estimates for the scaled recombination rate ρ = 4rN e and scaled mutation rate θ = 4μN e . MSMC uses Watterson’s estimator98 of θ, which we note has variance θ + θ 2. MSMC gave the scaled mutation rate as θ d = 6.69107 × 10−5 and θ t = 2.6563 × 10−5 for the Tasmanian devil and the thylacine, respectively. Provided the “–fixedRecombination” parameter is not used, MSMC estimates ρ during each Baum-Welch iteration, which is very accurate after 50 iterations92. MSMC gave the scaled recombination rate at the 50th iteration as ρ d = 2.34114 × 10−4 and ρ t = 5.78165 × 10−5 for the Tasmanian devil and the thylacine, respectively. For consistency with the scaled recombination rate estimate, we also infer the demography from the 50th iteration.

By multiplying the ratio of the scaled mutation rate to scaled recombination rate by the per-generation recombination rate implied by the karyotype we were able to estimate a per-generation mutation rate for scaling the MSMC output. The final mutation rate estimates are:

Tasmanian devil: μ d = rθ d / ρ d = 1.167813 × 10−9 mutations per base per generation.

Thylacine: μ t = rθ t / ρ t = 1.877287 × 10−9 mutations per base per generation.

Using the same strategy to estimate the human mutation rate from 1,000 Genomes sample HG00096, we obtained μ = 6.975339 × 10−8 mutations per base per generation. This value is faster than published rates, which range between 1 × 10−8 (pedigree estimates) and 2.5 × 10−8 (long-term divergence estimates) mutations per base per generation99. Therefore, we expect that the true thylacine and Tasmanian devil mutation rates may be slower than those we have estimated above. Mutation rate scaling is linear, and the effect of scaling with a slower mutation rate is to shift the demographic curve farther back in time and higher up on the population size axis. Consequently, the reconstructions for each species in Fig. 2b represent lower bounds for both timing and population size (that is, events may have occurred earlier and population sizes may have been larger).

Bootstrapping

A typical MSMC bootstrap is performed by first concatenating MSMC input files for all chromosomes, then sampling chunks with replacement until the sum of the sampled chunks is an equivalent genome size. A script for this can be found in the msmc-tools git repository. By default, this script samples 20 5 Mb chunks, placing them sequentially to make a 100 Mb pseudo-chromosome. 30 such pseudo-chromosomes are constructed with default settings. This bootstrap strategy accounts for: (1) variability in demographic inference along the chromosomes; and (2) variability due to a small rate of false recombinations, induced at chromosome boundaries and at chunk boundaries. For a chromosome level reference assembly, the impact of (2) should be minor. The scaffolds we have used are roughly exponentially distributed, ranging between 1–5 Mb (remembering that 1 Mb was our lower cutoff for a scaffold’s inclusion). Thus, the impact of (2) for our assembly is likely to be major, and not representative of falsely detecting recombinations from our data. We instead choose to sample scaffolds, with replacement, until the total original data size is reached. Because (1) should also include variability in erroneously detected recombinations, we believe this resampling strategy is more appropriate for data mapped to a scaffold level assembly. We did 100 bootstrap replicates for each of the Tasmanian devil and the thylacine. As new estimates for θ and ρ were obtained in each replicate, rescaling into real time was done separately for each replicate according to the procedure above.

Comparative genomic analyses

Orthologue identification and alignment

Thylacine and wild canid protein-coding DNA sequences (CDS) were extracted from the reference-based assembly using gffread100 and the Ensembl 84 annotations of the Tasmanian devil and dog genomes respectively. All CDS sequences for high-confidence one-to-one orthologues to Tasmanian devil genes from 21 vertebrates (inclusive of the devil) were downloaded from Ensembl 84 BioMart (Supplementary Table 8) and filtered for the best available CDS sequence per gene using custom scripts101. For each Ensembl gene ID, the longest complete CDS was selected. If no complete CDS was available, the longest sequence beginning with a start codon was accepted. Otherwise the longest available partial sequence was selected. Filtered sequences were grouped by orthology using ParaAT1.0102, then aligned and translated using MACSE version 1.01b103 with default parameters. In total, 11,254 groups of orthologous genes were retained for comparative genomic analyses (Supplementary Table 9). We used published phylogenies of all included species to construct a fixed tree topology104,105 used in all comparative genomic analyses (Supplementary Fig. 8).

Testing the frequency of homoplasious substitutions

To test the frequency of homoplasious substitutions across mammals, we collected all genes from our previous analyses which contained representative sequences for all therian mammals and the platypus as an outgroup. These genes were filtered to alignments with at least 50% average pairwise identity. Gene trees for the genes retained after filtering were used to construct a consensus tree with Mesquite v3.10 (build 765)106 using default parameters. We then counted the observed number of homoplasious substitutions for all pairwise comparisons between all nodes that do not share a most recent ancestral node in the fixed tree topology. We then compared it to the number of expected homoplasious substitutions under the JTT-F site model29, using a python script generously provided by Z. Zou and J. Zhang (University of Michigan). R values (observed/expected) were calculated for all pairs of species that did not share a most recent ancestral node and plotted against phylogenetic distance (Fig. 4 and Supplementary Table 10). Pairwise comparisons with expected values much lower than 1 resulted in artificially inflated R values, thus only comparisons with at least one expected homoplasious substitution were included in the plot.

Testing for individual homoplasious substitutions

To identify genes containing convergent and parallel substitutions between the thylacine and the canids, ancestral protein sequences were reconstructed using CodeML31 (Supplementary Table 25) using the fixed tree topology described above (Supplementary Fig. 8). Alignments with less than 50% average pairwise identity were excluded from analysis. The thylacine was compared to the ancestral sequence for all canids, as this represents the node in which adaptive substitutions are likely to have resulted in features shared by all canids. We used previously described definitions for parallel and convergent changes29. Briefly, for each amino acid position, the thylacine and the reconstructed ancestral sequence for the extant canids were compared to each other and to their most recent ancestors respectively. Positions that were identical between the thylacine and canid sequence, but different from both of their respective ancestors, were deemed to be homoplasious. Homoplasious substitutions in which the ancestors shared an identical residue were considered to be parallel, while residues that differed between the ancestors were considered to be convergent.

Testing for positive selection

The branch-site model implemented in CodeML (PAMLv4.7)31 was used to test for genes containing positively selected sites. The modified model A (Supplementary Table 26) was compared to the null (Supplementary Table 27) assuming the fixed tree topology described above (Supplementary Fig. 8). Four likelihood ratio tests were performed on nucleotide alignments of orthologous genes using different foreground branches; the thylacine, the branch leading to the subfamily Caninae, the Tasmanian devil and the branch leading to the family Bovidae (the latter two acting as controls). The total number of genes analysed for each foreground branch varies slightly (thylacine: 10,770; Caninae: 10,766; Tasmanian devil: 10,770; Bovidae: 9,544) due to the lack of representative sequences available for a given foreground branch in some alignments, and because PAML discards all alignment columns containing gaps. Additionally, alignments with average pairwise identity less than 50% were excluded from analyses. P values were obtained using a 50:50 mixture of a point mass 0 and a χ 2 distribution with 1 degree of freedom and were corrected for false discovery rate (FDR) using the Benjamini–Hochberg33 method with a cutoff of 0.05. Genes under positive selection in multiple lineages were found by comparing positively selected genes from each foreground branch (Supplementary Tables 12–13, 15–18 and 28–31). KEGG pathway enrichment among positively selected genes in thylacine and Caninae were each analysed using the Enrichr web server (http://amp.pharm.mssm.edu/Enrichr/). HGNC gene symbols for reference Tasmanian devil genes were downloaded from Ensembl 85 BioMart. Genes without HGNC gene symbols were excluded (highlighted in red, Supplementary Tables 12 and 13). The resulting lists of HGNC symbols for the thylacine and canids were then separately analysed for enrichment of KEGG 2016 pathways, the results of which are contained in Supplementary Fig. 9 and Supplementary Tables 19 and 20.

Life Sciences Reporting Summary

Further information on experimental design is available in the Life Sciences Reporting Summary.

Code availability

Custom PERL scripts were used for data parsing and to identify convergent and parallel amino acid substitutions. Code is available at https://github.com/charlesfeigin/thylacine_genome_project_scripts. Custom python scripts for implementing the JTT-F site model are available at https://github.com/ztzou/convCal.

Data availability

Genome sequence reads have been deposited in the Sequence Read Archive as BioProject PRJNA354646, BioSample SAMN06049672.