Contemporary Demographic Reconstruction Methods Are Robust to Genome Assembly Quality: A Case Study in Tasmanian Devils

Improved Reference Genome Assembly

We produced an improved genome assembly for the Tasmanian devil using the Chicago library method (Dovetail Genomics, Santa Cruz, CA). Hereafter, we refer to the genome of Murchison et al. (2012) as the 2012-Assembly, and the genome presented herein as the 2019-Assembly. The Dovetail pipeline produced a vastly improved genome assembly according to several metrics (table 1). Whereas assembly sizes are approximately equal (2012-Assembly: 3,174.71 Mb, 2019-Assembly: 3,177.42 Mb) the number of scaffolds was reduced from 35,975 to 10,010, and the L90 was reduced from 3,263 to 422 scaffolds. The N90 improved from 0.101 Mb in the 2012-Assembly to 2.164 Mb in the 2019-Assembly (table 1). A summary of the assembly contiguities may be found in supplementary figure S1, Supplementary Material online. BUSCO completeness using the mammalia obd9 database for each assembly is shown in supplementary table S1, Supplementary Material online. This new reference genome has been made publicly available and has been deposited at RefSeq under BioProject PRJNA504904. We compare these two assemblies to other assemblies for both model and nonmodel organisms in table 2.

Whole-Genome Resequencing

We resequenced the genomes of 12 devils trapped in northwestern Tasmania (fig. 1). The 12 resequenced genomes were subsequently aligned to the two alternative reference genomes. Across the 12 resequenced samples, using GATK, we identified 2,822,764 raw SNPs when aligning to the 2012-Assembly (Murchison et al. 2012) and 2,790,031 raw SNPs when aligning to the 2019-Assembly. However, following filtering with vcftools, we recovered 1,468,034 SNPs using the 2012-Assembly, as compared with 1,468,044 SNPs using the 2019-Assembly. The per-base mutation rate per-generation was estimated at 1.394808×10^–08. For a description of how this estimate was obtained, see Materials and Methods (see Calculation of the Mutation Rate). Mean coverage across the 12 resequenced samples aligned to the 2012-Assembly was 26.27 ± 2.68. In contrast, mean coverage for the 2019-Assembly was 28.47 ± 2.90. A more detailed summary of per-sample read mapping and depth of coverage may be found in supplementary table S2, Supplementary Material online.

Demographic Inference

Genome Assembly Quality Does Not Impact Inference of Population Size History

Surprisingly, we found each method to be broadly robust to genome assembly quality (figs. 2 and 3). Importantly, this is true for methods based on the SMC, which based on prior expectations could be sensitive to variation in scaffold length distributions. Further, bootstrap estimation variance did not increase with increasing genome fragmentation for any method as assessed using empirical data (fig. 2).

Nonparametric bootstrapping is used by all methods discussed herein to assess parameter estimate uncertainty. This assumes that the method is estimating without error, given a sample (in this case a collection of sequences or SFS). We argue that future efforts would benefit from quantifying and reporting this estimation uncertainty. One possible remedy could be Bayesian sampling of a posterior distribution of local TMRCAs, leading to a subsequent posterior distribution of estimates of N_e at each time step. Alternatively, using methods derived from the SMC, the distribution of estimated local TMRCAs and associated population sizes could be instead represented as an empirical distribution of N_e through time for a given sample. We also are not aware of any method that implements an approach to test model adequacy. Parametric bootstrapping, through simulation under the fitted model and resampling of estimates could prove useful to accomplish this.

