Genome of the Tasmanian tiger provides insights into the evolution and demography of an extinct marsupial carnivore

The Tasmanian tiger or thylacine (Thylacinus cynocephalus) was the largest carnivorous Australian marsupial to survive into the modern era. Despite last sharing a common ancestor with the eutherian canids ~160 million years ago, their phenotypic resemblance is considered the most striking example of convergent evolution in mammals. The last known thylacine died in captivity in 1936 and many aspects of the evolutionary history of this unique marsupial apex predator remain unknown. Here we have sequenced the genome of a preserved thylacine pouch young specimen to clarify the phylogenetic position of the thylacine within the carnivorous marsupials, reconstruct its historical demography and examine the genetic basis of its convergence with canids. Retroposon insertion patterns placed the thylacine as the basal lineage in Dasyuromorphia and suggest incomplete lineage sorting in early dasyuromorphs. Demographic analysis indicated a long-term decline in genetic diversity starting well before the arrival of humans in Australia. In spite of their extraordinary phenotypic convergence, comparative genomic analyses demonstrated that amino acid homoplasies between the thylacine and canids are largely consistent with neutral evolution. Furthermore, the genes and pathways targeted by positive selection differ markedly between these species. Together, these findings support models of adaptive convergence driven primarily by cis-regulatory evolution. The Tasmanian tiger is an extinct carnivorous marsupial. By sequencing the genome of a preserved specimen the authors show long-term population decline and reveal the genetic basis of the phenotypic convergence between Tasmanian tigers and canids.

T he Tasmanian tiger or thylacine (Thylacinus cynocephalus) was a large, carnivorous marsupial and the only species within the family Thylacinidae to survive into the modern era. Historically the thylacine was broadly distributed across Australia before becoming extinct on the mainland around 3,000 years ago 1 . A Tasmanian population became isolated by rising sea levels 2 approximately 14,000 years ago and persisted until the early twentieth century. European settlers deemed the thylacine a threat to the Tasmanian sheep industry and the government aggressively targeted it for eradication by offering a £1.00 bounty for each animal killed 1 . Consequently, the remaining population was rapidly exterminated and the last known thylacine died at the Hobart Zoo in 1936. The species was eventually declared extinct in 1982.
Apart from having a pouch in which its young developed and a conspicuous striped pattern on its hindquarters from which it derived its common name, the marsupial thylacine was phenotypically almost indistinguishable from a eutherian canid 1 (Fig. 1a,b). Because of their striking similarities and ancient divergence time (~160 million years ago (Ma)) 3 , the thylacine and canids are a widely recognized example of convergent evolution. Many cases of convergent evolution are thought to arise in response to similar selective pressures driving similar adaptations in different species 4 . The phenotypic resemblance between the thylacine and canids is strongly associated with a primarily carnivorous feeding ecology 5 and, correspondingly, both lineages exhibit extraordinary similarity in cranial shape 5,6 . However, the thylacine's extinction and the resulting paucity of molecular data has thus far prevented analyses that could clarify the genetic basis of this adaptive phenotypic convergence. Furthermore, our current understanding of several key aspects of the thylacine's evolutionary history, including its placement within the carnivorous marsupials and its genetic diversity before extinction, are built on mitochondrial data and only a small number of nuclear loci [7][8][9] . Therefore, we sequenced the thylacine genome to improve our understanding of this enigmatic marsupial predator.

Results and discussion
Thylacine genome sequencing. We extracted DNA from the soft tissue of a 108-year-old, alcohol-preserved thylacine pouch young specimen (Fig. 1c) from Museums Victoria, Australia. DNA fragments between 300-600 bp were isolated and sequenced on Illumina platforms (Supplementary Fig. 1 and Supplementary  Table 1), generating ~188 gigabases (Gb) of sequence data that were subsequently filtered for quality and length. Contaminating sequences were depleted by mapping to a database containing microbial and fungal genomes, leaving ~151 Gb (Supplementary Table 1). We assessed the quality of the data by mapping reads to the genome of the Tasmanian devil (Sarcophilus harrisii) 10 , the most closely related species available. Approximately 89.3% of the retained reads mapped to the complete unmasked Tasmanian devil assembly. Additionally, ~61.6% of reads mapped to the repeatmasked reference assembly, covering ~84.4% of unmasked bases to an average depth of ~42.8× (Supplementary Table 2). For comparison, we mapped the thylacine data to two other marsupial genomes, the tammar wallaby (Macropus eugenii) 11

