Massively parallel sequencing of complete mitochondrial genomes from hair shaft samples

Forensic Science International: Genetics 15 (2015) 8–15

Contents lists available at ScienceDirect

Forensic Science International: Genetics journal homepage: www.elsevier.com/locate/fsig

Original Research Paper

Massively parallel sequencing of complete mitochondrial genomes from hair shaft samples Walther Parson a,b,*, Gabriela Huber a, Lilliana Moreno c, Maria-Bernadette Madel a, Michael D. Brandhagen c, Simone Nagl a, Catarina Xavier a, Mayra Eduardoff a, Thomas C. Callaghan c, Jodi A. Irwin c,* a

Institute of Legal Medicine, Innsbruck Medical University, Innsbruck, Austria Penn State Eberly College of Science, University Park, PA, USA c FBI Laboratory, Quantico, VA, USA b

A R T I C L E I N F O

A B S T R A C T

Keywords: Next generation sequencing mtDNA mtGenomes Heteroplasmy Forensic science

Though shed hairs are one of the most commonly encountered evidence types, they are among the most limited in terms of DNA quantity and quality. As a result, DNA testing has historically focused on the recovery of just about 600 base pairs of the mitochondrial DNA control region. Here, we describe our success in recovering complete mitochondrial genome (mtGenome) data (16,569 bp) from single shed hairs. By employing massively parallel sequencing (MPS), we demonstrate that particular hair samples yield DNA sufficient in quantity and quality to produce 2–3 kb mtGenome amplicons and that entire mtGenome data can be recovered from hair extracts even without PCR enrichment. Most importantly, we describe a small amplicon multiplex assay comprised of sixty-two primer sets that can be routinely applied to the compromised hair samples typically encountered in forensic casework. In all samples tested here, the MPS data recovered using any one of the three methods were consistent with the control Sanger sequence data developed from high quality known specimens. Given the recently demonstrated value of complete mtGenome data in terms of discrimination power among randomly sampled individuals, the possibility of recovering mtGenome data from the most compromised and limited evidentiary material is likely to vastly increase the utility of mtDNA testing for hair evidence. ß 2014 Elsevier Ireland Ltd. All rights reserved.

1. Introduction Mitochondrial DNA testing in forensic casework has proven extremely useful in those cases for which nuclear DNA is too scarce or degraded to yield forensically informative STR profiles. For particular sample types such as shed hairs – which by some estimates represent up to 90% of the hair samples collected at crime scenes [1] – mtDNA testing is often the only DNA-based testing option. Historically, mtDNA typing has focused on the 600 bp of hypervariable segments (HVS) I and II of the 1100 bp mtDNA control region (CR). Due to the high substitution rates and the resulting inter-individual variability, the chance that any two randomly sampled individuals match in these regions is generally less than 1% (depending on the population [2]). Nevertheless, there are particular situations in which

* Corresponding authors. E-mail addresses: [email protected] (W. Parson), [email protected] (J.A. Irwin). http://dx.doi.org/10.1016/j.fsigen.2014.11.009 1872-4973/ß 2014 Elsevier Ireland Ltd. All rights reserved.

HVS-I/HVS-II or even CR data do not provide sufficient discriminatory information to resolve distinct maternal lineages. This is the case when common HVS-I/HVS-II and CR haplotypes are encountered (e.g. up to 7% and 4.4%, respectively, for U.S. Caucasians, empop.org [3]) or when questioned and known sequences differ by only a single base. Under these circumstances, examination of CR data alone restricts the power of mtDNA testing in general. A number of studies over the past decade have demonstrated that by extending mtDNA analysis from the CR to the complete mtGenome, the power of mtDNA testing for identification purposes can be greatly improved. With complete mtGenome data, the most common U.S. Caucasian, African American and Hispanic HVS-I/HVS-II types can be reduced from, in turn, 7 to 2%, 3.6 to 1% and 2.3 to 0.5%, respectively [4,5]. In addition, mtDNA coding region SNPs have been useful for resolving missing persons cases in which more than one reference family share the same CR haplotype [5,6], the sorting and re-association of commingled skeletal remains [5], and increasing statistical support when exclusionary references are unavailable [7].

W. Parson et al. / Forensic Science International: Genetics 15 (2015) 8–15

