Estimation of the Spontaneous Mutation Rate in Heliconius melpomene

1Understanding the process of spontaneous mutation is central for many of the most important questions in evolutionary genetics. The neutral nucleotide diversity expected within a species (θ) is proportional to the product of the spontaneous mutation rate per nucleotide site (μ) and the effective population size (Ne). Variation in the mutation rate is therefore expected to contribute to the large range of variation in neutral nucleotide diversity that has been observed in natural populations. Conversely, if nucleotide diversity for putatively neutral sites of a population has been estimated, and the mutation rate is known, it is possible to estimate Ne. Effective population size is an important factor determining the effectiveness of natural selection, and selection effects on diversity at linked sites could limit diversity in the genome. Species split times can be estimated using between-species neutral nucleotide divergence, which is also expected to be proportional to the mutation rate. This can be useful if fossil evidence-based dates of species divergence are not available. Estimates of the mutation rate for a range of species across the tree of life are therefore needed in order to better understand patterns of diversity in relation to Ne and the influence of natural selection on variation. However, at present, only a handful of direct mutation rate estimates are available, for a small number of model species.

2Mutation rate estimation has until recently depended on assaying rates of mutation at specific loci or on the between-species nucleotide divergence at putatively neutral sites, such as synonymous sites. There are, however, drawbacks to these approaches, including uncertainty about species divergence dates and nonneutral synonymous site evolution. This has led to efforts to directly estimate the mutation rate by sequencing mutation accumulation (MA) lines or outbred parents and their offspring. The MA line approach suffers from potential difficulties, however. For example, recessive mutator alleles might become fixed by inbreeding in the MA line progenitor. Furthermore, its practical applicability is limited, because inbred lines cannot be produced for most species. Sequencing parents and their offspring and searching for de novo mutations in the offspring are more generally applicable, but to date experiments have only been carried out in humans and Drosophila melanogaster.

3Both humans and D. melanogaster have “finished” genome sequences, but the genomes of most sequenced species are incomplete or “draft” and it is unknown whether parent–offspring sequencing can be applied in such cases. Widely used sequencing technologies produce short reads and spurious base calls arise due to the mismapping of paralogs. Here, we apply parent–offspring genome sequencing to a tropical butterfly species, Heliconius melpomene, whose genome sequence is currently draft. Heliconius melpomene has become a focal organism for genome-based studies of speciation and hybridization, and an accurate estimate of the mutation rate will have several immediate applications. We use deep sequencing of parents and 30 offspring to produce an estimate of the mutation rate that is close to one recently obtained by similar means in D. melanogaster.

4Among the 30 offspring, we sequenced 13 focal offspring at a high depth (mean = 26.3, SD = 7.1; supplementary table S1, Supplementary Material online) and 17 “bait offspring” at a lower depth (mean = 12.7, SD = 3.9). Bait offspring were used to remove regions prone to alignment errors that generate false positives by excluding sites at which any of these individuals had an alternate base call. Mutations were not called in the bait offspring nor did they contribute to the number of site-individuals. In 4,309 scaffolds of the draft genome assembly, there are 2.70 × 108 sites, of which 1.23 × 108 (46%) were estimated to be callable, yielding 1.60 × 109 site-individuals.

5We used the Genome Analysis Toolkit (GATK; ) to call mutations. We assume that reads having the alternate allele at a site present in multiple individuals are due to mismapping. This arises when a paralogous locus present in the sample, but not in the reference genome, contributes reads that map to the wrong place. We assume that such mismapping is equally likely to occur at mutated and unmutated sites.

6We applied the mutation calling rules described in Materials and Methods to the GATK genotype calls, yielding 15 candidate mutations appearing as de novo heterozygotes in up to two focal offspring (supplementary table S2, Supplementary Material online). We first examined each candidate using the Integrative Genomics Viewer to determine whether there are single nucleotide polymorphisms (SNPs) in complete association with the alternate base calls at the candidate sites, a characteristic of mismapped paralogs. Two closely linked candidates 17 bp apart on contig HE668478 (individual 110; supplementary fig. S1, Supplementary Material online) and candidates HE671028 and HE669561 (individual 110; supplementary figs. S2 and S3, Supplementary Material online) met this criterion. Furthermore, in each case reads containing the alternate allele are truncated and have unmapped mate pairs. We judged these four candidates as likely false positives caused by mismapping. We then attempted to verify the 11 remaining candidates by Sanger sequencing. Ten gave clearly interpretable chromatograms, confirming that eight are genuine mutations, and that candidates HE671439 and HE672001 are false positives (supplementary table S1, Supplementary Material online). Further attempts at sequencing the remaining candidate (HE671010) were unsuccessful, but mutant-bearing reads are well aligned and have aligned mate pairs (supplementary fig. S4, Supplementary Material online), suggesting that it is genuine. Thus, there are nine apparently genuine de novo mutations, eight confirmed by Sanger sequencing (table 1), each one affecting a single site and present in a single focal individual.