One methodological consideration that could impact our empirical results is the manner in which each method handles gaps. Because Dovetail Genomics inserts a gap of 100 N’s for every join it makes, our improved assembly (2019-Assembly) has many more gaps than does the 2012-Assembly (table 1), with each gap corresponding to a join between two of the 2012-Assembly’s contigs. Without a priori specification of gaps, this could lead SMC methods to interpret that contigs comprising the input scaffolds are closer than they are. Further, these gaps may be interpreted as runs of homozygosity, leading to the estimated population sizes being artificially depressed. Interestingly, we see little evidence that the greater number of gaps in the 2019-Assembly had an impact on our results, as inferences obtained using each assembly are extremely similar. Further, inferences are qualitatively similar between those obtained using simulated and empirical data. This instills further confidence that our results are robust to the handling of gaps, as simulated data were generated assuming a demography approximately equal to that inferred using empirical data.

Application of each method to simulated genomes that were progressively fragmented to increasing degrees only further supported these findings. Accuracy of parameter estimates by methods leveraging the SMC worsened for recent (∼333 gbp) and more ancient (>30 kgbp) time periods beginning at ∼5,000 fragments. Most methods reliant on the SMC exhibit a high degree of precision between 50 kgbp and 333 gbp (but see MSMC), whereas precision subsequently decreased for more recent estimates.

Use of SMC- or SFS-Based Methods Recovers Alternative Population-Size Histories

Whereas inferences were similar among genomes within each method, incongruent results were observed when comparing methods using either the SMC or SFS using both empirical and simulated data. Specifically, inferences made by the Stairway Plot were largely discordant with those made by PSMC, MSMC, and SMC++. In contrast, SMC++, which leverages both the SMC and the SFS, recovered results intermediate to the other SMC methods and the Stairway Plot toward the present, likely due to the conflicting signals picked up by these two approaches. These results were consistent with previous work showing methods dependent on the SMC suffer from difficulties in estimating recent population history, but excel in the distant past (Li and Durbin 2011; Liu and Fu 2015; Terhorst et al. 2017). In contrast, methods exclusively dependent on the SFS excel at inferring more recent history, but perform less well when estimating more ancient population size history (Liu and Fu 2015; Terhorst and Song 2015; Lapierre et al. 2017; Baharian and Gravel 2018; Rosen et al. 2018).

An important caveat regarding the use of the SFS is its sensitivity to sample size. Specifically, numerous studies have demonstrated that the accuracy and precision of parameter estimates obtained using the SFS positively correlates with sample size (Adams and Hudson 2004; Keinan and Clark 2012; Robinson et al. 2014; Liu and Fu 2015; Terhorst and Song 2015). Likewise, it has been shown that for any sample size, parameter estimates are most error-prone for ancient demographic events. Importantly, this tendency is exacerbated by small sample size. However, increasing sample size has been shown to not only increase parameter estimate accuracy and precision for ancient demographic events but also for recent ones (Robinson et al. 2014; Liu and Fu 2015). In this context of our study, the discordance between inferences obtained using the Stairway Plot, which exclusively relies on the SFS, and the other methods may be explained by our low sample size. However, limited sample sizes for nonmodel organism whole-genome resequencing projects may be common (table 2).

Geologic, Climatic, and Anthropogenic Drivers of the Demographic History of the Tasmanian Devil

We found evidence broadly supporting several hypothesized biogeographic, anthropogenic, and climatic drivers of the Tasmanian devil’s demographic history. First among these was a population decline during the penultimate deglaciation in which Tasmania was separated from mainland Australia (135–43 kybp Lambeck and Chappell 2001). Only SMC++ and PSMC detected this decline, with a larger reduction in N_e by ∼50,000 (fig. 2) detected by SMC++. Given the persistence of devils on the mainland following the re-emergence of the Bass Strait land bridge ∼43 kybp following the penultimate deglaciation, admixture among previously isolated mainland and island lineages may have masked some of the impacts of this event on genetic diversity (Chikhi et al. 2010).

Although ethnographic accounts indicate that Aboriginals may not have eaten carnivorous marsupials such as the Tasmanian devil (Owen and Pemberton 2005), it is possible that their presence may have negatively influenced devil population density and size via persecution and competition. Consequently, we hypothesize that the penultimate deglaciation and the colonization by Aboriginals each may have independently contributed to the observed population decline, or that it may have been their shared influence that initiated and maintained this decline. Although this hypothesis is compatible with our results, we cannot at this time tease apart the relative contributions of the anthropogenic and climatic/geologic factors.

