Increasing phylogenetic resolution at low taxonomic levels using massively parallel sequencing of chloroplast genomes

Background Molecular evolutionary studies share the common goal of elucidating historical relationships, and the common challenge of adequately sampling taxa and characters. Particularly at low taxonomic levels, recent divergence, rapid radiations, and conservative genome evolution yield limited sequence variation, and dense taxon sampling is often desirable. Recent advances in massively parallel sequencing make it possible to rapidly obtain large amounts of sequence data, and multiplexing makes extensive sampling of megabase sequences feasible. Is it possible to efficiently apply massively parallel sequencing to increase phylogenetic resolution at low taxonomic levels? Results We reconstruct the infrageneric phylogeny of Pinus from 37 nearly-complete chloroplast genomes (average 109 kilobases each of an approximately 120 kilobase genome) generated using multiplexed massively parallel sequencing. 30/33 ingroup nodes resolved with ≥ 95% bootstrap support; this is a substantial improvement relative to prior studies, and shows massively parallel sequencing-based strategies can produce sufficient high quality sequence to reach support levels originally proposed for the phylogenetic bootstrap. Resampling simulations show that at least the entire plastome is necessary to fully resolve Pinus, particularly in rapidly radiating clades. Meta-analysis of 99 published infrageneric phylogenies shows that whole plastome analysis should provide similar gains across a range of plant genera. A disproportionate amount of phylogenetic information resides in two loci (ycf1, ycf2), highlighting their unusual evolutionary properties. Conclusion Plastome sequencing is now an efficient option for increasing phylogenetic resolution at lower taxonomic levels in plant phylogenetic and population genetic analyses. With continuing improvements in sequencing capacity, the strategies herein should revolutionize efforts requiring dense taxon and character sampling, such as phylogeographic analyses and species-level DNA barcoding.


