Advancements in next-generation sequencing (NGS) technologies have permitted the assembly of large, genome-scale data sets that have shed light on the evolutionary history of many taxa (e.g., Parks et al., 2009; Moore et al., 2010; Xi et al., 2012; Eaton and Ree, 2013; Tennessee et al., 2013). For plant phylogeny, there has been a major focus on methods for chloroplast phylogenomics (e.g., Parks et al., 2009; Moore et al., 2010), although methods for collecting phylogenomic data sets from the nuclear and mitochondrial genomes have also been developed (e.g., Straub et al., 2012; Eaton and Ree, 2013). Stull et al. (2013) developed a custom RNA probe set designed to capture angiosperm plastomes via solution-based hybridization. While their capture system was broadly successful, Stull et al. (2013) found that the most variable spacer regions were often captured at much-reduced coverage compared to more conserved regions, and were sometimes missed entirely if the target taxon was phylogenetically divergent from one of the 22 plastomes used in the bait design. Moreover, the current cost of the capture probes makes this method most efficient for projects dealing with hundreds of species. Another commonly employed method for plant phylogenomic studies is genome skimming (Straub et al., 2012), which takes advantage of the fact that organellar DNA and nuclear ribosomal DNA are present at high copy numbers in genomic DNA. However, a significant limitation of this method for systematic studies is that only high-copy number regions are recovered consistently across all samples, whereas regions with lower representation are only recovered in some samples and missed completely in others (Straub et al., 2011). This can be problematic for molecular systematic studies where missing data may result in misleading phylogenetic results (Lemmon et al., 2009). Moreover, being limited to high-copy regions in the genome becomes restrictive for experimental design as it excludes putatively highly informative regions in the genome such as single-copy nuclear genes (e.g., the single-copy orthologous genes [COSII] and the pentatricopeptide repeat [PPR] gene family; Wu et al., 2006, and Yuan et al., 2009, respectively).

As an alternative, we present an NGS approach that combines long PCR and Illumina sequencing to strategically compile phylogenomic data sets for molecular systematic studies. Long PCR, or long-range PCR, uses a combination of two polymerases—a nonproofreading polymerase at high concentration and a proofreading polymerase at a lower concentration—to amplify DNA fragments that range between 3 and 15 kilobases (kb), although cases of extremely large fragments (22–42 kb) have been reported (e.g., Cheng et al., 1994). Long PCR has been used extensively in human genome projects (e.g., Craig et al., 2008) and to sequence complete mitochondrial genomes (e.g., Knaus et al., 2011; Alexander et al., 2013), using both Sanger sequencing and NGS technologies. Here, we use long PCR to generate chloroplast DNA templates for systematic studies using NGS. While we focus on whole chloroplast amplification, this approach is directly translatable to targeted studies where only particular regions of the plastome are of interest (e.g., the inverted repeat or the small single-copy region). In addition, long PCR could also be very useful for the enrichment of mitochondrial and/or nuclear regions where intron sizes are large or unknown, as well as for regions that are difficult to assemble bioinformatically, such as repetitive regions.

Table 1.

List of species included in this study, with voucher information, tissue sources, and NGS assembly statistics when available.^a

continued

Our focus on the chloroplast genome is driven by its phylogenetic informativeness at essentially all taxonomic scales and its relative ease of amplification (e.g., Downie and Palmer, 1992; Graham and Olmstead, 2000; Moore et al., 2007; Parks et al., 2009; Moore et al., 2010), which have made the chloroplast the workhorse of molecular plant systematics since the beginning of the field. Moreover, the availability of a large number of angiosperm plastome sequences had facilitated the design of potentially universal PCR primers. To test this approach, we amplified the chloroplast genomes of 30 species (17 genera) across angiosperms using a set of 58 chloroplast PCR primers that were designed to potentially be universal in angiosperms and that may work in some gymnosperm lineages.

METHODS AND RESULTS