During the LGM, sea surface temperatures were ∼4 °C cooler than at present, and climate in Tasmania was cooler and drier than prior to the LGM, with landcover being dominated by grassland and open shrub (Barrows and Juggins 2005). These climatic changes potentially restricted distributions of prey of the Tasmanian devil (Colhoun and Shimeld 2012; Brüniche-Olsen et al. 2014), and consequent submersion of the Bass Strait ∼14 kybp (Lambeck and Chappell 2001) likely restricted the distribution of devils. Indeed, we find a steepening of the decline of N_e across all methods consistent with these events. Following this most recent isolation of Tasmania from the continent and extinction of devils on mainland Australia ∼3.5 kybp (White et al. 2018), alleviation of bottlenecks in island devils by emigration from mainland populations was no longer possible. Furthermore, the TMRCA of Tasmanian and southern mainland populations of devils has been estimated at ∼26,400 ya using whole mtDNA genomes (Brüniche-Olsen et al. 2018). Notably, this timing is coincident with the LGM and corresponds to the large-scale decline in N_e estimated by most methods (fig. 2). Jointly, these factors lend additional support to the interpretation that the separation of Tasmania from mainland Australia during this time led to population subdivision and subsequent population decline. The severity of these declines is further explained by the documented coincident change in climate, and landcover. An alternative explanation for these declines is that for more ancient time periods (i.e., prior to the LGM), the population history being estimated could be for that of mainland populations. Following the LGM, inferred population size histories may be for Tasmanian population.

Extreme ENSO events 5–3 kybp led to periods of extended drought (Donders et al. 2008), further affecting devil population sizes via prey population declines (Brüniche-Olsen et al. 2014, 2018). Although the onset of ENSO events, increased populations of Aboriginals and the introduction of dingos have been implicated in the mainland extinction of devils, dingos were never introduced on Tasmania, and Tasmanian Aboriginal populations remained at lower densities than on the mainland (Brüniche-Olsen et al. 2018). Thus, only the impacts of the ENSO events are shared between Tasmanian and mainland populations of S. harrisii (Brüniche-Olsen et al. 2018), mirroring the drivers of mainland extinction of the thylacine (White et al. 2018). However, our results show that the ENSO events occurred near the end of the population decline (fig. 2B), suggesting that they are unlikely to have been the sole driver of the inferred decline.

As recently as 1830, European colonists made concerted efforts to reduce, if not exterminate, devils from rural areas, considering them as pests that threatened livestock (Owen and Pemberton 2005). Continued efforts to cull populations of devils resulted in the further reduction of population sizes. Not surprisingly, only the Stairway Plot detected the reduction in N_e around this time (fig. 2) due to the increased sensitivity of the SFS to recent demographic history. Methods dependent on the SMC detected to some degree an increase in N_e toward the present, with PSMC and MSMC unrealistically inferring present day population sizes in the millions (fig. 2). We believe these results to be an artifact of the SMC’s well-documented difficulty in estimating N_e toward the present (Sheehan et al. 2013; Liu and Fu 2015; Nadachowska-Brzyska et al. 2016). The limited increase inferred by SMC++ is likely due to an interaction between the SMC and SFS.

Recommendations for Practical Application and Interpretation

