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

1Molecular 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. 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. Similarly, broad genome sampling is necessary to offset the low level of genetic divergence among individuals of recent co-ancestry and to overcome low phylogenetic signal to noise ratios characteristic of rapid radiations. 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. 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. To date, comparatively few plant or animal evolutionary genetic analyses have utilized MPS, 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; harnessing the capacity of these platforms to sequence targeted regions in multiplex; streamlining sample preparation and improving throughput; and developing accurate de novo assemblers that reduce reliance upon a predefined reference sequence.

2In this paper we demonstrate the feasibility and effectiveness of MPS-based chloroplast phylogenomics for one-third of the world's pine species (Pinus), a lineage with numerous unresolved relationships based on previous cpDNA-based studies. 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 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. Three species are also represented by two chloroplast genomes each (P. lambertiana, P. thunbergii, P. torreyana).

3Assemblies 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 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.

4The 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, 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.

5Full 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 comparisons 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). 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 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).

6Of 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).

7In 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*log10 [matrix size, bp]; r2 = 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 combined - showed significantly higher resolution than the genome-wide 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.

8Previous cpDNA based estimates of infrageneric relationships in Pinus 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 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.

9From our sampling, infrageneric analyses in plants published from 2006 to 2008 were typically based on 2574 aligned bp (95% bootstrap confidence interval: 2,292, 2,864) of sequence data, evaluated 31.7 ingroup species (95% bootstrap confidence interval: 20.2, 43.2), and resolved 22.6% of nodes at ≥ 95% bootstrap support (95% bootstrap confidence interval: 18.6, 26.5). Regression analysis shows that the proportion of highly resolved nodes in these studies is significantly and positively correlated with matrix length (F1,96 = 18.032; r2 = 0.149; P < 0.0001) but not the number of included taxa (F1,97 = 0.546; r2 = 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.

10Our 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 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. 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. An unequivocal resolution of this species will likely require the inclusion of multiple nuclear loci.

11When 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 and recently in Pinus subsection Ponderosae. 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 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), 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.

12It is reasonable to ask whether increased resolution is worth the effort of assembling whole plastomes. Considering the conservative nature of bootstrap measures, systematists often accept bootstrap values of ≥ 70% as reliable indicators of accurate topology. Simulation studies, 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. 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 properties of the coalescent process. Nonetheless, phylogenetic reconstruction based on complete organellar sequences may facilitate the detection of such phenomena, by reducing errors and uncertainty due to insufficient sampling of DNA sequence.

13Plastome 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., 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. 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 upon two locus chloroplast barcode for plants claims only 72% unique identification to species level. 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. Results in this area (as well as in phylogenetic and phylogeographic analyses) will be impacted particularly if advances in target isolation and enrichment and streamlining sample preparation prove globally effective.

14DNA extraction, amplification and sequencing are described in and followed Cronn et al., with 4 bp multiplex tags, replacing the original 3 bp tags (Table 1). For one sample, P. ponderosa, additional reads from three non-multiplexed lanes of genomic DNA were also included.

15Sequence assembly and alignment are described in and followed Whittall et al.. An analysis of interspecific recombination was conducted using RDP(Recombination Detection Program) v. 3.27. 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. Each subgenus was checked for recombination events using standard settings for several recombination-detection strategies, including: RDP, GeneConv, Chimaera, MaxChi, BootScan, and SiScan. 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. 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 GenBank AY228468) is apparently erroneous. This region was removed and reference guided analysis was rerun for this amplicon.

16Aligned sequences were annotated using DOGMA (Dual Organellar Genome Annotator) with manual adjustments to match gene predictions from GenBank and the Chloroplast Genome Database http://chloroplast.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.

17Sequence 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.

18Maximum 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) and automatically determined numbers of bootstrap replicates. Bayesian inference analyses (BI) were performed using MrBayes v. 3.1.2 using the GTR+G+I model, which was selected using MrModelTest v. 2.3 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 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.