7Among the nine de novo mutations, the number of transitions exceeds the number of transversions, as is usually observed in eukaryotes. The number of mutations divided by twice the number of callable site-individuals yields an estimated mutation rate (uncorrected for false negatives) of 2.8 × 10−9 (95% confidence interval = 1.3 × 10−9–5.3 × 10−9, assuming that the number of mutations is Poisson distributed).

8To estimate the frequency of false negatives, we simulated synthetic mutations by modifying sequencing reads for randomly selected sites in the focal offspring. We realigned and analyzed the modified data using the same procedures as for the real data. Among 1,000 synthetic mutations, 456 occurred at callable sites where all other focal offspring, both parents and all bait offspring were pure. Of the callable mutations, 436 (96%) were called. The small proportion of uncalled mutations presumably reflects mutant-bearing reads mapping less frequently than reference reads (fig. 1). A corrected estimate for the mutation rate is therefore 2.9 × 10−9 (95% confidence interval = 1.3 × 10−9–5.5 × 10−9).

9We estimated the mutation rate per base pair by genome sequencing of parents and offspring in H. melpomene. The incomplete state of the genome causes difficulties in identification of de novo mutations, because paralogous reads map more frequently to the wrong location, often yielding false heterozygote calls. Disregarding impure sites affecting any bait offspring and more than two focal offspring effectively addressed this problem. It has been estimated that approximately 20% of spontaneous recessive sex-linked lethal mutations in D. melanogaster males occur as premeiotic clusters. In the present experiment, we detected no mutation clusters affecting two focal offspring, but in view of the small number of mutation events detected, the sequences of many more individuals will be needed to accurately estimate the rate of premeiotic cluster mutations in H. melpomene. The draft state of the genome precluded the detection of large-scale de novo variants, such as rearrangements and duplications, which are particularly sensitive to mismapping.

10Autosomal nucleotide diversity (π) at 4-fold degenerate sites in H. melpomene is approximately 2.4% (Martin SH, unpublished data). Assuming neutral synonymous sites evolution, and equating π to 4Neμ, Ne for the species is therefore approximately 2 million. This will be an underestimate if selection reduces diversity at 4-fold sites. However, estimates of Ne for D. melanogaster based on this approach are of similar magnitude, but they are orders of magnitude smaller than both species’ census population sizes. Similar diversities and effective population sizes are consistent with the small range of genetic diversity across eukaryotes, suggesting a role for processes such as genetic draft limiting diversity.

11Estimates of μ can also be used to date species divergences, assuming that neutral nucleotide divergence d = 2μt, where t is the divergence time in generations. For example, synonymous divergence between H. melpomene and its Müllerian comimic H. erato corrected for diversity within H. melpomene is 14% (Martin SH, unpublished data), yielding t = 23 million generations, which will be an underestimate if selection reduces synonymous site divergence. Assuming four generations per year, the divergence date is approximately 6 Ma, which is somewhat more recent than that estimated from a fossil-calibrated phylogeny of 10–13 Ma. Although our data suggest that current estimates of the age of the Heliconius radiation are approximately correct, further work will be required to reconcile these estimates. It remains to be seen whether the hypothesis that the early radiation of Heliconius coincided with a time of rapid uplift in the Andes about 10 Ma is supported.

12This is only the second direct estimate of the mutation rate per base pair in insects, and the first in Lepidoptera. There have been several direct estimates in D. melanogaster, by Denaturing High Performance Liquid Chromatography, by whole-genome sequencing of MA lines, and most recently by parent–offspring sequencing. There is significant variation among these estimates, but most are close to 3 × 10−9, which is remarkably close to our estimate of 2.9 × 10−9 for Heliconius. We have demonstrated that it is possible to estimate the mutation rate by offspring–parent genome sequencing for the case of a draft genome sequence. It should soon be possible to address the question of whether this lack of variation in the mutation rate extends to other arthropod groups whose draft genome sequences are now becoming available.

13The cross was previously used to produce chromosomal scaffolds of the H. melpomene genome (, supplementary material section S4). After four generations of inbreeding, a male H. melpomene melpomene from the same lineage as for the H. melpomene reference genome was crossed with a female H. melpomene rosina from a laboratory strain established from Gamboa, Panama. DNA from two F1 parents and 30 of their F2 offspring was extracted using the DNeasy Blood and Tissue Kit (Qiagen). Illumina TruSeq libraries (300 bp insert size) were sequenced using 100-bp paired-end reads on an Illumina HiSeq2500.