Given that the SFS tends to perform better when inferring recent population history, and the SMC performs better when estimating more ancient history, we propose the following recommendation. Based on the results of the simulation portion of our study, PSMC is most reliable (i.e., least sensitive to fragmentation) between ∼100 kgbp and ∼333 gbp. In contrast, MSMC will be most reliable between ∼100 kgbp and ∼1,000 gbp, however, accuracy is low, with parameter estimates consistently overestimated. Based on our results, we recommend that results of MSMC for times older than 10 kgbp should not be interpreted. SMC++ will be most reliable between ∼100 kgbp and ∼166 gbp. However, PSMC is still more accurate overall across this time interval than SMC++. Broadly speaking results of SMC methods for times more recent than 333–166 gbp should be interpreted with caution. For more ancient inferences, we recommend cautious interpretation of estimations later than 10 kgbp for MSMC, 17 kgbp for SMC++, and 33 kgbp for PSMC. Although present-day parameter estimates are close to their true values (albeit reliably overestimated), the population size trends are either uninformative or misleading. The Stairway Plot is the only method tested that accurately reconstructs the shape of very recent (∼16 gbp–present) population size fluctuations. However, due to inconsistency between the results of Stairway Plot analyses conducted using simulated and empirical data, we recommend careful interpretation of these results.

Summary of Devil Demographic History

The Tasmanian devil has experienced staged population declines since the penultimate deglaciation, ∼143 kybp, likely due to a combination of: repeated isolation and reconnection to the Australian mainland, changing climate that caused periods of cooling and later drought that restricted Tasmanian devil and prey distributions, competition and potential persecution by Aboriginal Australians, with recent culling by European colonists. Although inferences made by the Stairway Plot more closely reflect recent population size trends, parameter estimation should be interpreted with caution. Specifically, the results of our simulations indicate that the Stairway Plot tends to underestimate contemporary population sizes, despite most-accurately estimating the shape of recent trends. Regardless, it appears that the long history of decline in devils may have led to increased susceptibility to transmissible cancers (DFTD, DFT2) due to a loss of standing genetic diversity.

Importantly, sampling design may impact demographic reconstruction. Previous work has explored the factors that may confound the estimation of N_e by genomic methods such as those that rely on the SMC and SFS. For instance, population structure and gene flow among demes has been repeatedly shown to bias demographic inference by generating spurious patterns of population expansion (Nielsen and Beaumont 2009; Chikhi et al. 2010; Paz-Vinas et al. 2013; Mazet et al. 2015, 2016). Specifically, sampling two demes with ongoing gene flow tends to generate a pattern of population expansion (Mazet et al. 2016). In contrast, sampling only one of these demes leads to the erroneous inference of a population bottleneck (Nielsen and Beaumont 2009; Mazet et al. 2016). In response, some methods have attempted to account for this by estimating independent demographic histories and timing of population splits for samples obtained from different present-day populations (e.g., MSMC, SMC++). However, we know of no method/approach that accounts for the reverse of this situation, in which ancient populations merged into a single contemporary population.

Our sampling of devils is not spatially representative of the whole of the species’ range. Our 12 samples are limited to the North-West-central portion of Tasmania, whereas devils are broadly distributed across the extent of the island, excepting the South-eastern third. Importantly, our sampling occurs within a single population according to previous population genetic studies (Hendricks et al. 2017). Consequently, there is little concern that we have sampled individuals from two contemporary demes. The potential impacts of this sampling design on our results are 2-fold. First, this limits the confounding effects of conflicting demography on our inferences. Second, it limits our ability to detect the demographic consequences of DFTD, which was first detected in the north-eastern corner of Tasmania ∼22 years (∼7 generations) ago (Hawkins et al. 2006). As the disease only arrived in the region from which our samples were taken ∼8 years (<3 generations) ago, we are unable to detect the consequences of this disease on genetic diversity in our samples, despite well described population declines (Lachish et al. 2007; Hamede et al. 2015). However, this fact also means that our results are unlikely to be confounded by recent changes driven by DFTD.

With the exception of the one sample inferred by PSMC to have a slightly different history (fig. 2B), estimates for each sample (faint lines) largely overlap. As is, we have no hypothesis as to why this individual, one of the seven samples obtained from Takone (the most north-westerly population sampled) might have a different demographic history. The close geographic proximity of the remaining 11 samples likely explains the minimal variation in demographic history among them however.