Representatives of 17 different genera (30 spp.) spanning 12 orders of angiosperms sensu APG III (Angiosperm Phylogeny Group, 2009) were chosen to test this approach (Table 1). Special focus was given to three genera in Orobanchaceae: Lamourouxia Kunth (one species), Bartsia L. (two species), and Castilleja Mutis ex L. f. (12 species). High-quality genomic DNA was extracted from ca. 0.02 g of silica gel–dried or herbarium tissue using a modified 2× cetyltrimethylammonium bromide (CTAB) method (Doyle and Doyle, 1987), yielding 30–70 ng/µL of DNA per sample. Using the 83 plastid gene angiosperm alignments of Moore et al. (2010;  Appendix S1 (apps-d-13-00063_apps1.nex)), we developed 58 primers with a goal of maximizing universality across angiosperms (Table 2). Conserved regions for primer design were identified by eye, and the primers were tested with IDT OligoAnalyzer tools (Integrated DNA Technologies, Coralville, Iowa, USA) to ensure that melting temperatures (T _m) were greater than 50°C and that there were no significant hairpins or self-dimerization problems. From these, 16 overlapping primer combinations were chosen to amplify the entire chloroplast genome in appropriately sized, overlapping fragments, making sure to allow at least 100 bp of overlap between regions (Fig. 1, Table 2) to minimize the decrease in sequencing depth usually associated with the ∼30 bp immediately adjacent to the primer sites (Cronn et al., 2008; Harismendy and Frazer, 2009; Cronn et al., 2012).

PCRs were performed using a combination of two high-quality Taq polymerases—QIAGEN Taq DNA Polymerase (5 units/µL) and QIAGEN HotStar HiFidelity DNA Polymerase (2.5 units/µL) (QIAGEN, Valencia, California, USA)—to obtain amplification of fragments between 5 kb and 12 kb. The QIAGEN HotStar HiFidelity DNA Polymerase was diluted to 0.2 units/µL by combining 0.1 µL of 5× QIAGEN HotStar HiFidelity PCR buffer, 0.36 µL of double-deionized water (ddH₂O), and 0.04 µL of QIAGEN HotStar HiFidelity DNA Polymerase (2.5 units/µL). Each PCR had a total volume of 25 µL, was prepared on ice, and contained the following reagents: 2.5 µL of 10× PCR buffer (QIAGEN CoralLoad or colorless, with 15 mM MgCl₂), 1.0 µL MgCl₂ (QIAGEN 25 mM), 0.75 µL of deoxyribonucleotide triphosphates (dNTPs, each at 10 mM), 5.0 µL of 5× QIAGEN Q solution, 2.5 µL of both forward and reverse primers (each at 5 µM), 0.25 µL (1.25 units) of QIAGEN Taq DNA Polymerase, 0.5 µL of the diluted QIAGEN HotStar HiFidelity DNA Polymerase solution, 9 µL of ddH₂O, and 1.0 µL of DNA template. Long PCR profiles were as follows: preheat at 93°C, initial denaturation at 93°C for 3 min followed by 35 cycles of denaturation at 93°C for 15 s, annealing at 48–68°C (depending on the primer pair) for 30 s, and extension at 68°C for 5–12 min (1 min/kb of target). To assess amplification, 2 µL of the final reactions were examined on a 1% agarose gel with appropriate size standards and the final products were kept at 4°C. The complete, step-by-step long PCR protocol can be found in Appendix 1.

Table 2.

Universal angiosperm primers used for chloroplast genome amplifications. The 16 primer combinations chosen for this study are in bold with approximate amplicon sizes in kilobases (kb) indicated.^a

Fig. 1.

The final annotated chloroplast genome assembly of Bartsia inaequalis with the 16 overlapping primer combinations indicated. Note that the primer combinations for regions 11, 12, 13, and 16 amplify both inverted repeat A and B in a single reaction. Photos by Simon Uribe-Convers.