and the gray short-tailed
The Tasmanian tiger or thylacine (Thylacinus cynocephalus) was the largest carnivorous Australian marsupial to survive into the modern era. Despite last sharing a common ancestor with the eutherian canids ~160 million years ago, their phenotypic resemblance is considered the most striking example of convergent evolution in mammals. The last known thylacine died in captivity in 1936 and many aspects of the evolutionary history of this unique marsupial apex predator remain unknown. Here we have sequenced the genome of a preserved thylacine pouch young specimen to clarify the phylogenetic position of the thylacine within the carnivorous marsupials, reconstruct its historical demography and examine the genetic basis of its convergence with canids. Retroposon insertion patterns placed the thylacine as the basal lineage in Dasyuromorphia and suggest incomplete lineage sorting in early dasyuromorphs. Demographic analysis indicated a long-term decline in genetic diversity starting well before the arrival of humans in Australia. In spite of their extraordinary phenotypic convergence, comparative genomic analyses demonstrated that amino acid homoplasies between the thylacine and canids are largely consistent with neutral evolution. Furthermore, the genes and pathways targeted by positive selection differ markedly between these species. Together, these findings support models of adaptive convergence driven primarily by cis-regulatory evolution.

Articles
Nature ecology & evolutioN opossum (Monodelphis domestica) 12 . As expected, the percentage of reads that mapped decreased with increased phylogenetic distance between the reference genome and thylacine ( Supplementary Fig. 2). The thylacine sequence data was then used to assemble contigs de novo. The resulting assembly was comparable to other draft marsupial genomes [10][11][12] , with an N50 contig length of ~3.2 kilobases (Kb) and a total assembly size of ~3. 16 Gb, suggesting the thylacine genome size was similar to previously sequenced marsupials. The G + C content of the de novo contigs (~36%), was similar to the thylacine's close relative the Tasmanian devil (36.4%), and considerably lower than human (46.1%) and mouse (51.2%) (Supplementary Table 3). MapDamage 13 analysis of reads mapped to the thylacine de novo contigs failed to detect cytosine deamination, a standard characteristic of historical and ancient DNA ( Supplementary Fig. 3). This is likely to be due to the relatively young age of the sample and its method of preservation.
Phylogenetic placement of the thylacine. The placement of the thylacine within the order of carnivorous marsupials, Dasyuromorphia, has been a subject of much contention for over a century 8,9,14 . Previous analyses based on mitochondrial sequences supported the thylacine as the basal lineage in Dasyuromorphia, followed by the numbat as the sister lineage to Dasyuridae (the family comprising the Tasmanian devil and its close relatives) 8 . Unlike nucleotide substitutions, retroposon insertions are nearly always permanent and independent, and parallel insertions are exceptionally rare. This makes retroposon insertion patterns virtually homoplasy free, and therefore reliable indicators of phylogenetic history. Here, we screened the genome sequences of the Tasmanian devil and thylacine for 25 diagnostic phylogenetic retroposon markers (SINEs; short interspersed elements) previously identified in dasyuromorph species 15 . In addition, we PCR amplified and sequenced an additional 225 potentially informative loci in the numbat, resulting in a total of 250 loci. Of these, 11 supported a tree in which the thylacine represents the basal lineage in Dasyuromorphia followed by the numbat as the direct sister group to Dasyuridae (Fig. 2a). However, our multidirectional screening (see Methods) also revealed four markers supporting a potential sister group relationship of Dasyuridae and thylacine, and one marker uniting numbat and thylacine. Multidirectional KKSC 16 analysis revealed significant support for the basal position of thylacine in Dasyuromorphia (11:4:1 markers; P < 0.005). This represents the first time that the genome-wide presence/absence patterns of retroposons have been used to clarify phylogenetic relationships of an extinct species. The presence of five conflicting signals in a multidirectional bioinformatics and experimental screening suggests that ancestral incomplete lineage sorting of retroposed elements, frequently detected in many other vertebrate groups 17,18 including ancestral marsupial lineages 19 , also occurred in Dasyuromorphia.
Demographic history of the thylacine. We assessed the recent demographic history of the thylacine and compared it to the  Phylogenetic tree of Dasyuromorphia based on presence/absence insertion patterns of SINE elements. Green and red circles represent diagnostic insertions of WALLSI1A and WSINE1/1a, respectively, leading to the significantly supported tree topology (solid branches). Black circles represent cases of incomplete lineage sorting (ILS), indicating non-significantly supported alternative tree topologies (dashed branches). b, Demographic history of the thylacine (yellow) and Tasmanian devil (black), both scaled using a generation time of 3 years and mutation rates calibrated from the conserved karyotypes (see Methods). One hundred bootstrap replicates for each species are shown as thinner lines of the same colour. Vertical grey lines indicate the approximate timing of the arrival of ancestral Australian Aboriginals in Australia (~48,800 years ago, based on synthesis of archaeological evidence 107 ) and the isolation of Tasmania from the Australian mainland (~14,000 years ago, based on shoreline reconstructions 2 ). N e denotes effective population size.