Note also that our estimates represent a lower bound for timing of the population size history, because they are derived from a genome-wide estimate of the substitution rate between the Tasmanian devil and three other closely related marsupials (Macrophus eugenii, Phascolarctos cinereus, Monodelphis domestica). Consequently, the true value of the mutation rate may be even greater, pushing the timing of population size shifts toward the present. In fact, such a shift would result in even greater concordance between biogeographic/anthropogenic events and observed population size history.

Importantly, we observe a number of consistencies between results obtained using empirical and simulated data. First, MSMC infers a spurious, enormous population size fluctuation between ∼100–30 kybp. MSMC also consistently estimates larger population sizes than other methods. We can confidently say that this is not due to the approach used for genotyping (samtools), as we also used these data for analysis with PSMC, yet PSMC performed much more consistently and in a manner consistent with results obtained using simulated data. Were the genotyping approach to be the cause of MSMC’s poor performance, we would expect similarly poor performance from PSMC, but we do not observe this. Second, SMC++ appears to consistently (using both empirical and simulated data) estimate a population size decrease and recovery between 500 ybp and the present day. In the case of estimates obtained using simulated data, this decline is inconsistent with the history used to simulate the input sequences. Thus, we do not interpret literally these recent declines inferred by SMC++ using empirical data. Lastly, the Stairway Plot infers a demography incompatible with those inferred using other methods. Again, this holds for both inferences obtained using empirical and simulated data.

Genome Assembly

We extracted ∼4 μg of high-quality genomic DNA from a captive, adult female Tasmanian devil from the San Diego Zoo. Note that, this was a different individual than was sequenced for the 2012-Assembly, however, the individuals came from the same geographic regions (North-eastern Tasmania) and a single interbreeding population (Hendricks et al. 2017). Consequently, we anticipated minimal differences outside of slight variation in the number/distribution of variable sites when aligning resequenced genomes to the two assemblies. Two sequencing libraries were constructed using the Chicago library protocol by Dovetail Genomics (Santa Cruz, CA) and run on a single lane of paired-end 150 bp sequencing on an Illumina HiSeq X machine. Scaffolding of the assembly was conducted with the Murchison et al. (2012) genome as an input assembly in the Dovetail HiRise assembly pipeline.

Sampling and Whole-Genome Sequencing

Individuals were nonlethally trapped as previously described (Hamede et al. 2015; Margres, Jones, et al. 2018). We sequenced the genomes of 12 S. harrisii from northwestern Tasmania (fig. 1), 10 of which are included in Margres, Ruiz-Aravena, et al. (2018). Genomic DNA was extracted from ear biopsies. Whole-genomes for each individual were sequenced with 150 bp paired-end on an Illumina HiSeq X platform to $\sim 30\times$ coverage; sequencing was performed at the Northwest Genomics Center at the University of Washington (Seattle, Washington) and GENEWIZ (South Plainfield, NJ).

Alignments and Variant Calling

To obtain SNPs across all samples, data were processed as previously described Margres, Ruiz-Aravena, et al. (2018). Briefly, raw reads were merged with FLASH2 (Magnoc and Salzberg 2011) and trimmed with Sickle (Joshi and Fass 2011). Reads were aligned to the previously published reference genome (Murchison et al. 2012) (downloaded from Ensembl January 2016) as well as the new genome assembly described in this study using the BWA-MEM algorithm (Li and Durbin 2009). Duplicate reads were marked and removed using Picard (Broad Institute 2018). Joint genotyping of all 12 samples was performed in GATK according to their best practices (McKenna et al. 2010; DePristo et al. 2011) using HaplotypeCaller for each reference genome independently as previously described (Margres et al. 2017). Briefly, only variants matching any of the following criteria were retained: quality by depth (QD) $\ge$ 2.0, Phred-scaled P value (FS) $\le$ 60.0 and root mean square of the mapping quality (MQ) $\ge$ 40.0. Variants with mapping quality rank sum test approximation (MQRankSum) of 12.5, and a read position rank sum test approximation (ReadPosRankSum) of 8 were removed. We identified 2,822,764 raw SNPs and 1,194,776 raw indels when aligning to the previously published reference genome (Murchison et al. 2012). We identified 2,790,031 raw SNPs and 1,189,816 raw indels when aligning to the new genome assembly described in this study.