More recently, efforts to develop reference population data are beginning to reveal the true value of complete mtGenome data for forensic purposes in terms of inter-individual variation among randomly sampled individuals [8,9]. In a recent study examining individuals from across the United States, HVS-I and HVS-II information provided resolution of 64–76% of the sampled haplotypes (depending on population) and CR data increased those values to 74–80%. However, complete mtGenome data permitted nearly full resolution in the three populations analyzed (94% in Hispanics, 98.5% in Caucasians and 98.8% in African Americans [9]. Despite these developments, coding region information is still not routinely employed in forensic casework. While this may be at least partially due to the policy issues that would need to be addressed before coding region information could be used for identification purposes, there have also been significant technical barriers to generating complete mtGenome sequences from the types of evidentiary material commonly analyzed for mtDNA. Because mtDNA data are primarily sought when the genetic material is severely limited and/or compromised, and because the acquisition of sequence data in most mtDNA casework (i.e. hair and calcified tissue specimens) is generally dependent on the recovery of numerous small amplicons, often in independent amplification reactions, there simply is not enough evidentiary material in most forensic cases to support amplification and Sanger sequencing of the entire mtGenome. Furthermore, the fact that such an approach would be cost and labor prohibitive with currently employed Sanger technologies has discouraged the development of mtGenome assays. However, as predicted nearly five years ago [10], massively parallel sequencing technologies (MPS) are now finding their way into the forensic arena [8,11–20]. These technologies largely remove the greatest existing technical and cost barriers to mtGenome sequencing of low DNA quality and quantity specimens. In fact, MPS has already been paired with assays developed by the ancient DNA community [21,22] to develop complete mtGenome haplotypes from some of the most compromised unidentified human remains [23,24]. Here, we describe the development of entire mtGenome sequence data from single hair shafts and/or shed hairs lacking follicular tissue. We demonstrate that particular hair samples yield DNA sufficient in quantity and quality to produce larger mtGenome amplicons and that entire mtGenome data can be recovered from hair extracts even without PCR enrichment. Most importantly, we describe a small amplicon multiplex assay that can be routinely applied to the types of compromised hair samples typically encountered in forensic casework. 2. Materials and methods 2.1. Large fragment method 2.1.1. Sample preparation To both demonstrate the general feasibility of recovering entire mtGenomes from hair samples via PCR enrichment, and explore the potential recovery of mtGenome data from hair extracts that have not undergone enrichment, single hair samples that included the root (but lacked visible tissue when observed under a stereomicroscope) were collected from two donors of west Eurasian ancestry for whom Sanger-based complete mtGenome sequence data were available. Three hairs in total were collected: Two dyed/treated head hairs representing one individual (Samples B1 and B2) and one un-dyed head hair representing the second individual (Sample A). Hair samples B1 and B2 had been stored at 4 8C for approximately 6 months prior to DNA extraction, while hair sample A had been stored at 4 8C for 3.5 years. All samples

9

were collected with fully informed consent. Approximately 4 cm of each hair’s root end were removed, placed into a tube containing UV-irradiated 5% Terg-a-ZymeTM and placed into a sonicating water bath at room temperature for 20 min. Hairs were then briefly rinsed by sequentially transferring the fragments (using sterile forceps) from the 5% Terg-a-ZymeTM tube to a new tube containing UV-irradiated 100% ethanol, and then to a new tube containing UVirradiated Molecular Biology Grade water. Each hair was extracted using a silica-based method adapted from a Qiagen1 userdeveloped protocol [25,26]. The method employs Qiagen1 digestion buffers (Qiagen, Germantown, MD) and the PrepFilerTM Forensic DNA Extraction Kit magnetic particles (ThermoFisher, Waltham, MA). Extract mtDNA yields were assessed using a custom real-time quantitative PCR (qPCR) assay specific for a 105 bp region of the human mtGenome [27] or the Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA) prior to any further processing. From these samples, extracts of samples A and B1 were selected for mtDNA enrichment via PCR amplification, while extract B2 was taken directly to library preparation with no intermediate enrichment step. 2.1.2. Large amplicon mtDNA enrichment Hair extracts that underwent amplification enrichment were processed according to [28]. This mtGenome amplification strategy employs eight overlapping amplicons that range in size from 1 to 3.5 kbp, with the CR amplicon being the only product smaller than 2 kbp. Based on qPCR quantitative values of the two hair extracts, mtDNA inputs per 50 ml amplification reaction were estimated to be 260,000 and 130,000 mtGenome equivalents (mtGEs). For the HL60 positive control, total genomic DNA template input was 50 pg, or 6000 mtGE based on mtDNA qPCR values. In order to gain some insight into how experimental controls may need to be handled in a forensic MPS workflow, the extraction blanks and PCR negatives processed alongside the hair samples were also taken through the subsequent library preparation and sequencing steps. 2.1.3. Library preparation Amplification products were quantified using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA). PCR amplicons were purified using Exo-SAP IT (Affymetrix, Santa Clara, CA). All samples, mDNA enriched or otherwise, were prepared for sequencing on the Illumina MiSeq instrument using the NexteraXT kit (Illumina, San Diego, CA). The eight amplicons for each sample were normalized to 0.2 ng/ ml and then pooled prior to library preparation. Library preparation and sample barcoding were performed according to manufacturer recommendations. The samples were multiplexed in the MiSeq run with a positive control, and the respective reagent blanks (RB) and negative amplification controls (NC). In order to maintain the integrity of the controls in terms of them representing the maximum possible contaminant concentration that may be present in an associated sample, the amplification products for the RBs and NCs were pooled based on the largest volume required for each DNA sample amplicon at the normalization/pooling stage. In other words, if 6 ml and 9 ml of a given amplicon from two different samples were required for downstream normalization and pooling, 9 ml (the max volume for that particular amplicon) were used for the pooled RB and NC library. 2.1.4. Shotgun genomic DNA sequencing The single hair extract selected for shotgun sequencing was concentrated down to 5 ml for input into the Nextera XT library preparation kit. Based on Qubit quantification of the original 50 ml

10

W. Parson et al. / Forensic Science International: Genetics 15 (2015) 8–15

extract, genomic DNA concentration was determined to be 0.0286 ng/ml. Thus, total genomic DNA input was estimated to be 1.43 ng. The manufacturer’s recommended guidelines for library preparation were followed with one modification: 17 instead of 12 cycles were used for library amplification. Sequencing was performed on an Illumina MiSeq Sequencer, using 2  251 cycles according to the manufacturer’s recommendations.

3500xL capillary electrophoresis sequencer (Applied Biosystems, Foster City, CA, USA), POP6 and a 50 cm capillary array. Runs comprised an injection time of 40 s, an injection voltage of 2 kV, a run voltage of 11.6 kV and a run temperature of 50 8C. Sequencing raw data were imported into Sequencher (Vs. 5.1, GeneCodes, Ann Arbour, USA) and aligned to the rCRS [29] following the updated guidelines for mtDNA typing [33,34].