Nature ecology & evolutioN
Tasmanian devil, another marsupial carnivore that underwent isolation on Tasmania. The Tasmanian devil has experienced a recent catastrophic population decline due to an extreme lack of genetic diversity, causing widespread susceptibility to devil facial tumour disease 10 . Previous studies of mitochondrial DNA suggested that thylacine populations in Tasmania were even less genetically diverse at the time of their extinction than modern Tasmanian devils 7 . We performed pairwise sequentially Markovian coalescent (PSMCʹ ) analyses, which revealed that the effective size of the thylacine population largely mirrored that of the Tasmanian devil, with both species undergoing an apparent steep decline in diversity beginning around ~70-120 ka (thousand years ago) (Fig. 2b). The population decline appears to have begun before the human colonization of Australia (> ~50 ka) 20 and overlaps with climate changes associated with the beginning of the penultimate glacial cycle 21 . Mutation rates for the thylacine and Tasmanian devil were conservatively estimated (see Methods), and we expect that the true thylacine and Tasmanian devil mutation rates may be slower. Therefore, reconstructions for each species (Fig. 2b) are likely to represent lower bounds for both timing and population sizes. However, PSMCʹ analyses can also be impacted by marked population structure, and simulations of a constant sized metapopulation with a large number of demes and low migration 22 generated similar trends to those in Fig. 2b. Such population structuring is plausible for a large, geographically dispersed Australian carnivore.
Quantifying phenotypic convergence. We quantified the degree of cranial shape convergence between the thylacine and canids using 3D landmark-based geometric morphometrics on a broad range of modern and extinct mammalian taxa within a phylogenetic framework (Supplementary Fig. 5 and Supplementary Tables 4 and 5). Despite a substantial clustering of metatherians and eutherians in cranial morphospace ( Fig. 3 and Supplementary Figs. 6 and 7), diet was a significant predictor of cranial shape (Procrustes ANOVA, F = 8.899, P < 0.0001), even after accounting for phylogeny (phylogenetic ANOVA, F = 2.199, P = 0.008; Supplementary Table 6).  Supplementary Fig. 10 projected into shape space using squared-change parsimony. More than half of the total variance was contained in the first two components: PC1 describes shortening or lengthening of the face and jaws and relative height of the sagittal crest, while PC2 describes variation in skull width and snout shape (tubular versus blunt). The thylacine (blue star) displays striking cranial convergence with eutherian canids, falling closer in shape space to the red fox (Vulpes vulpes; shown) and grey wolf (Canis lupus) than to its own closest relatives. Using the pairwise distance measure C 1 23 on the full phenotypic dataset (PC1-31; 99% of the total shape variance), we estimate that approximately 34% of the ancestral phenotypic space between the thylacine and red fox has been closed by convergent evolution and 40% has been closed between the thylacine and grey wolf. These values are significantly greater than would be expected by chance (convratsig C 1 , P < 0.05, n = 113). In contrast, the thylacine did not exhibit cranial convergence with other carnivorous marsupials or eutherian felids (Supplementary Table 7). b, Locations of 30 homologous landmarks used to compare cranial shapes shown in black (dorsal and lateral) and grey (ventral; expanded in Supplementary Fig. 5).