These raw SNPs were subsequently filtered more stringently using vcftools v0.1.14 (Danecek et al. 2011), retaining only biallelic sites matching the following criteria: genotype quality $\ge$ 60 (minGQ 60), mean depth of coverage $\ge$ 15 (min-meanDP 15), present in at least 75% of samples (max-missing 0.75). Filtering in this way reduced the number of SNPs recovered to 1,468,034, and 1,468,044 in the 2012- and 2019-Assemblies, respectively.

To obtain sample-specific SNP data sets for use with PSMC and MSMC, we used the samtools mpileup and bcftools call functions, removing sites with mean map and base qualities <20 (samtools mpileup -q 20 -Q 20 -C 50 -u). Samtools differs from GATK in that uses individual samples for genotype calls, and thus does not incorporate population allele frequencies into its genotype calls. Likewise, it makes no assumptions regarding Hardy–Weinberg equilibrium. We chose this approach, as both MSMC and PSMC are applied in our study to individual samples, and so we sought to maximize the number of unique sites per sample. Further, this samtools workflow has become standard for studies using PSMC or MSMC (see Nadachowska-Brzyska et al. 2016).

To test for differences between SNP calling procedures, we used the bcftools isec command(Li et al. 2009) to identify the number of shared and unique sites among approaches for each sample. These results are summarized in supplementary table S3, Supplementary Material online. To further test whether our results are robust to our genotype calling approaches, we applied PSMC to data sets for one sample (Sample 1, table 2) obtained using both GATK and Samtools (supplementary fig. S4, Supplementary Material online). Results were qualitatively similar and thus not discussed further. Note that results for samtools plotted in supplementary figure S4, Supplementary Material online, are the same as in figure 2.

Calculation of the Mutation Rate

To obtain an estimate of the per-generation mutation rate, a parameter needed a priori by all demographic reconstruction methods used in this study, we followed a two-stage procedure. First, genome-wide distributed estimates of the synonymous substitution rate (dS) were obtained from the 2012-Assembly of S. harrisii, and the orthologous coding sequences of the closely related marsupials for which whole genome data are available (Macrophus eugenii, Phascolarctos cinereus, Monodelphis domestica) using the method of Nei and Gojobori (1986). From this distribution of dS values, we obtained a point estimate as the median of this distribution. This rate (per-site) was rescaled to units of time by pruning the time-calibrated phylogenetic tree of Mitchell et al. (2014) to these four taxa to obtain the branch-length (in number of years) from S. harrisii to the other three taxa. This was subsequently rescaled to generations by dividing this value by three, the generation time of S. harrisii (Pemberton 1990). By dividing our median estimate of dS by this branch length (in units of generations), we obtained an estimate of the per generation substitution rate (1.394808e-08). We use this as a conservative lower bound estimation of the mutation rate.

Inference of the SFS

The Stairway Plot method requires as input a SFS. To infer a folded SFS (in which only minor alleles are considered) we used ANGSD (Nielsen et al. 2012). ANGSD implements a two-stage procedure to infer the SFS. Briefly, ANGSD first calculates individual genotype likelihoods from the.bam assemblies assuming Hardy–Weinberg Equilibrium (HWE: option -doSaf 1) and uses these to calculate the sample allele frequencies at each site. Using these genotype likelihoods, a maximum likelihood estimate of the SFS is obtained through subsequent iterative application of an Expectation Maximization (EM) algorithm (program realSFS; Nielsen et al. 2012). We filtered such that bases with meaningless map quality and poor Qscores were removed (-minMapQ 1 -minQ 20). We repeated this procedure using the resequenced genomes aligned to both the 2012- and the 2019-Assemblies.