Background
Molecular phylogenetic and phylogeographic analyses are typically limited by DNA sequencing costs, and this forces investigators to choose between dense taxon sampling with a small number of maximally informative loci, or genome-scale sampling across a sparse taxon sample [1][2][3][4]. Balancing these choices is particularly difficult in studies focused on recently diverged taxa or ancient rapid radiations, as taxon sampling needs to be sufficiently large to define the magnitude of intraspecific variation and the phylogenetic depth of shared alleles [5,6]. Similarly, broad genome sampling is necessary to offset the low level of genetic divergence among individuals of recent coancestry and to overcome low phylogenetic signal to noise ratios characteristic of rapid radiations [6]. Next generation DNA sequencing is poised to bring the benefits of affordable genome-scale data collection to such studies at low taxonomic levels (genera, species, and populations). Massively parallel sequencing (MPS) has increased per instrument sequence output several orders of magnitude relative to Sanger sequencing, with a proportional reduction in per-nucleotide sequencing costs [7,8]. In principle this could allow the rapid sequencing of large numbers of entire organellar genomes (chloroplast or mitochondria) or nuclear loci, and result in greatly increased phylogenetic resolution [9]. To date, comparatively few plant or animal evolutionary genetic analyses have utilized MPS [10][11][12], due to associated costs and the technical challenge of assembling large contiguous sequences from micro-reads. These barriers have been largely eliminated through four innovations: development of strategies for targeted isolation of large genomic regions [9,[13][14][15]; harnessing the capacity of these platforms to sequence targeted regions in multiplex [9,14,16]; streamlining sample preparation and improving throughput [17]; and developing accurate de novo assemblers that reduce reliance upon a predefined reference sequence [18,19].
In this paper we demonstrate the feasibility and effectiveness of MPS-based chloroplast phylogenomics for onethird of the world's pine species (Pinus), a lineage with numerous unresolved relationships based on previous cpDNA-based studies [20][21][22]. We also highlight the broad applicability of our approach to other plant taxa, and remark on the potential applications to similar mitochondrial-based studies in animals and plant DNA barcoding. Using multiplex MPS approaches, we sequenced nearly-complete chloroplast genomes (120 kilobases (kb) each total length) from 32 species in Pinus and four relatives in Pinaceae. Our sampling of Pinus includes both subgenera (subg. Pinus, 14 accessions; subg. Strobus, 21 accessions) and species exemplars chosen from all 11 taxonomic subsections [21] to evenly cover the phylogenetic diversity of the genus. Taxon density is highest for a chosen subsection (subsect. Strobus) as representative of a species-rich clade lacking phylogenetic resolution in previous studies [5,[21][22][23]. Three species are also represented by two chloroplast genomes each (P. lambertiana, P. thunbergii, P. torreyana).

Genomic Assemblies and Alignment
Assemblies in subgenus Strobus averaged 117 kb, with an estimated 8.8% missing data (compared to P. koraiensis reference); subg. Pinus assemblies averaged just less than 120 kb (6% estimated missing data, compared to P. thunbergii reference). Outgroup assemblies averaged just over 119 kb (10.4% average estimated missing data compared to P. thunbergii reference). Median coverage depth for determined positions was variable but typically high (range 21 to 156×) ( Table 1, [also see additional file 1]). Full alignment of all assemblies was 132,715 bp in length, including 62,298 bp from exons encoding 71 conserved protein coding genes (20,638 amino acids), 36 tRNAs and 4 rRNAs. A high degree of co-linearity is inferred for these genomes due to the absence of major rearrangements within de novo contigs, and by the overall success of the polymerase chain reaction-based sequence isolation strategy (indicating conservation of the order of anchor genes containing primer sites). However, minor structural changes (a tandem duplication in two species [24] and the apparent loss of duplicate copies of psaM and rps4 in P. koraiensis) could not be confirmed. No evidence of interspecific recombination was detected, consistent with the rarity of recombination in plant plastomes [25].
The aligned matrix contained 7,761 parsimony informative ingroup substitutions (4,286 non-coding positions and 3,475 coding positions) ( Table 2). Over one-half of parsimony informative sites (55.0%) in protein coding regions resided in ycf1 and ycf2, two large genes of uncertain function [26], that accounted for 22% of all exon sequence ( Figure 1A, B). No other exons in the pine plastome exhibit such a disproportionate number of parsimony informative sites ( Figure 1C). These loci have an elevated nonsynonymous substitution rate (Table 3) and appear to have a substantial number of indels in Pinus, although it was not possible in many cases to confidently score indels in these loci due to the inherent limitations of reference-guided assembly of short reads in length variable regions. Start codon position, overall length and stop codon positions were nonetheless largely preserved in these loci across the genus. In addition to substitutions in exons, 48 ingroup exon indels and 23 ingroup stop codon shifts were identified in 26 loci.

Phylogenetic Resolution in Non-Random and Randomized Data Partitions
Full alignment partitions yielded a higher proportion of highly supported nodes, with 88 to 91% (29 to 30/33) of ingroup nodes resolved with bootstrap support ≥ 95% in likelihood analysis. The four largest data partitions tested (full alignment and concatenated exon nucleotides, both with and without ycf1 and ycf2) yielded results that were topologically identical with the exception of four taxa (P. albicaulis, P. krempfii, P. lambertiana N, P. parviflora) (Figures 2 and 3). In addition, support for the branching order of P. cembra, P. koraiensis and P. sibirica was low in full alignment partitions. Topological differences were found to be significant according to Shimodaira-Hasegawa com-parisons of the full alignment topology to two of the other major partitions (full alignment and exon nucleotides without ycf1 and ycf2). Trends in significance were most strongly influenced by the two alternative positions of P. krempfii (Figure 2 vs. Figure 3A, C; Table 4). With the exception of P. krempfii, areas of topological uncertainty reside in a single clade that historically has lacked internal resolution (subsection Strobus) [20][21][22]. Coalescent estimations suggest that these poorly resolved subsection Strobus haplotypes diverged in rapid succession relative to the age of their shared nodes (0.009 to 0.44 coalescent units, or ca. 90,000 to 450,000 years) ( Table 5). A putative chloroplast capture event in P. lambertiana previously documented [5] was also supported with whole-plastome results. Substantial resolution was achieved in analyses of ycf1 and ycf2 data partitions, however we observed several topological differences from the full alignment with high support (primarily involving the species discussed above) ( Figure 4).
Of the 71 exon coding indels and stop codon shifts identified, 35 mapped unambiguously to monophyletic groups (that is, no accessions in a group were missing data for that event) ( Figures 5 and 6). All of these groups had strong support in nucleotide-based phylogenetic analyses (100% likelihood and parsimony bootstrap support). The remainder of these events were primarily either putatively monophyletic (missing data in one or more members of a clade) or showed strong evidence of homoplasy ( Figures  5 and 6). In parsimony analyses of variable-sized jackknife samples of our full alignment, nodal support showed a strong positive correlation with the length of the nucleotide matrix (proportion nodes ≥ 95% = -1.0808 + 0.38497*log 10 [matrix size, bp]; r 2 = 0.915, P < 0.0001) ( Figure 7A). Resolution of full alignment and exon nucleotide partitions was indistinguishable from random jackknife samples of comparable size, indicating similar phylogenetic content of these partitions and corresponding similar-sized random genomic subsamples. Partitions consisting of ycf1 and ycf2 -in particular ycf1, and ycf1 and ycf2 combinedshowed significantly higher resolution than the genomewide average ( Figure 7A). The concatenated partition ycf1 + ycf2 (13.1 kb; 77.4% nodes ≥ 95% bootstrap support) yielded only slightly less phylogenetic resolution than all exons combined (62.3 kb; 80.6% nodes ≥ 95% bootstrap support) in parsimony analysis.

Comparisons to Previous Pinus Phylogenies
Previous cpDNA based estimates of infrageneric relationships in Pinus [20][21][22] sampled the same species and/or lineages as our study, and inferred relationships using 2.82 to 3.57 kb of chloroplast DNA. Results of these studies are largely consistent with our results, although highly supported nodes (≥ 95%) accounted for only 13 to 23% of the total ingroup nodes (23% to 42% if [20,21] adjusted to match our species composition). The empirical results of these studies fell within or close to the 95% prediction intervals established from our jackknife resampling response from our full genome alignment ( Figure  7A), indicating that the loci used in prior studies (primarily rbcL and matK) are similarly informative as a comparable sample of random nucleotides from the chloroplast genome. . Regression analysis shows that the proportion of highly resolved nodes in these studies is significantly and positively correlated with matrix length (F 1,96 = 18.032; r 2 = 0.149; P < 0.0001) but not the number of included taxa (F 1,97 = 0.546; r 2 = 0.006; P = 0.461), although there was a negative trend in the latter ( Figure 7B, C). Our current sample size is typical in the number of taxa sampled, but both matrix length (132.7 kb) and the proportion of highly bootstrap-supported nodes (84.8% parsimony, 90.3% maximum likelihood) were substantially higher.

Discussion
Our results highlight that whole plastome sequencing is now a feasible and effective option for inferring phylogenies at low taxonomic levels. Compared to previous chloroplast-based phylogenetic analyses in Pinus, our data matrix contained approximately 60 times more phylogenetically informative characters resulting in an approximately two-to four-fold increase in the proportion of highly resolved nodes (after adjusting results of previous studies to match our species composition) ( Figure 8, Table 2). An important question arising from these comparisons is whether the difference in resolution is entirely attributable to the increase in nucleotides, or whether the genomic partitions sequenced in prior studies were less informative on average than the rest of the genome. In fact, the resolution provided by loci used in previous Pinus studies is indistinguishable from or slightly greater than that of comparably sized random genomic subsamples from our full alignment. Combined with the strong correlation between resolution and the size of random genomic subsample, this suggests that the increase in resolution in this study is primarily due to the increase in Data from Gernandt et al. [21] and Eckert and Hall [20] pruned to include only ingroup species and outgroup genera common to our study. (PI = parsimony informative.) matrix length. This is further supported by a significant relationship between resolution and matrix length in a broad sampling of chloroplast-based infrageneric phylogenies. Based on these results, we predict that whole-plastome analysis will yield similar gains in phylogenetic resolution not only in the genus Pinus but for most land plant genera. On the other hand, it is apparent that even the entire chloroplast genome may be insufficient to fully resolve the most rapidly radiating lineages. In this regard, our results are reflective of previous analyses of ancient rapid radiations wherein nodal resolution does not scale proportionately to the length of sequence analyzed [27,28]. Notably, the position of P. krempfii was significantly different between the four largest data partitions (Table 4), even though this species does not appear to be associated with a rapid radiation (Table 5). This result is not completely unexpected, as this species has previously been difficult to place phylogenetically [29,30]. An unequivocal resolution of this species will likely require the inclusion of multiple nuclear loci [30].
When considering recent divergence, the disproportionately high mutation rate in ycf1 (and ycf2, to a lesser extent) demonstrated here is of importance, and mirrors findings in other plant taxa [31,32] and recently in Pinus subsection Ponderosae [33]. These loci should be informative for phylogenetic studies in recently-diverged clades or in population-level studies in a range of plant species. Discretion is advised, however, as ycf1 (and possibly ycf2) appears to be a target of positive selection at least in Pinus and may reflect adaptive episodes rather than neutral genealogies. In likelihood analyses of ycf1 and ycf2, we Length and information content of 71 exons common to Pinus accessions sampled in this study  observed several topological differences from the full alignment at the subsectional level, further demonstrating that caution must be taken in drawing phylogenetic conclusions from these two loci. Although we were able to confidently score small structural changes (indels and stop codon shifts) for all other exons, it was not possible to score indels for ycf1 and ycf2 due to the apparent high rate of indel formation in these loci. In all other loci examined, small structural changes only delineated clades with concurrent high support from nucleotide-based analyses (both in present study and [20][21][22]), and thus are likely to be of limited use in species or population level discrimination. It is not clear whether this will also be the case in ycf1 and ycf2.
It is reasonable to ask whether increased resolution is worth the effort of assembling whole plastomes. Considering the conservative nature of bootstrap measures [34][35][36][37], systematists often accept bootstrap values of ≥ 70% as reliable indicators of accurate topology [36]. Simulation studies [34], however, have demonstrated greatly increased accuracy (approximately 42×) with bootstrap values ≥ 95% versus ≥ 70%, and the initial formulation of the phylogenetic bootstrap used ≥ 95% as the threshold for topological significance [38]. Our results similarly support using a 95% bootstrap support cutoff for conclusive evidence as in both areas of topological differences, more than one clade received bootstrap support ≥ 70% by analysis of alternate data partitions. It is probable that conflicting topologies with ≥ 70% but < 95% bootstrap support accurately reflect data partitions yet may not represent the plastome phylogeny, and here the use of entire organelle genomes makes it possible to adopt more conservative criteria of nodal support. There are further biological reasons why an organellar phylogeny (essentially a single-gene estimate) may not accurately represent the organismal phylogeny; these include interspecific hybridization, incomplete lineage sorting, and stochastic proper-

P. krempfii topologies P. albicaulis, P. lambertiana N, P. parviflora topologies P-value
Phylogenetic relationships of 35 pines and four outgroups as determined from ycf1 and ycf2 partitions Phylogenetic distribution of exon coding indel mutations in sampled Pinus accessions Figure 5 Phylogenetic distribution of exon coding indel mutations in sampled Pinus accessions. Exon names given above boxes, size of indel (bp) and polarity ("+" = insertion, "-" = deletion) given below boxes. Polarity of events determined by comparison to most distant outgroups. Due to the apparent high rate of indel formation in ycf1 and ycf2, these loci were not able to be confidently scored for indels and are not included in this diagram. Events for only the first copy of psaM are reported. Branching order of tree corresponds to RAxML analysis of complete alignment. Diagonal lines represent putative reversals of indel events. * indicates missing data for one or more accessions of clade. Thin internal branches correspond to ML bootstrap support < 95% or topological difference in four largest data partitions (full alignment and exon nucleotides, with and without ycf1 and ycf2).  Figure 6 Phylogenetic distribution of stop codon mutations in sampled Pinus accessions. Exon names given above boxes, amino acid shift relative to stop codon position in outgroups given below boxes. Polarity of events determined by comparison to most distant outgroups; "+" signifies extension of coding region due to stop codon mutation, "-" signifies shortening. The value of zero for the psbH-and psaM-associated events corresponds to events that alter the original stop codon without altering the total number of codons in the locus. Events for only the first copy of psaM are reported. Diagonal line represents a putative reversal in psaJ of P. parviflora. Branching order of tree corresponds to RAxML analysis of complete alignment. * indicates missing data for one or more accessions of clade. Thin internal branches correspond to ML bootstrap support < 95% or topological difference in four largest data partitions (full alignment and exon nucleotides, with and without ycf1 and ycf2).

Conclusion
Plastome sequencing is now a reasonable option for increasing resolution in phylogenetic studies at low taxonomic levels and will continue to become an increasingly simple process. As sequencers evolve to even higher capacity and multiplexing becomes routine in the near future, this will allow more extensive taxon and genomic sampling in phylogenetic studies at all taxonomic levels. It is estimated that sequencing capacity on next generation platforms will approach 100 gigabase pairs per sequencing run by the end of 2009. For perspective, this is sufficient sequence capacity to produce all 100 genus-level data sets used in our meta-analysis (including ours) at greater than 100× coverage depth in a single sequencing run. Based on the estimates of Cronn et al. [9], this sequencing capacity would also allow the simultaneous sequencing of several thousands of animal mitochondria, which could greatly benefit low-level taxonomic or population-based studies in animals that currently tend to rely on relatively short sequences from many individuals [39]. It is also clear that these improvements could enable other pursuits that are currently hindered by limited sequencing capacity, such as identification of plants by diagnostic DNA sequences (DNA barcoding). The recently agreed Relationships between matrix size and resolution in current study and meta-analysis of published studies Comparative phylogenetic resolution of Pinus species used in this study Figure 8 Comparative phylogenetic resolution of Pinus species used in this study. Resolution from A) two chloroplast loci [21] and B) our complete alignment. Distance bar corresponds to 100 nucleotide changes, and is scaled for either tree. * indicate branches with < 95% (likelihood) bootstrap support in B) (likelihood and parsimony topologies were completely congruent).
upon two locus chloroplast barcode for plants claims only 72% unique identification to species level [40]. Based on results herein, whole plastome sequences have the potential to be more highly discriminating and efficient plant DNA barcodes; in fact, the possibility of plastome-and mitome-scale barcodes has been raised previously [41]. Results in this area (as well as in phylogenetic and phylogeographic analyses) will be impacted particularly if advances in target isolation and enrichment [13][14][15] and streamlining sample preparation [17] prove globally effective.