Nature ecology & evolutioN
Using distance-based phenotypic measures on 99% of the total shape variation (PC1-31), we estimated that over one-third of the ancestral phylomorphospace between the thylacine and canids (such as Canis lupus and Vulpes vulpes) has been closed by convergent evolution ( Fig. 3 and Supplementary Table 7). This level of convergence is similar to that calculated using the same metric for Caribbean Anolis lizards, a classic example of convergent evolution 23 . However, in contrast to the closely related anoles, metatherians and eutherians have been separated for over 160 million years of independent evolution 3 . Thus, the thylacine-canid comparison represents an exceptional model of convergent evolution in distantly related species.
Genomic convergence. Selection is often suggested to repeatedly target a limited set of conserved developmental genes, constituting a genetic toolkit 24 . Under this assumption, molecular convergence in orthologous loci is suggested be an important driver of phenotypic convergence 25 . Recent genome-wide studies have attempted to identify homoplasious amino acids in protein-coding genes, notably among several marine mammal lineages and among echolocating mammals 26,27 . Both studies found instances of homoplasy in genes with plausible links to their convergent phenotypes 26,27 . However, similar amounts of homoplasy were also detected between non-convergent species 26 . Further analyses showed that homoplasious substitutions between marine mammals were largely consistent with relaxed purifying selection 28 . However, it could be argued that overall phenotypic convergence among marine and among echolocating mammals is not extensive. Therefore, it is unclear whether the scarcity of adaptive homoplasy is due to evolutionary constraints on adaptation via protein-coding changes, or simply limited phenotypic convergence in the species being examined. Furthermore, a recent study showed that the frequency of convergent and parallel substitutions in protein-coding genes declines with increasing phylogenetic distance 29 . It was suggested that accumulation of changes in the genomic background can modify epistatic interactions, limiting the opportunities for adaptive homoplasy 29,30 . While that study examined a broad sampling of mammalian diversity, the species included also showed little phenotypic convergence. In contrast, the thylacine and canids possess highly convergent phenotypes while also sharing one of the deepest evolutionary splits in mammals. Therefore, we interrogated the genomes of the thylacine and canids to determine if their extensive phenotypic similarities were reflected by adaptive amino acid homoplasy.
We reconstructed thylacine genes by mapping the thylacine sequence data to the Tasmanian devil genome 10 and generating a reference-based assembly. The Tasmanian devil genome annotations were then used to extract thylacine genes. We used a similar approach to assemble genes from publicly available sequence data for five wild canid species: wolf, golden jackal, coyote, red fox and arctic fox. The resulting thylacine and wild canid genes were grouped with high-confidence one-to-one orthologues from 21 other vertebrates including 14 eutherian mammals and 3 marsupials, with platypus and birds as outgroups (Supplementary Table 8). We assumed a species tree topology based on published phylogenies of the species represented in our orthologue dataset ( Supplementary Fig. 8). After filtering for high-quality sequences, we retained 11,254 groups of orthologous genes (Supplementary Table 9). Additionally, we performed ancestral state reconstructions using CodeML 31 that were used for subsequent analyses.
Genome-wide prevalence of amino acid homoplasy. We investigated the genome-wide prevalence of amino acid homoplasy in the thylacine and ancestors of canids, comparing the observed number of homoplasious amino acids with the number expected under the neutral JTT-F site model 29 . This relative frequency was further compared to that calculated for 691 additional pairwise comparisons of nodes in our phylogeny ( Fig. 4 and Supplementary Table 10). Our analysis recovered the reported negative relationship between phylogenetic distance and the frequency of homoplasious substitutions. Notably, the prevalence of homoplasy between the thylacine and canids did not defy expectations based on neutral evolution and was comparable to pairwise comparisons between non-convergent species with similar phylogenetic distances ( Fig. 4 and Supplementary Table 10).
Amino acid homoplasy and positive selection. If the complement of genes in the genetic toolkit is indeed limited, then phenotypic convergence may not be reflected by an enrichment of homoplasy. Rather, molecular changes in a small number of key developmental genes may be targeted by selection to drive convergent adaptations. To test this, we first compared individual amino acid positions between the thylacine and the ancestor of extant canids. We then identified positively selected protein-coding genes using the branchsite model 32 implemented in CodeML (PAMLv4.7) 31 . We performed likelihood ratio tests for genes on the branches leading to the thylacine and the subfamily Caninae respectively, and applied a false discovery rate (FDR) of 0.05 (Benjamini-Hochberg) 33 . The families Thylacinidae and Dasyuridae diverged around 26 Ma and early thylacine fossils in the Riversleigh deposits indicate the thylacine morphology was fully established by 20 Ma 34 . Thus, selection driving the thylacine's phenotypic convergence with canids was likely to have been confined to a 6 million year period, and these distant signatures of positive selection might be masked by subsequent sequence evolution. By allowing variable selective pressures among sites, the branch-side model has a greater power to detect historical episodes of selection over other codon-based methods 35 . We defined adaptively convergent genes as those that showed signatures of positive selection in both the thylacine and canids, and contained homoplasious amino acid substitutions.