2.1.5. Data analysis Read mapping and consensus variant calling were performed with the CLC Genomics Workbench 6.5.1 software (CLC Bio, Aarhus, Denmark). The default Genomics Workbench mapping and alignment parameters, which included similarity and length alignments of 90%, as well as insertion/deletion (indel) and mismatch costs of 3, were used for all samples. Alignments to the mtGenome were performed using the revised Cambridge Reference Sequence (rCRS, [29]). MtGenome variant calling was performed using the probabilistic variant caller, with the following parameters: 90% variant probability, both forward and reverse read coverage, and a minimum coverage of 10 reads. Given known difficulties with alignment of homopolymer stretches and other repeats, indels in these regions were ignored for the purposes of establishing concordance between Sanger and NGS haplotypes. For the shotgun sequencing experiment, alignments to the complete human genome were performed using the human reference genome sequence build 19, Hg19. The percentages of reads mapping to both the mtGenome and the nuclear genome were determined manually based on summary data retrieved from the MiSeq Reporter (Illumina, San Diego, CA) and CLC Genomics software packages.

2.2.2. Midi-sized amplicon mtDNA enrichment In order to address the fragmented nature of DNA routinely encountered in evidentiary hair samples, midi-sized amplicons, analogous to the previously developed midi-sized CR sequencing approach [35] spanning the entire mtGenome were developed. Primers were designed using Primer3 [36], Oligo Analyzer 3.1 [37] and UCSC in-silico PCR [38]. AutoDimer [39] was used to screen all primer pairs for potential cross-reactivity. After testing the primer pairs with single-plex PCR assays, the entire mtGenome was amplified using a total of 62 PCR amplicons (Figs. 1 and 2) organized in two multiplex PCR assays (MPX1, MPX2) each consisting of 31 non-overlapping fragments and 65 primers (one fragment of MPX1 was amplified with 3, one amplicon of MPX2 was amplified with 4 primers due to the requirement of alternative amplification primers). Primer sequences and concentrations are listed in Table S1. Both multiplexes were amplified in 50 ml reactions that included 1 BD Advantage PCR Buffer (Clontech Laboratories, Mountain View, CA, USA), 0.25 mg/ml BSA (Serva, Heidelberg, Germany), 200 mM each dNTP, 60 nM each primer, and 1 BD Advantage 2 Polymerase Mix (Clontech Laboratories, Mountain View, CA, USA). Thermal cycling was performed on a PTC240 Tetrad 2 (MJ Research, St. Bruno, Canada) starting with 95 8C for 2 min, followed by 40 cycles at 95 8C for 15 s, 55 8C for 15 s, 72 8C for 15 s, one cycle at 72 8C for 10 min and ending in a 4 8C hold. Amplicons were purified using AMPure XP beads (Beckman Coulter, Inc., Brea, CA, USA), as detailed in the Illumina protocol for

2.2. Midi-sized amplicon method 2.2.1. Sample preparation 2.2.1.1. Hair samples. Hair samples were collected with informed consent from seven volunteers for whom buccal swab reference samples were also available. Four of the hair samples were collected in 2009 and stored in a petri dish at room temperature until testing. The remaining three samples were collected recently (early 2014). The root section of each hair was excised, while the remaining hair shaft was cut into 2 cm segments. The proximal segment (the first 2 cm segment next to the root) of each hair shaft was used for mtGenome sequencing. Hair segments were incubated for 2 h at 56 8C in 900 ml hair lysis buffer [30] and 36 ml Proteinase K (20 mg/ ml) to clean the surface of potential exterior DNA contamination. Three subsequent washing steps were performed in 500 ml buffer. The hair segments were then completely digested in 180 ml hair lysis buffer, 10 ml Proteinase K and 10 ml DTT overnight at 56 8C. DNA was extracted using the EZ1 DNA Investigator Kit and the EZ1 Advanced DNA Investigator Card on the EZ1 Advanced XL (Qiagen, Hilden, Germany) according to the manufacturer’s protocol and then resuspended in 50 ml Tris buffer (10 mM Tris, 0.1 mM EDTA, pH 8.0). MtDNA was quantified using a modular real-time quantitative PCR (qPCR) approach as described in [31] using a 143 bp mtDNA fragment as target. Reactions were set up in a total volume of 20 ml and run on an AB 7500 Fast Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). 2.2.1.2. Buccal samples. MtGenomes were amplified from buccal sample extracts of the seven volunteers in two overlapping 8.5 kbp fragments following a long-range PCR protocol described in [32]. The entire mtGenome was Sanger sequenced following the protocol detailed in [32] with the exception of using the AB

Fig. 1. Scheme of the midi-sized fragments designed for the amplification of the entire mtGenome in two PCR multiplexes consisting of 31 non-overlapping amplicons each (MPX1 ‘‘inner amplicons’’, MPX2 ‘‘outer amplicons’’).

W. Parson et al. / Forensic Science International: Genetics 15 (2015) 8–15

11

Fig. 2. Fragment sizes of the amplicons targeted in the two midi-amplicon-sized multiplex reactions. PCR primers were designed to yield 300–500 bp amplicons.