DNA Extraction, Amplification and Sequencing
DNA extraction, amplification and sequencing are described in and followed Cronn et al. [9], with 4 bp multiplex tags, replacing the original 3 bp tags (Table 1). For one sample, P. ponderosa, additional reads from three nonmultiplexed lanes of genomic DNA were also included.

Sequence Assembly and Genome Alignments
Sequence assembly and alignment are described in and followed Whittall et al. [42]. An analysis of interspecific recombination was conducted using RDP(Recombination Detection Program) v. 3.27 [43]. Rather than using the full genomic alignment, which was too memory-intensive, concatenated nucleotide sequences for 71 exons common to all accessions were used (reflective of order on the plastome). Subgenera were investigated separately as members of opposing subgenera appear incapable of hybridization [44]. Each subgenus was checked for recombination events using standard settings for several recombination-detection strategies, including: RDP [45], GeneConv [46], Chimaera [47], MaxChi [48], BootScan [49], and SiScan [50]. A total of 24 putative recombination events were identified. On close investigation, all events involved one or more of the following: misalignment, autapomorphic noise coupled with missing data, and amplification of pseudogenes. In cases of misalignment, alignments were corrected prior to subsequent phylogenetic analyses. In cases of amplification of pseudogenes, the entire amplicon for the accession involved was turned to Ns. Inspection of the alignment also revealed that some amplicons in some accessions had failed to amplify, or amplified apparently paralogous loci (evidenced by substantially higher divergence). These regions were masked in affected accessions. The locus matK was determined to be a putative paralog in several accessions, and in four (P. armandii, P. lambertiana S, P. albicaulis, and P. ayacahuite) it was replaced with Sanger sequence [5]. We also replaced 2180 bp of poor quality sequence of the locus ycf1 in P. ponderosa with Sanger sequence. In all accessions amplified by PCR, the regions adjacent to primer sites typically had low coverage, while primers had very high coverage, thus primer-flanking regions (where problematic) and the primers were also excluded. It was also determined through Sanger sequencing that a 600 bp region of the previously published P. koraiensis plastome (positions 48808 to 49634 in Gen-Bank AY228468) is apparently erroneous. This region was removed and reference guided analysis was rerun for this amplicon.
Aligned sequences were annotated using DOGMA (Dual Organellar Genome Annotator) [51] with manual adjustments to match gene predictions from GenBank and the Chloroplast Genome Database http://chloro plast.cbio.psu.edu/. Exons were evaluated for reading frame and translations, and validity of exon mutations was judged based on presence in de novo sequence, effect on the resulting polypeptide sequence, and sequence coverage depth.