PSMC Analysis

PSMC was the first method developed that uses the SMC to infer historical population sizes. A single input file is used as input for a given analysis. Consensus fastq files were produced using the vcf2fq command implemented by the vcfutils.pl script distributed by bcftools-1.3.1 using the snps obtained by samtools before being converted to PSMC input files using the fq2psmcfa command. Estimates were conducted using the single input file per sample. All estimates were run using 50 EM iterations.

To conduct bootstrap estimates, we first fragmented long scaffolds using the splitfa function distributed with PSMC. The splitfa function fragments long scaffolds into 5 Mb long segments. Estimates were then conducted using the psmc -b option. Ten bootstrap estimates were made per sample, resulting in a total of 120 bootstrap replicates.

MSMC Analysis

Although MSMC allows for demographic inference using multiple individuals, it is computationally limited in the number of haplotypes it can take as input (maximum of eight, or four diploid individuals). Due to the large number of possible combinations of individuals and resultant immense number of necessary bootstrap replicates for our sample size, we chose to produce estimates using each of our 12 individuals. Although MSMC requires phased data when applied to more than one sample, it does not when applied to a single sample. As such, we did not attempt to phase our data for this analysis. MSMC requires individual input files per chromosome, or in our case, scaffold. Thus, we produced an MSMC input file per individual, per scaffold with a length greater than that of the N90. Input files were generated by converting the SNPs obtained using samtools to MSMC input format using the bamCaller.py script. Estimates were conducted using the full set of input files.

To conduct bootstrap estimation, scaffolds >5 Mb were broken up into a number of smaller segments (as done by default in PSMC) such that the resultant fragments were equal in size. This resultant set of fragmented scaffolds was then sampled with replacement until a resampled genome was produced of equivalent length. As our MSMC estimates were for each individual, we conducted 10 bootstrap estimates for each individual, producing a total of 120 bootstrap replicates. All estimates were conducted using 50 EM iterations.

SMC++ Analysis

Unlike the Stairway Plot, SMC++ leverages structural/linkage-based information in its estimates of historical demography, however, it does still incorporate the SFS in its estimate. As we lack positional data relative to chromosomes for the 2019-Assembly, we treated each assembly as a series of scaffolds. Thus, although unable to investigate demographic history at a chromosome by chromosome basis, we can be sure that any observed runs of homozygosity are real and not a consequence of missing data among scaffolds.

SMC++ leverages both the SFS as well as the SMC and requires a single input file per scaffold. Rather than calculating the SFS relative to the reference genome, SMC++ uses a single “distinguished individual” selected from the pool of samples against which all other samples are compared. As the choice of distinguished individual under this model may impact demographic inference, we generated 12 SMC++ files per scaffold, each differing in which individual is distinguished, leaving the other 11 as “undistinguished.” Consequently, use of all scaffolds would result in the need to specify as input ∼420,000 input files in the case of the 2012-Assembly. Attempts to do so resulted in the program crashing. Therefore, we chose instead to restrict our analysis to the N90 set of scaffolds. Doing so reduced the 2019-Assembly to 422 scaffolds, whereas the 2012-Assembly was only reduced to 3,263 scaffolds. To ensure consistency, we used the N90 set of scaffolds for all analyses.

We converted each scaffold into an SMC++ formatted input file using the vcf2smc script distributed with SMC++. All SMC++ input files were used together in a single run, resulting in a composite likelihood estimate, integrating across the variability generated by choice of distinguished individual. Runs were conducted assuming the previously described mutation rate, 40 knots (potential inflection points), thinning every 2,000 sites, and 50 EM iterations. Estimates were restricted to be made between 10 and 750,000 generations since present.

To bootstrap, we used the same sets of fragmented and resampled scaffolds as in our MSMC analysis. Bootstrap estimates were conducted using 25 such sets produced in this way. This number was chosen due to the extreme computational resources demanded by SMC++ (estimates for the 2019-Assembly took nearly 50 days to complete on a 28 core machine, and nearly double that for the 2012-Assembly).