For the three genera of Orobanchaceae in which PCR optimization was performed, amplification of the fragments was straightforward and had an average success rate of 89.7% (range = 73–100%). The most difficult regions to amplify were regions 2(trnQ^(UUG)-rpoC2), 9 (petA-psbB), 10 (psbB - rps3), and 14 (trnN^(GUU)-ndhA), which are among the largest fragments (10.3 kb, 9.8 kb, 10.9 kb, and 11.2 kb, respectively; Table 2). It was possible to split region 2 into two smaller fragments, 2a (trnQ ^(UUG)-atpH: 6.3 kb) and 2b (atpF-rpoC2: 4 kb), which facilitated its amplification in several taxa. This was not the case for regions 9, 10, and 14, for which multiple long PCR experiments using varying amounts of DNA template were necessary to obtain successful amplifications. Amplification outside of Orobanchaceae was highly variable, with an average success rate of 70.8% (range = 22–100%) with regions 5, 6, 9, 10, and 11 showing the lowest success. Importantly, the results for these taxa were obtained after just two rounds of PCR where the annealing temperatures were changed to either 48°C or 55°C. Although we did not optimize the long PCRs for each group, we are confident that optimization on a per group basis (e.g., increasing template volume, altering annealing temperatures, and/or long PCR profiles) and/or the use of fresh tissue for DNA extractions would improve success rates. Furthermore, if genomic rearrangements and/or primer mismatches are present in certain groups, primer combinations other than the 16 that were used here could be tested (Table 2). Nevertheless, we successfully amplified all 16 regions in seven species, whereas in the remaining 23 species it was only possible to amplify between six (1 sp.) and 15 (8 spp.) regions (Table 1). These results translate to 21 species having at least 12 regions amplified (114.7 kb based on potential amplicon size), representing ca. 74% of the chloroplast genome when considering only one copy of the inverted repeat. Even the species with the smallest number of amplified fragments (Castilleja arvensis Cham. & Schltdl.) was represented by ∼73 kb of data, exemplifying the effectiveness of this approach.

It is notable that many of the DNAs that were tested were extracted from herbarium tissues that ranged from five to 25 yr old when isolated. In addition, we tested these primers in several species of Abies Mill. (Pinaceae; Table 1) with surprising success, amplifying between six and nine regions without any PCR optimization. We caution that our long PCR protocol works best using recent DNA extractions that have not been through multiple freeze-thaw cycles. Ideally, long PCR should be conducted using new DNA extractions that are stored at 4°C while performing experiments. Additionally, discrete PCR bands were only obtained using high-quality Taq polymerases. When conventional polymerases were used (e.g., Go Taq [Promega Corporation, Madison, Wisconsin, USA] or Top Taq [QIAGEN]), the resulting PCR products were smears rather than discrete bands and were not used for sequencing.