Phylogenetic Analyses
Sequence data was analyzed using all genome positions and concatenated nucleotide sequence from 71 exons common to all pine accessions; both partitions were analyzed with and without the loci ycf1 and ycf2. A relatively short (approximately 630 bp) repetitive stretch of the locus ycf1 of subgenus Strobus accessions was masked in all analyses due to alignment ambiguity. The loci ycf1 and ycf2 (ca. 14 kb combined) were also analyzed individually and together.
Maximum Likelihood (ML) phylogenetic analyses were performed through the Cipres Web Portal http:// www.phylo.org/portal/Home.do using RAxML bootstrapping with the general model of nucleotide evolution (GTR+G) [52] and automatically determined numbers of bootstrap replicates. Bayesian inference analyses (BI) were performed using MrBayes v. 3.1.2 [53] using the GTR+G+I model, which was selected using MrModelTest v. 2.3 [54] under both Aikake Information Criterion and Hierarchical Likelihood Ratio Test frameworks. Each analysis consisted of two runs with four chains each (three hot and one cold chain), run for 1000000 generations with trees sampled every 100 generations. The first 25% percent of trees from all runs were discarded as burn-in. Unweighted maximum parsimony analyses (MP) of data partitions were conducted in PAUP* (Phylogenetic Analysis Using Parsimony (*and other methods)) v. 4.0b10 [55] by heuristic search with 10 replicates of random sequence addition, tree bisection and reconnection branch swapping and a maxtrees limit of 1,000. Non-parametric bootstrap analysis was conducted under the same conditions for 1,000 replicates to determine branch support.
Topological differences between the full alignment topology and each of the three other largest data partitions (full alignment without ycf1 and ycf2, and exon nucleotides