Stairway Plot Analysis

Using the SFS described above, we proceeded to infer historical population demography using the Stairway Plot. All settings described here were applied to analyses of the SFS obtained from the resequenced sample aligned to either the 2012- or 2019-Assembly. To obtain estimates of N_e through time, 1,000 random subsamples of two-thirds of all sites were obtained from the empirical SFS. Individual Stairway Plots were estimated from these subsamples using a random selection of break points ranging from 3 to 24 in intervals of 3. These break points may be thought of as points at which N_e changes in time. We report the median of these estimates as our final, as well as the upper (97.5%) and lower (2.5%) bounds as our pseudo-CI.

Simulation, Fragmentation, and Analysis of Genome-Scale Sequences

To further test the robustness of the methods used herein to genome assembly contiguity, we replicated the previously described analyses on simulated genomes. We assumed a demography intermediate to those inferred from empirical data (fig. 3, black lines). To accomplish this, we used MaCS (Chen et al. 2009). Specifically, we simulated six chromosomes of equal size (480 Mb), such that the sum of their lengths was approximately equal to 90% of the devil genome (length of sequence used in the empirical component of this study). The scaled mutation rate (θ  =  4Nμ) was established using the mutation rate derived in this study, whereas the scaled recombination rate (ρ  =  4Nr) was based on the estimate of per-generation recombination rate by Feigin et al. (2018). Importantly, MaCS requires that both μ and r are scaled by the length of the sequence being simulated (in our case 480 Mb). Thus, for our simulations, θ  =  0.00697404, and ρ  =  0.004573798. We simulated 12 diploid genomes in this way to maximize comparability between our simulated and empirical results. Simulated sequences were converted to VCF format using custom scripts.

We aimed to use these simulated genomes to test the impact of genome fragmentation on demographic inference. As such, we fragmented the simulated genomes to six levels, each varying in their degree of contiguity. We produced data sets comprised of 100, 500, 1,000, 5,000, and 10,000 fragments of equal length using SnpSift split (Cingolani et al. 2012), which splits a VCF file into fragments of specified length. As we split each chromosome into an equal number of fragments, rounding became necessary, leading the true number of fragments in each data set to equal 96, 498, 996, 4,998, and 9,996 fragments. This procedure produces VCFs in which variant sites retain their original position, regardless of the length of sequence actually contained within the file. Consequently, we rescaled the positions of the sites within each VCF proportional to the length of the sequence contained in that file. Failure to do so results in the methods used herein to assume long runs of homozygosity preceding the first site.

For each data set (including the unfragmented, chromosome level sequences), we produced input files for PSMC, MSMC, and SMC++. Note that, we only applied the Stairway Plot to the unfragmented data set because fragmentation of simulated data will not impact the SFS. In contrast, increasing fragmentation of empirical genome assemblies may lead to lower rates of alignment and thus fewer SNPs identified postfiltering. PSMC input files were produced using the smc2psmc.py script distributed with SMC++. Specifically, a single-sample SMC++ input file was produced and subsequently converted to PSMC format. MSMC input files were produced using the generate_multihetsep.py script distributed with MSMC. SMC++ input files were produced using the vcf2smc command in SMC++. Lastly, the folded SFS was calculated directly from the multisample VCF of all 12 simulated, unfragmented diploid genomes using the vcftools –freq command.

We replicated the empirical analyses using each of the aforementioned data sets. Note that, we did not conduct bootstrapping of simulated data sets, due to both time and computational limitations. To determine the extent to which each method was sensitive to genome assembly contiguity, we calculated the proportional error $(=\frac{\mathit{Estimated}-\mathit{Simulated}}{\mathit{Simulated}})$ through time for each method and each fragmentation level (fig. 4). We chose to use this metric it retains the sign of error (i.e., communicates whether parameters were over- or underestimated), and because it may be intuitively interpreted in terms of percent error.