14Reads for parents and offspring were aligned to version 1.1 of the H. melpomene genome (available on Ensembl and from http://butterflygenome.org, last accessed November 5, 2014) using SMALT version 0.7.0.1 with default options. Output sequence alignment/map (SAM) files were converted to binary format (BAM) files, sorted and annotated with read groups using Picard version 1.84 (http://picard.sourceforge.net/, last accessed November 5, 2014).

15Each individual’s BAM file was processed to remove duplicates using Picard tools, then to realign indels using GATK. SNPs were called using the GATK UnifiedGenotyper across all individuals simultaneously to produce a variant call format (VCF) file, assuming a heterozygosity parameter of 0.01. For high read depth, genotype calls are insensitive to this parameter.

16We processed the VCF by a similar algorithm as described previously filtering sites as follows: Not marked as low quality (GATK LowQual).Read depth of both parents ≥10, both homozygous references, containing no alternate allele reads.None of the 17 bait offspring contains alternate allele reads.The genotypes of all 13 focal offspring are defined (i.e., are called by GATK).At most two focal offspring are called as heterozygous by GATK.No other focal offspring contains alternate allele reads.

17Not marked as low quality (GATK LowQual).

18Read depth of both parents ≥10, both homozygous references, containing no alternate allele reads.

19None of the 17 bait offspring contains alternate allele reads.

20The genotypes of all 13 focal offspring are defined (i.e., are called by GATK).

21At most two focal offspring are called as heterozygous by GATK.

22No other focal offspring contains alternate allele reads.

23Our method excludes sites containing alternate alleles in either parent, which precludes the identification of mutations at polymorphic sites. Assuming that mutations are not more frequent at polymorphic sites, this should reduce the number of mutations and callable sites proportionally. There was no filtering carried out on read depth of the focal or bait offspring. Heterozygotes called among the focal offspring were marked as candidate mutations.

24We estimated the proportion of false negatives (genuine mutations we failed to call) by simulating mutations in the Heliconius data, running our pipeline, and calculating the fraction of simulated mutations called. Synthetic mutations were simulated by modifying the reads overlapping a random site in a focal offspring. We sampled the number of reads to be altered from empirical distributions of numbers of nonreference base calls at heterozygous sites (see below). For each synthetic mutation, we randomly sampled a genomic position b and a focal offspring. We sampled a random integer y from the frequency distribution of nonreference base number for the individual’s read depth (e.g., see fig. 1). We then changed y reference bases to a different randomly selected base by modifying reads overlapping position b in the individual’s BAM file.

25We generated 1,000 synthetic mutations in the BAM files of focal individuals, extracted all reads from the modified BAM files, and aligned the modified reads to the reference genome by the procedure used for the original data. We then applied the mutation-calling algorithm, exactly as described above, with the exception that filters were not applied to the focal offspring carrying the synthetic mutation. The fraction of callable simulated mutated sites estimates the fraction of callable sites in the genome. Uncallable sites will include, for example, sites of low mapping quality and sites where genotypes are undefined in one or more focal offspring.

26To produce the distributions used to generate synthetic mutations, we identified a set of sites heterozygous for natural polymorphisms, regardless of the genotypes called from the sequencing data, taking advantage of the lack of recombination in Heliconius females. F2 offspring receive whole chromosomes from the F1 mother, so SNPs from the same chromosome have identical segregation patterns in the offspring, and segregation patterns for each of the 21 H. melpomene chromosomes for this cross are known. We identified SNPs from each chromosome by compiling sites called as heterozygous in the F1 mother, homozygous in the F1 father, and matching one of the chromosome segregation patterns for the bait offspring. Heterozygous focal offspring could then be identified based on segregation pattern, without reference to their sequenced genotype. We used these heterozygous focal offspring to tabulate frequency distributions of numbers of nonreference base calls for read depths 1, … 100 (see fig. 1).

27With the exception of four candidates ruled out by inspection (see Results), we checked all candidates by Sanger sequencing. We sequenced the focal individual and one control individual on both strands. If the initial sequencing failed, an alternative primer pair was tried.

28Whole-genome sequence data for the parents and the 30 offspring from the mapping cross used for this study are available from the European Nucleotide Archive, study accession PRJEB7581.

29Supplementary figures S1–S4 and tables S1 and S2 are available at Molecular Biology and Evolution online (http://www.mbe.oxfordjournals.org/).