Nature ecology & evolutioN
In total, we identified 16 convergent and 67 parallel changes between the thylacine and canid ancestor, contained within 81 orthologous genes (Supplementary Table 11). Among these genes were BRPF1, which maintains the anterior-posterior axis of craniofacial bones in fish 36 ; DEGS1, a predicted target of RUNX2 37 , a known driver of craniofacial diversity in carnivorans 38 ; and LTBP1 and LRIG1, both of which are associated with snout length in mice 39,40 . After FDR correction, we identified 259 genes under positive selection in the thylacine and 185 in the branch leading to Caninae (Supplementary Tables 12 and 13), which were then intersected with the list of genes containing homoplasious amino acid substitutions. Interestingly, none of the 81 homoplasious genes were found to be adaptively convergent in the thylacine and canids. Five homoplasious genes were under positive selection in only the thylacine, including BRPF1 and DEGS1. However, we did not find evidence of adaptation in their canid orthologues (Supplementary  Table 14). Thus, the observed homoplasies do not appear to be driven by convergent selective pressures. We identified a further six non-homoplasious genes that showed signatures of positive selection in both the thylacine and Caninae (Supplementary Table 15). For comparison, we performed additional branch-site likelihood ratio tests on two non-convergent lineages with comparable divergence times: the Tasmanian devil (a dasyurid marsupial) and the eutherian Bovidae (a family consisting exclusively of ruminant herbivores). Despite the slightly smaller number of positively selected genes identified in the Tasmanian devil and Bovidae (209 and 40 genes, respectively; Supplementary Tables 16 and 17), this comparison returned five overlapping genes, nearly the same number observed between thylacine and canids (Supplementary Table 18). From this we concluded that positive selection has not targeted orthologous genes more frequently in the thylacine and canids than in non-convergent species with similar divergence times. This finding may suggest that protein-coding genes are not a major driver of convergent evolution in therian mammals. Alternatively, selection in these lineages may have targeted different genes involved in similar biological processes. Therefore, we performed a KEGG pathway analysis 41 of thylacine and canid positively selected genes using Enrichr 42 . We identified one pathway, inflammatory mediator regulation of TRP channels (hsa04750), among the top terms in both canids and the thylacine (Supplementary Tables 19 and 20). Among the genes identified in this pathway were phospholipase A2 orthologues (PLA2G6 and PLA2G4E in the thylacine and PLA2G4F in canids; Supplementary Fig. 9). This pathway is associated with thermal hypersensitivity and may be involved in persistent pain responses 41 , and does not appear to be associated with convergent adaptations between the thylacine and canids. Thus, it appears that selection on genes in conserved pathways is very limited and unlikely to be relevant to convergence.
In summary, our data demonstrate that observed amino acid homoplasies between thylacine and canid orthologous protein-coding genes are largely consistent with neutral evolutionary processes. Therefore, such molecular changes are unlikely to explain the convergent phenotypes between these species. Furthermore, our data demonstrate that there is surprisingly little concordance between positive selection on thylacine and canid protein-coding orthologues regardless of the presence of homoplasy. Given the extraordinary adaptive similarities between these distantly related taxa, we suggest that positive selection has targeted other genomic regions such as cis-regulatory elements (CREs) to drive their phenotypic convergence. CREs such as promoters and enhancers regulate spatial and temporal gene expression and have a modular structure. Consequently, many CREs may have greater evolutionary flexibility than their target genes and these elements have been suggested as preferred targets of natural selection 43 . Cis-regulatory evolution may be particularly important for driving morphological adaptations 43 , which involve changes in body patterning. We suggest that future genome-wide studies on the molecular basis of convergent phenotypic evolution should explore the contribution of cis-regulatory changes to this phenomenon.