To confirm that our long PCR approach was compatible with NGS and that our primers would yield complete chloroplast genomes, the amplicons from each of the 15 Orobanchaceae taxa were purified by precipitation in a 20% polyethylene glycol 8000 (PEG)/2.5 M NaCl solution and washed in 70% ethanol. The amplicons were sheared by nebulization at 30 psi for 70 s, yielding an average shear size of 500 bp as measured by a Bioanalyzer High-Sensitivity Chip (Agilent Technologies, Santa Clara, California, USA). DNA normalization is a critical step when pooling samples for multiplexing in NGS; however, due to the large number of plastomes per cell and the very few samples that were being sequenced in such a high-throughput sequencing platform, no DNA quantification was made and the sheared amplicons were pooled by species at equal volume ratios. Sequencing libraries were constructed using the Illumina TruSeq library preparation kit and protocol (Illumina, San Diego, California, USA) and were standardized at 2 nM prior to sequencing. Library concentrations were determined using the KAPA qPCR kit (KK4835; Kapa Biosystems, Woburn, Massachusetts, USA) on an ABI StepOnePlus Real-Time PCR System (Life Technologies, Grand Island, New York, USA). The resulting libraries were multiplexed in one Illumina HiSeq 2000 lane (∼187.5 million reads per lane [Glenn, 2011]) at the Vincent J. Coates Genomics Sequencing Laboratory at the University of California, Berkeley, yielding ∼12.5 million 100-bp single-end reads for each taxon (GenBank Sequence Read Archive accessions: SRR1023085, SRR1023089, SRR1023095, SRR1023112, SRR1023113, SRR 1023126, SRR 1023128–SRR1023136). Average depth of coverage of our sequencing experiment was ∼8333× (taking 150 kb as the average plastome size). The results obtained here clearly do not maximize the potential of the Illumina HiSeq 2000 for plastome sequencing. To take full advantage of the large amount of data produced by a HiSeq 2000 for plastome sequencing, it would be theoretically possible to sequence ∼4170 samples per lane and still reach the 30× minimum threshold generally regarded as ideal for plastome sequencing (Straub et al., 2012). However, high-level multiplexing in NGS with this or any other high-throughput method requires careful normalization of DNA concentrations across samples and sufficient adapter barcodes; commonly used commercial kits currently offer either 96 (NEXTflex DNA Barcode kit; Bioo Scientific, Austin, Texas, USA) or 386 (Fluidigm, San Francisco, California, USA). Alternatively, one could choose to perform this type of experiment on an NGS platform that yielded a lesser amount of data, e.g., 1 million 250-bp paired-end reads on an Illumina MiSeq Reagent Nano Kit version 2, which would allow a 30× sequencing depth for 96 samples (or 50× sequencing depth for 64 samples).

Because of the high depth of coverage of our sequencing experiment, reads were cleaned at high stringency (minimum quality = 30/40, maximum number of low-quality bases per read = 5, maximum number of duplicate reads = 10, minimum number of duplicate reads = 2) and assembled against a reference genome (Sesamum indicum L., GenBank accession no. JN637766) using the Alignreads pipeline version 2.25 (Straub et al., 2011) with the following options: percent identity = medium, minimum coverage depth = 5, and single nucleotide polymorphism (SNP) minimum coverage depth = 25 with 80% of those reads supporting the SNP. The resulting assemblies had an average depth of ∼700×, an average of 0.79% bases that were masked for not reaching the minimum sequencing depth of 5×, and an average N50 of 35,053 bp (Table 1 ; contigs and ACE files deposited in the Dryad Digital Repository:  http://doi.org/10.5061/dryad.kc75n; Uribe-Convers et al., 2014). We noticed a small decrease in sequencing depth in regions immediately adjacent to some primer sites, which is a phenomenon that has been reported in the past (Whittall et al., 2010; Knaus et al., 2011; reviewed in Cronn et al., 2012). Given that our shortest overlap between amplicons is 135 bp (between regions 9 and 10; Table 2), with the rest spanning hundreds of base pairs (Table 2), and that our experiment yielded a high sequencing depth, we had no problems calling bases unambiguously (99.99% on average, Table 1). The Bartsia inaequalis Benth. assembly (Fig. 1; GenBank accession no. KF922718) was annotated using DOGMA (Wyman et al., 2004) and visualized in GenomeVx (Conant and Wolfe, 2008).

CONCLUSIONS

We present an alternative approach for systematic studies that combines long PCR and NGS to strategically compile phylogenomic data sets for molecular systematic studies. This approach is on par with genome skimming in terms of costs, but it has the advantage of being a targeted approach and has the potential to produce data more uniformly across samples, i.e., minimizing missing data across taxa. Although this approach was only tested with chloroplast data, we emphasize that the long PCR amplicons can be generated using DNA from any genome, expanding the possibilities of long PCR and NGS for molecular systematic studies. This last point is important for studies targeting the mitochondrion or low-copy regions of the genome that otherwise might be missed or not shared across all samples using genome skimming approaches. For example, this approach may be particularly useful for the enrichment of nuclear regions, where intron sizes are large or unknown.