human mtDNA D-loop Hypervariable Region Sequencing (Illumina, San Diego, CA, USA). PCR product quantities from hair and buccal swab samples were determined and monitored using the 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany) according to the manufacturer’s recommendations. DNA was fragmented and tagged with adapter sequences using the Nextera XT DNA Sample Preparation Kit (Illumina, San Diego, CA). Indices were added to tagmented DNA by limited-cycle PCR amplification. The library DNA was purified using AMPure XP beads and normalized according to Bioanalyzer quantitative results to pools of 2 nM. DNA libraries were denatured using NaOH and hybridization buffer was added to obtain a final concentration of 12 pM. 5% of 2 nM PhiX were added to DNA library pool as an internal quality control. MPS was performed using the Illumina MiSeq Desktop Sequencer. Libraries were loaded onto the cartridge and sequenced using 2  151 cycles with 2 base-index reads following the manufacturer’s recommendations. 2.2.3. Data analysis MiSeq-generated sequences were analyzed using the NextGENe software (SoftGenetics, State College, PA, USA) and a noncommercial analysis pipeline. The latter used BWA [40] for alignment (settings: clipping penalty settings of 0, mismatch penalty of 5 and minimum seed length of 25 and a Z-dropoff of 40), picardtools for BAM file creation and sorting [41] and GATK [42] for Indel realignment and variant calling. GATKs probabilistic UnifiedGenotyper tool was used with default settings (expect genotype_likelihoods_model set to ‘‘both’’ for SNP and Indel detection and downsampling_type to ‘‘none’’) to generate VCF files containing the haplotype information. BAM files were visualized using the Integrative Genome Viewer (Broad Institute, Cambridge, MA, USA; [43].

3. Results and discussion 3.1. Large fragments In an effort to demonstrate the feasibility of recovering entire mtGenome data from hair samples, mitochondrial DNA enrichment from the extracts of single shed hairs (these included the root, but lacked follicular material) was first performed using an eight amplicon approach that primarily targets fragments of 2– 3 kb in size [28]. While this assay is generally intended for higher quality specimens harboring pristine DNA, it has been shown to be highly effective on low quality and quantity samples [44]. As a result, two hair extracts with relatively high qPCR-determined mtDNA quantities were selected for enrichment with this method. PCR product was generated from both extracts for all amplicons. In one case, amplification was robust for all reactions, while in the second case, only six of the eight reactions yielded PCR product concentrations >1.0 ng/ml. The remaining two amplicons yielded product concentrations of 0.5 ng/ml and 0.34 ng/ml. These results were generally consistent with the qPCR values, in that the lower quantity hair extract was also the one that yielded less mtDNA PCR product overall, across all amplicons. Although the 105 bp qPCR amplicon cannot necessarily be expected to reflect the quantity of large fragments in any given extract (this is particularly true for hair, where the quality and quantity of DNA can vary dramatically both within and between individual hairs), in our experience it has been useful in providing a loose gauge of large amplicon PCR success. Generally speaking, the PCR yield in the samples reported here tended to correspond with amplicon size (i.e. the largest amplicon, at 2758 bp, produced the lowest yields for both samples tested). However, reaction efficiency was clearly also a factor.

12

W. Parson et al. / Forensic Science International: Genetics 15 (2015) 8–15

From this sequencing run of 7 samples (2 hair samples, 2 positive control DNAs, 1 Reagent Blanks and 2 PCR Negatives), over 17 million reads were generated. 16 million (96.3%) of these passed the initial MiSeq quality filter and of these, over 94% had quality values of Q30 or greater. For hair sample B1, 98.95% of the 8.5 million reads mapped to the human mtGenome, with average coverage of 68,651. For hair sample A, 99.0% of the 9.1 million reads mapped to the rCRS, with average coverage of 73,697 (Fig. 3). Though average read depth was quite high in both cases, coverage across the genome was variable, ranging from 5,132 to 175,465 reads for sample 1 (a 35-fold difference), and 13,292 to 248,514 for sample two (a 19-fold difference). In both cases, depth of coverage was greatest across the control region. This can likely be traced back to the normalization step, where the size of the amplicon (half that of the others) was not factored into the calculation. Coverage patterns across the genome were consistent between the two hair samples (Fig. 3) and also consistent with coverage maps from pristine, control DNA samples amplified and prepared for MPS using the same strategy (data not shown). Reduced coverage was apparent around polycytosine stretches at positions 309 and 573, as well as between positions 3500–3600 and 13700– 13900. While the C-tract coverage can almost certainly be attributed to known issues with the sequencing and alignment of homopolymer regions [34], the reasons underlying poor coverage in other areas are not as clear. Low coverage around position 3500 has been previously reported in samples amplified with different strategies, yet prepared and sequenced using similar methods [8,15]. But despite the consistent pattern across studies, a sufficient explanation has yet to be found. Analyses of GC content and other motifs known to cause sequencing and alignment problems revealed nothing unusual about this region [15]. In the current data, one of the reduced coverage areas may be explained by the amplification strategy itself. Low coverage near the beginning of the CR coincides with the start of one amplicon (15971) and the end of another (16402). This region does not appear to exhibit low depth of coverage when different amplification strategies are used (data not shown, [15]); and this is the only amplicon overlap region of the eight total that clearly exhibited

this feature. As with some of the other reduced coverage regions mentioned above, it is difficult to know whether these data features are the results of library preparation or data alignment/ analysis. Given the complexity of the overall MPS process in terms of laboratory steps (target enrichment, library preparation) and bioinformatics (e.g. read filtering, read alignment), it will become increasingly important moving forward to better understand which features of the final data can be attributed to sample preparation and chemistry (data production) versus informatics (data analysis). More experiments are required to address this issue. Despite this variability in read depth across the mtGenome, the NGS derived haplotypes for both samples, as well as the two positive controls, were concordant with the Sanger-generated haplotypes. The only anomaly that could not be traced back to repeat region alignment issues was the reporting of position 3107 for all samples, which should have been reported as a deletion in the NGS data but was not. This is due solely to the fact that an ‘‘N’’ is recorded at this position in the rCRS to correct for an erroneous insertion in the original Anderson reference sequence [45] and serve as a placeholder for numbering. As a result, deletions at this position in the NGS data were not considered differences by the software. Otherwise, the 34 expected variants for the first sample and the 10 expected variants for the second sample were accurately called. 3.1.1. Reagent blank and negative controls Experimental controls for the mtDNA enriched samples were assessed by comparing the total number of reads, and total number of mtGenome mappable reads, for each of the controls to the total number of reads for the samples. When considering all of the reads passing filter, 96.26% of the reads represented the four DNAs sequenced in the run (Sample B1, Sample A and two amplification positive controls), 3.70% were not associated with a sample, and the remaining 0.03% could be attributed to the reagent blank and negative controls. The two PCR negatives produced 130 and 3308 reads, of which very few – only 29 and 28, respectively – aligned to the mtGenome. In contrast, the

Fig. 3. MtGenome coverage maps for two hair extracts amplified using 1–3 kb mtGenome amplicons. Libraries were prepared with Nextera XT and sequenced on the MiSeq. For hair B1 (top panel), the average coverage across the genome was approximately 68,000 reads, with a minimum of 5,132 and a maximum of 175,465. For hair A, read coverage across the genome averaged 73,697 and ranged from 13,292 to 248,514.

W. Parson et al. / Forensic Science International: Genetics 15 (2015) 8–15

extraction blank yielded 6,156 reads, of which 5,664 (92%) aligned to the mtGenome. Of those 5,664 reads, 4,489 (79%) aligned to the CR specifically, and were clearly evident as a CR amplicon contaminant in the MPS alignment. As the mtGenome amplifications were prepared in a pre-PCR laboratory, it would seem that the contamination event most likely occurred following the enrichment PCR, at some point during library preparation. However, the sequence could not be sourced. It did not reflect the samples processed, nor any profiles of the scientists routinely working in the laboratory. Despite these clearly visible reads in the extraction blank, it is important to note that the total number of reads for the reagent blank was 1000 times less than the lowest yielding sample. Furthermore, due to the final normalization step of the Nextera XT process, the recovery of DNA in the controls is inflated in such a way that the resulting data would tend to indicate more contamination than is actually present in the corresponding samples. Because the normalization step is based on bead saturation, the contaminant molecules in the controls tend to have more opportunity to bind than that same number of contaminant molecules in a successful sample, purely as a result of lack of competition. Final sequencing of the controls may not ultimately be necessary. However, this example highlights the extreme sensitivity of MPS, and underscores the need to further explore appropriate treatment of reagent blanks/negative control, and define tolerable levels of contamination.

13

Table 1 Summary of mtDNA quantities, number of mtDNA Genome Equivalents (mtGE) used for the PCR assay using the midi-sized amplicons and mean coverage values obtained after massively parallel sequencing using the MiSeq (Illumina). Hair shaft 1 2 3 4 5 6 7

(2009) (2009) (2009) (2009) (2014) (2014) (2014)

mtGE/‘

mtGE/ml for PCR

mean CV (min/max)

3,896 1,382 646 2,183 1,768 1,352 969

19,480 6,910 6,460 21,830 8840 6760 4845

17,145 (4/34,286) 6,245 (4/12,486) 9,195 (32/18,358) 21,793 (3/43,583) 25,598 (40/51,156) 19,766.5 (7/39,526) 18,360 (23/36,697)

reference samples produced via Sanger-type sequencing (Table S2). The following exceptions were observed (Table S2):  In hair sample 4 a coverage dip in the HVS-II C-tract (303–316) did not allow proper reading of this region. In the same sample both NextGENe and the non-commercial BWA indicated the presence of a 3010A (number of reads: 1937A/48G) that could not be confirmed by Sanger-type sequencing of the hair shaft sample in that region. The reasons for this discrepancy are unknown.  A coverage dip was observed in HVS-II of sample 6; therefore positions 73 to 316 were not reliably readable.

3.2. Midi-sized amplicon multiplexes While the recovery of 2–3 kb mtDNA amplicons demonstrated that it is feasible and practical to sequence the entire mtGenome from a single hair via MPS, this type of large amplicon strategy would likely not be viable for most hair samples routinely encountered in casework. Though some samples may harbor enough DNA for such an approach (and mtDNA quantification of the extract might provide some indication of which samples those would be up front), routine casework implementation demands an approach that can accommodate degraded DNA. Reviews of casework data in our respective laboratories indicate that 4050% of hair shaft samples received and mtDNA tested produce amplicons of 350 bp or greater based on mtDNA quantification values. Although amplicons smaller than 350 bp would likely produce even greater success, the large number of primer sets required for such an approach would further complicate assay optimization and performance, particularly considering that primer design is difficult along the mtGenome with its high number of polymorphic sites. As a result, the multiplex assays presented here for degraded DNA are based on amplification products in the range of 300– 500 bp (midi-sized amplicons) as earlier introduced by [29] (Figs. 1 and 2, Table S1). The average amplicon sizes were 381.61  29.46 bp for MPX1 and 388.24  38.38 for MPX2. The seven hair samples were sequenced in two runs, containing 26 and 29 samples in total, producing about 9.5 and 20.8 million reads, respectively. Reads passing filter amounted to 95% for both runs and 70% of these reads reached Q scores above 30. Mean coverage values ranged between 6,245 and 25,598 (Table 1). Although nearly full mtGenome sequence coverage was observed for all seven hair shaft samples investigated in this study, some regions showed low coverage values down to a few reads only (Fig. 4, Table S2). As mentioned earlier, the observed imbalances are likely due to different reasons. For the midi-primer multiplexes reported here, further optimization of primer location and primer concentration may reduce imbalance and lead to more even overall sequence coverage. The MPS haplotypes derived from the hair shaft samples were generally concordant with those determined from the buccal swab

In total seven PHPs were observed in the hair shafts, four of which were found in hair sample 1. This lies in the ballpark of earlier findings were up to three PHPs were observed in the CR of hairs hafts of the same size (data not shown). In four instances PHPs differed between hair and reference samples. The use of different alignment programs resulted in comparable results in general, however slight sequence coverage differences stemming from different read filtering were observed between the two. 3.3. Shotgun sequencing of hair extracts We chose to take our efforts one step further and assess whether or not complete mtGenomes could be recovered from hair extracts without enrichment. Under the assumption that the vast majority of the DNA present in hairs is mitochondrial, and given that the majority of DNA recovered from hair is endogenous (unlike degraded bone, for instance, which is heavily contaminated with microbes), we predicted that it might be possible to recover a complete mtGenome sequence without employing an enrichment step. For those hair samples that fail to yield midi-sized PCR products, such an approach could be a straightforward means of recovering entire mtGenome data without resorting to 100 or more amplicons. Of the reads passing filter, approximately 66% of them mapped to the human reference sequence. We did not attempt to characterize the unmappable reads; but based on other studies of large tufts of hair, they were likely largely microbial [46]. For our single hair, 5,787 of the human mappable reads aligned to the mtGenome. Coverage of any given position ranged from 1 to 133, with an average of 48.63 (Figs. S1 and S2). It is important to note that due to a technical glitch during the MiSeq run, paired end sequences were not recoverable. Thus, the data presented here represent only half of the mtGenome data one might expect from a successful paired-end sequencing run. Despite this relatively low coverage (as compared to the results from the multiplexed amplicon sequences, for instance), 31 of the 34 expected mtDNA variants were called (Table S3). The three variants that were

14

W. Parson et al. / Forensic Science International: Genetics 15 (2015) 8–15

Fig. 4. MtGenome coverage maps for two of the seven hair shaft samples amplified using a midi-sized amplicon approach.

missed did not meet the filters manually set in the software’s variant caller. These filters required both forward and reverse coverage, and a minimum of 10 reads covering the position. Since positions 250 and 263, and 315.1 were covered only by reverse reads, they were not captured by the variant caller. For the 31 variants that were called, the average coverage was 48 reads, the average Phred scores for the haplotype variants was 33, and the ratio of forward to reverse reads averaged 0.40. The only exception to this was variant 10398, whose coverage was heavily biased towards reads in the reverse direction (0.13 F/R). While the recovery of complete mtGenome data directly from a single hair extract is itself encouraging, even more interesting is the fact that the 5,787 reads mapping to the mtGenome were only 0.07% of the reads that mapped to the human genome reference. In other words, 99.93% mapped to the human nuclear genome. An earlier experiment based on an alternative library preparation method (data not shown) produced similar results (0.12% of the human mappable reads aligned to the mtGenome and 99.88% aligned to the nuclear genome); however, further confirmation is still required. A somewhat higher percentage (0.2%) of mitochondrial to nuclear DNA content has been described in MPS studies of ancient hair tufts [47], but the difference could be due to sample type and quantity, to include the natural variability within and between hairs. The single hair tested in the current study included the root but no visible tissue; since the root is likely to contain more nuclear DNA than the shaft, this may be once source of discrepancy. Though we have not yet definitively established the authenticity of the nuclear DNA reads, the fact that the correct mtDNA sequence was obtained would suggest that the nuclear DNA reads are also authentic. Presumably, if contaminant DNA were present, we would have also detected it in the mtGenome alignments. 4. Conclusions We have demonstrated that via Massively Parallel Sequencing (MPS) it is now possible to recover entire mtDNA sequence information from hair samples typically encountered in forensic casework. Though mtGenomic DNA, and even nuclear genomic data, have previously been retrieved from large milligram-to-gram weight tufts of ancient hair [46–48], to the best of our knowledge, this is the first time complete mtGenomes have been developed from single hair specimens–samples that pose significant additional challenges from the standpoint of sample and DNA quantity.

We demonstrate that particular hair samples yield sufficient DNA to produce larger mtGenome amplicons and that entire mtGenome data can be recovered from hair extracts even without PCR enrichment. The shotgun data suggest that in the most difficult hair cases, those for which only the smallest of fragments (<300 bp) remain, shotgun sequencing may be a viable option (though cost prohibitive on a broad scale) for recovering complete mtGenome haplotypes. Most importantly, we demonstrate repeated success in the recovery of entire mtGenome sequence data from only 2 cm of hair shaft material by employing two small amplicon multiplexes. With all methods (large amplicon, small amplicon and shotgun), the haplotypes recovered via MPS accurately reflected the mtGenome haplotypes generated via traditional Sanger sequencing. While MPS data coverage, alignment, analysis and interpretation parameters still require optimization, our data demonstrate that it is now technically feasible to develop entire mtGenomes from the types of compromised samples that have historically yielded only HVS-I/HVS-II data or, at best, control region data. Given the demonstrated value of complete mtGenome data in terms of discrimination power, the possibility of recovering mtGenome haplotypes from the most compromised and limited evidentiary material is likely to vastly increase the utility of mtDNA testing in forensic casework. Acknowledgements The research leading to this publication was funded in part by the Austrian Science Fund (FWF) [P22880-B12] and TR L397, as well as by the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 285487. The authors would like to thank Odile Loreille, Rebecca Just and Harald Niedersta¨tter for useful discussion, as well as Anthony Onorato, Eric Pokorak, Alice Isenberg and David Valle for logistical and administrative support. This is FBI Laboratory publication number 15-01. The opinions and assertions contained herein are solely those of the authors and are not to be construed as official or as views of the US Department of Justice, its branches, the Federal Bureau of Investigation or the U.S. Government. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.fsigen.2014.11.009.

W. Parson et al. / Forensic Science International: Genetics 15 (2015) 8–15

References [1] K. Bender, P. Schneider, Validation and casework testing of the BioPlex-11 for STR typing of telogen hair roots, Forensic Sci. Int. Genet. 161 (2006) 52–59. [2] M.M. Holland, T.J. Parsons, Mitochondrial DNA sequence analysis–Validation and use for forensic casework, Forensic Sci. Rev. 11 (1999) 21–50. [3] W. Parson, A. Du¨r, EMPOP - A forensic mtDNA database, Forensic Sci. Int. Genet. 1 (2007) 88–92. [4] M.D. Coble, R.S. Just, J.E. O’Callaghan, I.H. Letmanyi, C.T. Peterson, J.A. Irwin, et al., Single nucleotide polymorphisms over the entire mtDNA genome that increase the power of forensic testing in Caucasians, Int. J. Legal Med. 118 (2004) 137–146. [5] R.S. Just, M.D. Leney, S.M. Barritt, C.W. Los, B.C. Smith, T.D. Holland, et al., The use of mitochondrial DNA single nucleotide polymorphisms to assist in the resolution of three challenging forensic cases, J. Forensic Sci. 54 (2009) 887–891. [6] K.A. Sturk, M.D. Coble, S.M. Barritt, T.J. Parsons, R.S. Just, The application of mtDNA SNPs to a forensic case, Forensic Sci. Int. Genet. Suppl. Ser. 1 (2008) 295–297. [7] J.A. Irwin, S.M. Edson, O. Loreille, R.S. Just, S.M. Barritt, D.A. Lee, T.D. Holland, T.J. Parsons, M.D. Leney, DNA identification of Earthquake McGoon 50 years postmortem, J. Forensic Sci. 52 (2007) 1115–1118. [8] J.L. King, B.L. LaRue, N. Novroski, M. Stoljarova, S.B. Seo, X. Zeng, et al., High-quality and high-throughput massively parallel sequencing of the human mitochondrial genome using the Illumina MiSeq, Forensic Sci. Int. Genet. 12 (2014) 128–135. [9] R. Just, S. Fast, M. Scheible, K. Sturk-Andreaggi, A. Ro¨ck, J. Bush, et al., Full mtGenome reference data: development and characterization of 588 forensic-quality haplotypes representing three U.S. populations, Forensic Sci. Int. Genet. (under review). [10] J.A. Irwin, W. Parson, M.D. Coble, R.S. Just, mtGenome reference population databases and the future of forensic mtDNA analysis, Forensic Sci. Int. Genet. 5 (2011) 222–225. [11] C. Van Neste, F. Van Nieuwerburgh, D. Van Hoofstat, D. Deforce, Forensic STR analysis using massive parallel sequencing, Forensic Sci. Int. Genet. 6 (2012) 810–818. [12] D.M. Bornman, M.E. Hester, J.M. Schuetter, M.D. Kasoji, A. Minard-Smith, C.A. Barden, et al., Short-read, high-throughput sequencing technology for STR genotyping, BioTechniques (2012) 1–6. [13] W. Parson, C. Strobl, G. Huber, B. Zimmermann, S.M. Gomes, L. Souto, et al., Evaluation of next generation mtGenome sequencing using the Ion Torrent Personal Genome Machine (PGM), Forensic Sci. Int. Genet. 7 (2013) 543–549. [14] M. Mikkelsen, R.F. Hansen, A.J. Hansen, N. Morling, Massively parallel pyrosequencing 454 methodology of the mitochondrial genome in forensic genetics, Forensic Sci. Int. Genet. 12 (2014) 30–37. [15] J.A. McElhoe, M.M. Holland, K.D. Makova, M.S. Su, I.M. Paul, C.H. Baker, S.A. Faith, B. Young, Development and assessment of an optimized next-generation DNA sequencing approach for the mtgenome using the Illumina MiSeq, Forensic Sci. Int. Genet. 13 (2014) 20–29. [16] M.K. Scheible, O. Loreille, R.S. Just, J.A. Irwin, Short tandem repeat typing on the 454 platform: Strategies and considerations for targeted sequencing of common forensic markers, Forensic Sci. Int. Genet. 12 (2014) 107–119. [17] S. Dalsgaard, E. Rockenbauer, A. Buchard, H.S. Mogensen, R. Frank-Hansen, C. Borsting, et al., Non-uniform phenotyping of D12S391 resolved by second generation sequencing, Forensic Sci. Int. Genet. 8 (2014) 195–199. [18] E. Rockenbauer, S. Hansen, M. Mikkelsen, C. Borsting, N. Morling, Characterization of mutations and sequence variants in the D21S11 locus by next generation sequencing, Forensic Sci. Int. Genet. 8 (2014) 68–72. [19] J. Weber-Lehmann, E. Schilling, G. Gradl, D.C. Richter, J. Wiehler, B. Rolf, Finding the needle in the haystack: differentiating identical twins in paternity testing and forensics by ultra-deep next generation sequencing, Forensic Sci. Int. Genet. 9 (2014) 42–46. [20] B.J. Bintz, G.B. Dixon, M.R. Wilson, Simultaneous detection of human mitochondrial DNA and nuclear-inserted mitochondrial-origin sequences (NumtS) using forensic mtDNA amplification strategies and pyrosequencing technology, J. Forensic Sci. 59 (2014) 1064–1073. [21] A.W. Briggs, J.M. Good, R.E. Green, et al., Targeted retrieval and analysis of five Neandertal mtDNA genomes, Science 325 (2009) 318–321. [22] J. Krause, A.W. Briggs, M. Kircher, et al., A complete mtDNA genome of an early modern human from Kostenki, Russia, Curr. Biol. 20 (2010) 231–236. [23] O. Loreille, H. Koshinsky, V.Y. Fofanov, J.A. Irwin, Application of next generation sequencing technologies to the identification of highly degraded unknown soldiers’ remains, Forensic Sci. Int. Genet. Suppl. Ser. 3 (2011) e540–e541.

15

[24] J.E. Templeton, P.M. Brotherton, B. Llamas, J. Soubrier, W. Haak, A. Cooper, et al., DNA capture and next-generation sequencing can recover whole mitochondrial genomes from highly degraded samples for human identification, Investig. Genet. 4 (2013), 26-2223-4-26. [25] QIAGEN1 User-Developed Protocol. Purification of total DNA from nails, hair, or feathers using the DNeasy1 Blood & Tissue Kit. [26] E.S. Burnside, B.J. Bintz, M.R. Wilson, Improved extraction efficiency of human mitochondrial DNA from hair shafts and its implications for sequencing of the entire mtGenome from a single hair fragment,, in: Proceedings of the American Academy of Forensic Sciences 65th Annual Meeting, Washington, DC, February 18–23 2012. [27] M.F. Kavlick, H.S. Lawrence, R.T. Merritt, C. Fisher, A. Isenberg, J.M. Robertson, B. Budowle, Quantification of human mitochondrial DNA using synthesized DNA standards, J. Forensic Sci. 56 (2011) 1457–1463. [28] E.A. Lyons, M.K. Scheible, K. Sturk-Andreaggi, J.A. Irwin, R.S. Just, A high-throughput Sanger strategy for human mitochondrial genome sequencing, BMC Genomics 14 (2013), 881-2164-14-881. [29] R.M. Andrews, I. Kubacka, P.F. Chinnery, R.N. Lightowlers, D.M. Turnbull, N. Howell, Reanalysis and revision of the Cambridge reference sequence for human mitochondrial DNA, Nat. Genet. 23 (1999) 147. [30] A. Hellmann, U. Rohleder, H. Schmitter, STR typing of human telogen hairs - a new approach, Int. J. Legal Med. 114 (2001) 269–273. [31] C.M. Bauer, H. Niedersta¨tter, G. McGlynn, H. Stadler, W. Parson, Comparison of morphological and molecular genetic sex-typing on mediaeval human skeletal remains, Forensic Sci. Int. Genet. 7 (2013) 581–586. [32] L. Fendt, B. Zimmermann, M. Daniaux, W. Parson, Sequencing strategy for the whole mitochondrial genome resulting in high quality sequences, BMC Genomics 10 (2009) 139. [33] H.-J. Bandelt, W. Parson, Consistent treatment of length variants in the human mtDNA control region: a reappraisal, Int. J. Legal Med. 122 (2008) 11–21. [34] W. Parson, L. Gusmao, D.R. Hares, J.A. Irwin, W.R. Mayr, N. Morling, E. Pokorak, M. Prinz, A. Salas, P.M. Schneider, T.J. Parsons, DNA Commission of the International Society for Forensic Genetics: revised and extended guidelines for mitochondrial DNA typing, Forensic Sci. Int. Genet. 13 (2014) 134–142. [35] C. Berger, W. Parson, Mini-midi-mito: Adapting the amplification and sequencing strategy of mtDNA to the degradation state of crime scene samples, Forensic Sci. Int. Genet. 3 (2009) 149–153. [36] Primer 3: bioinfo.ut.ee/primer3-0.4.0/. [37] Oligo Analyzer 3.1: eu.idtdna.com/analyzer/Applications/OligoAnalyzer/Default.aspx. [38] UCSC In-silico PCR: genome.ucsc.edu/cgi-bin/hgPcr. [39] P.M. Vallone, J.M. Butler, AutoDimer: a screening tool for primer-dimer and hairpin structures, Biotechniques 37 (2004) 226–231. [40] H. Li, R. Durbin, Fast and accurate short read alignment with Burrows–Wheeler transform, Bioinformatics 25 (2009) 1754–1760. [41] picardtools: http://picard.sourceforge.net. [42] M.A. DePristo, E. Banks, R. Poplin, K.V. Garimella, J.R. Maguire, C. Hartl, et al., A framework for variation discovery and genotyping using next-generation DNA sequencing data, Nat. Genet. 43 (2011) 491–498. [43] H. Thorvaldsdo´ttir, J.T. Robinson, J.P. Mesirov, Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration, Briefings Bioinform. 14 (2013) 178–192. [44] R.S. Just, M.K. Scheible, S.A. Fast, K. Sturk-Andreaggi, J.L. Higginbotham, E.A. Lyons, et al., Development of forensic-quality full mtGenome haplotypes: Success rates with low template specimens, Forensic Sci. Int. Genet. 10 (2014) 73–79. [45] S. Anderson, A.T. Bankier, B.G. Barrell, M.H. de Bruijn, A.R. Coulson, J. Drouin, et al., Sequence and organization of the human mitochondrial genome, Nature 290 (1981) 457–465. [46] M. Rasmussen, Y. Li, S. Lindgreen, J. Pederson, A. Albrechtsen, I. Moltke, M. Mestpalu, et al., Ancient human genome sequence of an extinct Palaeo-Eskimo, Nature 463 (2010) 757–762. [47] M.T.P. Gilbert, L.P. Tomsho, S. Rendulic, M. Packard, D.I. Drautz, A. Sher, et al., Whole-genome shotgun sequencing of mitochondria from ancient hair shafts, Science 317 (2007) 1927–1930. [48] W. Miller, D.I. Drautz, A. Ratan, B. Pusey, J. Qi, A.M. Lesk, et al., Sequencing the nuclear genome of the extinct woolly mammoth, Nature 456 (2008) 387–390.