Methods
Geometric morphometric analyses. Taxon sampling. Taxon sampling was largely based on a previous study 6 , which itself was an expansion of the 3D cranial dataset assembled by Wroe & Milne 5 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 species 44 . Among these are the well-known large-bodied carnivores such as lions, wolves and hyenas, but also omnivores (some bears), insectivores (aardwolf and bateared 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 hyenasized and were the most successful and abundant predators of the early Tertiary period, an ecological niche later dominated by carnivorans 45 . Among metatherians, carnivory (here defined as the killing and comminuting of vertebrate prey, sensu 46 ) 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 study 5 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 study 6 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 dataset 6 by an additional 35 species, including members of Canis and Vulpes (the eutherian genera to which the thylacine is most commonly compared 46-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 noncarnivorous 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 eutherians 54 . 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 dataset 6 . 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 analysis 6 , 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 crest 6 . 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 Nature ecology & evolutioN 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 fit 55 . 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 modularity 56 , 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 approach 57 . 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.5 58 . 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 published 59 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 nodes 61 . 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 ANOVA 62 . 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 analysis 63 . Due to the non-independence of related species 64 , 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 contrasts 65 , 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 time 66 . Although D-PGLS is currently limited to BM, comparisons of phenotypic patterns derived from BM and other processes such as the Ornstein-Uhlenbeck model 67 may be more appropriate in the future 65 . 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 other 23 . Unlike other convergence methods like SURFACE 68 , 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 resembles 46,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 evolution 23 . 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.16 69 . 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.06d 55 , except for K mult and ANOVA, which were estimated in the geomorph package of R 70 , and convergence measures C 1 -C 3 calculated in the R package convevol 23 . 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 ethanolpreserved 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 1994 71 . 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 previously 72 . 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 expected 73 , 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.32 74 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 Nature ecology & evolutioN 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 pipeline 75 . Briefly, residual adaptor sequences were trimmed and reads shorter than 25 nucleotides were removed, and overlapping paired reads were collapsed using AdapterRemoval v.2.0 76 . Mapping of collapsed and paired reads was performed using the mem algorithm with default parameters in bwa v.0.7.13 77 .
Reads with mapping quality below 25 were removed. Finally, mapDamage 2.0 13 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, SOAPdenovo2 82 (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 vcf2fq 87 . 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 markers 15 . 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 test 16 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 blast 89 ) 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.0b10 90 ; 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.5 91 ; 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 markers 16 . 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 MSMC 92 .
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 mem 77 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 loci 93 . 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 chromosomes 94 , 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 MSMC 92 input format with samtools 87 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).

Nature ecology & evolutioN
Parameter estimation. MSMC produces an output that is scaled by the pergeneration 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 species 95 .
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 threeyear 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 arms 96 . 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 genome 97 . 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 estimator 98 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 iterations 92 . 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 generation 99 . 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 gffread 100 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 scripts 101 . 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.0 102 , then aligned and translated using MACSE version 1.01b 103 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 topology 104,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 model 29 , 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 CodeML 31 (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 changes 29 . 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-Hochberg 33 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

Replication
Describe whether the experimental findings were reliably reproduced.
Boostrapping procedure of PSMC' analysis of the genome is detailed in the methods.

Randomization
Describe how samples/organisms/participants were allocated into experimental groups.
Not applicable.

Blinding
Describe whether the investigators were blinded to group allocation during data collection and/or analysis.

Not applicable.
Note: all studies involving animals and/or human research participants must disclose whether blinding and randomization were used.

Statistical parameters
For all figures and tables that use statistical methods, confirm that the following items are present in relevant figure legends (or in the Methods section if additional space is needed).

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement (animals, litters, cultures, etc.) A description of how samples were collected, noting whether measurements were taken from distinct samples or whether the same sample was measured repeatedly A statement indicating how many times each experiment was replicated The statistical test(s) used and whether they are one-or two-sided (note: only common tests should be described solely by name; more complex techniques should be described in the Methods section) A description of any assumptions or corrections, such as an adjustment for multiple comparisons The test results (e.g. P values) given as exact values whenever possible and with confidence intervals noted A clear description of statistics including central tendency (e.g. median, mean) and variation (e.g. standard deviation, interquartile range)

Clearly defined error bars
See the web collection on statistics for biologists for further resources and guidance.