19Topological differences between the full alignment topology and each of the three other largest data partitions (full alignment without ycf1 and ycf2, and exon nucleotides both with and without ycf1 and ycf2) were tested for significance using the Shimodaira-Hasegawa test with resampling estimated log-likelihood (RELL) bootstrapping (1,000 replicates) under the GTR+G model of evolution. To further determine which topological differences were most influential, tests were repeated with the positions of topology-variable accessions alternately modified to match the full alignment topology. In total, the full alignment data set was compared to nine different topologies.

20Exon indels and stop codon shifts were mapped onto the topology determined by ML analysis of the full alignment by parsimony mapping using Mesquite v. 2.6 (Maddison and Maddison, http://mesquiteproject.org). Tests of selection for exons were performed in MEGA v. 4.0 using the codon-based Z-test for selection, with pairwise deletion and the Nei-Gojobori (P-distance) model; variance of the differences were computed using the bootstrap method with 500 replicates.

21Divergence times for four nodes with topological uncertainty (P. albicaulis - P. lambertiana N - P. parviflora, P. sibirica - P. cembra - P. koraiensis, P. krempfii-section Quinquefoliae of subgenus Strobus) were estimated according to Pollard et al.. Chloroplast mutation rate was estimated by averaging maximum and minimum mutation rates for Pinaceae chloroplast genomes from two previous studies and assuming a generation time of 50 years. Two estimates were calculated for each node using either low (10,000) or high (100,000) effective population size.

22Variable-size random subsamples of the full alignment were tested under the parsimony criteria using PAUP* v. 4.0b10 (the faststep option was used for all but the two smallest partitions due to time considerations). Eleven partition sizes were tested (2.5, 5, 10, 20, 30, 40, 50, 60, 80, 100 and 120 kb) in five replicates each, with resolution measured as the percentage of ingroup nodes produced with ≥ 95% jackknife support. Relationships between partition size and ingroup resolution were estimated using least squares regressions, and 95% confidence limits for individual points were estimated based on linear regression using SAS JMP 7.0.1 (S.A.S. Institute, Inc., http://www.jmp.com/). Our full alignment, exon nucleotides and ycf1/ycf2 partitions were analyzed under the same parsimony criteria for comparison, as were the alignments of. Accessions from Gernandt et al. and Eckert et al. were pruned to include only taxa common to our sampling; the original analysis of Wang et al. was used since this data matrix was not available for alternative phylogenetic analyses.

23We evaluated 99 phylogenetic analyses from 86 studies published between 2006 and 2008 in Systematic Botany, Systematic Biology, American Journal of Botany, Taxon, Molecular Phylogenetics and Evolution, and Annals of the Missouri Botanical Garden [see additional file 2]. Analyses were selected based on: 1) the presented phylogeny was based solely on chloroplast DNA sequence; 2) the analysis included ≥ 10 species from a monophyletic genus; 3) there were more inter- than intra-specific taxa analyzed within the genus; 4) parsimony-based bootstrap or jackknife values were presented. Ingroup branches with bootstrap support ≥ 95%, the number of ingroup taxa and the aligned base pairs used in the analysis were recorded for each case. The authors' taxonomic interpretations were accepted in instances of taxonomic uncertainty. Conspecific clades were treated as one taxon unless clearly differentiated from one another, and internal bootstrap values were disregarded. The number of branches with bootstrap support ≥ 95% was regressed both on the number of aligned base pairs and the number of taxa (both log-transformed to meet assumptions of normality and equal variances).

24Illumina sequencing reads and quality scores have been deposited in the NCBI SRA database as accession SRA009802. New sequences have been deposited in GenBank as accessions FJ899555-FJ899583.

25[GenBank FJ899555-FJ899583, EU998739-EU998746, SRA009802]

26BI: Bayesian Inference; bp: base pairs; cpDNA: chloroplast DNA; kb: kilobase; ML: maximum likelihood; MP: maximum parsimony; MPS: massively parallel sequencing.

27MP obtained and assembled the plastome data, conducted the phylogenetic and statistical analyses, conducted the meta-analysis and drafted the manuscript. RC and AL conceived the study and contributed to data collection, data analysis and manuscript writing. All authors read and approved the final manuscript.