The plastid genomes of the majority of photosynthetic angiosperms are highly conserved in structural organization, gene arrangement, and gene content (Palmer 1991; Raubeson and Jansen 2005; Wicke et al. 2011; Jansen and Ruhlman 2012; Ruhlman and Jansen 2014). Their hallmark is the presence of two large duplicate regions in reverse orientation known as the inverted repeat (IR), which separate the remainder of the genome into large single copy (LSC) and small single copy (SSC) regions. Variation in size of this molecule is due most typically to the expansion or contraction of the IR into or out of adjacent single-copy regions and/or changes in sequence complexity due to insertions or deletions of novel sequences. Of the two equimolar structural isomers existing for plastid DNA (ptDNA; Palmer 1983), the structure most commonly presented follows the convention used for Nicotiana tabacum L. (tobacco) in which one copy of the IR is flanked by genes ycf1 and trnH-GuG (and is designated as IR_A) and the other copy is flanked by genes rps19 and ndhF (and is designated as IR_B; Shinozaki et al. 1986; Yukawa et al. 2005). The junctions between the LSC region and each of these IR copies are designated as J_LA (LSC/IR_A) and J_LB (LSC/IR_B), and the junctions flanking the SSC region are designated as Jsa (SSC/IR_A) and J_SB (SSC/IR_B; Shinozaki et al. 1986). In most non-monocot angiosperm ptDNAs, J_LB lies within the ribosomal protein S10 operon in a more or less fixed position within or near the rps19 gene and J_LAa lies just downstream of LSC gene trnH-GUG. The IRs of angiosperm plastomes can fluctuate greatly in size, although small contractions and expansions of <100 base pairs (bp) are most frequent, and the positions of all four IR/single-copy junctions can vary even among closely related species (Goulding et al. 1996; Plunkett and Downie 2000). Large IR expansions (>1,000 bp) occur less frequently and outnumber large contractions (Plunkett and Downie 2000; Raubeson and Jansen 2005; Hansen et al. 2007). Because major changes in position of IR junctions can accompany structural rearrangements elsewhere in the plastid genome (Palmer 1991; Boudreau and Turmel 1995; Aii et al. 1997; Cosner et al. 1997; Perry et al. 2002; Chumley et al. 2006; Haberle et al. 2008; Guisinger et al. 2011; Wicke et al. 2011), the availability of sequence data for these genomes can help elucidate the mechanisms responsible for these large-scale IR junction shifts and other major genomic changes.

Previous ptDNA restriction site mapping studies have shown that Apiaceae (Umbelliferae) exhibit unprecedented variation in position of their LSC/IR boundary regions (Palmer 1985a; Plunkett and Downie 1999, 2000). While most umbellifer species surveyed possess a J_LB indistinguishable from that of Nicotiana tabacum and the vast majority of other eudicots, at least one expansion and seven different contractions of the IR relative to the N. tabacum J_LB were detected in 55 species, each ranging in size from ∼1–16 kilobase pairs (kb; Plunkett and Downie 2000). As examples, all examined members of the “Aegopodium group,” representatives of tribes Careae and Pyramidoptereae (Downie et al. 2010), have a ∼1.1 kb expansion of J_LB relative to that of N. tabacum. Careae and Pyramidoptereae are monophyletic sister groups (Banasiak et al. 2013), indicating that IR junction shifts have the potential to demarcate major clades within the family. Coriandrum sativum L. and Bifora radians M. Bieb., both of tribe Coriandreae, have the most contracted IRs, approximately 16.1 kb smaller than that of N. tabacum. Coriandrum also has a ∼5.7 kb insertion of unknown composition in the vicinity of the 16S rRNA gene near the terminus of the IR (Plunkett and Downie 2000), a region subsequently identified through entire plastome sequencing as a duplication of genes trnH-GUG and psbA (Peery et al. 2006). Several of these junction shifts in Apiaceae are parsimony-informative and all are restricted to the “apioid superclade” of Apiaceae subfamily Apioideae, a large, morphologically heterogeneous group of umbellifers comprising 14 tribes and other major clades of dubious relationship (Plunkett and Downie 1999, 2000; Downie et al. 2010). Plastid genome sequencing has also revealed unique, noncoding DNA in the region bounded by J_LA and trnH-GUG in several related species of Apiaceae that does not show significant similarity to any other known ptDNA sequence (Peery et al. 2007, 2011), and BLAST searches of these sequences suggested a mitochondrial to plastid transfer of DNA (Peery et al. 2011). Subsequent studies have confirmed other instances of mitochondrial DNA transfer into the plastome of Apiaceae (Goremykin et al. 2009; Iorizzo et al. 2012a,b). Thus, the Apiaceae provide a model system in which to study the mechanisms leading to large-scale expansions and contractions of the IR, with concomitant novel insertions near J_LA either the cause or consequence of these major IR junction shifts.

The Apiaceae have been the subject of numerous phylogenetic studies using ptDNA and while these studies have contributed greatly to a broad understanding of its evolutionary history, uncertainties remain with regard to many specieslevel relationships, particularly those of the largest genera within the “apioid superclade” (reviewed in Downie et al. 2001, 2010). These uncertainties are the result of a paucity of variable DNA sequence characters obtained from the relatively few plastid loci examined. The earliest of these studies focused on gene (matK exons, rbcL) and intron (rpl16, rps16, rpoC1) sequences (Downie et al. 1996, 2000; Plunkett et al. 1996a,b; Downie and Katz-Downie 1999), whereas subsequent studies examined primarily the five noncoding regions between genes trnK and psbI, and intergenic spacer regions trnH-psbA, trn D-trnT, and trnF-trnL-trnT (Calviño and Downie 2007; Downie et al. 2008; Degtjareva et al. 2009; Nicolas and Plunkett 2009; Sun and Downie 2010). Rates of nucleotide change between these regions, as well as between other plastid loci, can vary tremendously within a given taxonomic group (Shaw et al. 2005); therefore, it is useful to investigate the relative utility of previously unexplored or underutilized noncoding regions among Apiaceae plastomes for their potential in resolving phylogenetic relationships at low taxonomic levels within the family. The most variable noncoding regions uncovered can also be considered for intraspecific phylogeographic and population genetic studies.

To characterize plastome structural organization and IR evolution in the Apiales, we report the complete genome sequences of Anthriscus cerefolium (L.) Hoffm. (chervil), Crithmum maritimum L. (sea samphire), Hydrocotyle verticillata Thunb. (whorled marshpennywort), Petroselinum crispum (Mill.) Fuss (parsley), and Tiedemannia filiformis (Walter) Feist & S. R. Downie subsp. greenmanii (Mathias & Constance) Feist & S. R. Downie (giant water cowbane), and compare the results obtained to previously published plastomes of Daucus carota L. subsp. sativus Schübl. & G. Martens (domesticated carrot; Ruhlman et al. 2006) and Panax schin-seng T. Nees (Korean ginseng; Kim and Lee 2004). Anthriscus and Daucus represent tribe Scandiceae subtribes Scandicinae and Daucinae, respectively, Crithmum is classified in tribe Pyramidoptereae, Tiedemannia represents tribe Oenantheae (Feist et al. 2012), and Petroselinum is placed in tribe Apieae (Downie et al. 2010). Crithmum and Petroselinum are members of the “apioid superclade” (Downie et al. 2010), and Hydrocotyle and Panax are classified in the closely related family Araliaceae (Nicolas and Plunkett 2009). To identify the most variable noncoding loci, we compare the plastomes of the Apiaceae genera Anthriscus, Crithmum, Daucus, Petroselinum, and Tiedemannia. The seven species considered herein represent disparate clades throughout the Apiales; they also represent plastomes exhibiting major IR expansions and contractions (as assessed previously through restriction site mapping studies), as well as those possessing IR junctions in positions very similar to those of many other unrearranged angiosperm plastid genomes. This information on Apiales plastome structural organization, IR junction shifts, and noncoding loci variability will help guide subsequent systematic and evolutionary studies of this economically, ecologically, and medicinally important group of flowering plants.

Materials and Methods

Plastid Isolation, Amplification, and Sequencing —Fresh leaves of Anthriscus and Petroselinum were obtained at a local grocery store, and living plants of Crithmum (UCONN 198501242), Hydrocotyle (UCONN 198501441), and Tiedemannia (as Oxypolis greenmanii Mathias & Constance; UCONN 200202456) were obtained from the Ecology and Evolutionary Biology Greenhouse Collections at the University of Connecticut. Leaves from multiple individual plants were combined to provide enough material for the plastid isolations. Vouchers are deposited at the University of Illinois herbarium (ILL).

Approximately 5–10 g of fresh, young leaf tissue was used in each plastid genome isolation. Plastids were isolated using the sucrose stepgradient centrifugation method (Palmer 1986) and the extraction procedures outlined in Jansen et al. (2005). The plastid pellet was resuspended in wash buffer to a final volume of 2 ml. Whole plastid genome amplification was performed with rolling circle amplification (RCA) using a REPLI-g midi kit (Qiagen Inc., Valencia, California). The manufacturer's protocol for amplification of genomic DNA from cells was followed, with the following modification to improve ptDNA amplification (M. Guisinger, unpubl.): 1.5 µl of alkaline lysis solution (10 µl of solution A in Jansen et al. 2005, activated with 1.0 µl 1 M DTT) was added to 1.0 µl of isolated plastids and 4.0 µl of PBS and incubated on ice for 10 min. Subsequently, 3.5 µl of stop solution was added to this mixture. Each RCA reaction consisted of 0.5–1.0 µl of lysate, 23.0–23.5 µl of REPLI-g midi reaction buffer, and 1.0 µl of REPLI-g midi DNA polymerase. The reaction was incubated at 30°C for 16 hr and terminated at 65° for 3 min. After confirmation that the RCA was successful by running out 2 µl of the product on a minigel, another 2 µl of RCA product was then digested with 20 units of EcoRI in a 20-µl reaction, run out on a 1% agarose gel with a low DNA mass ladder, and stained and visualized with ethidium bromide. The bright, discrete banding patterns relative to the nuclear DNA background ensured that the amplification product contained sufficient quality and quantity of ptDNA to proceed with genome sequencing.

The amplified genomic DNA was sent to the W. M. Keck Center for Comparative and Functional Genomics at the University of Illinois at Urbana-Champaign, sheared by nebulization, subjected to 454 library preparation, and then sequenced using the Roche/454 Genome Sequencer (GS) FLX platform (454 Life Sciences Corp., Branford, Connecticut), following the manufacturer's protocol. For two genomes that were sequenced at two later dates (Hydrocotyle, Tiedemannia), GS FLX Titanium series reagents and plates were used.

Genome Assembly, Finishing, and Annotation —The obtained nucleotide sequence reads were assembled de novo into contiguous sequences (“contigs”) using Newbler 2.3, a DNA sequence assembly software package designed specifically for the 454 GS FLX sequencing platform, using default settings (90% minimum overlap identity). To identify contigs that were plastid sequences, all contigs were searched against NCBI's Nucleotide Collection database using BLAST'S megablast option (Altschul et al. 1990). Putative gene identifications in each plastid contig were then made using the annotation program DOGMA (Dual Organellar GenoMe Annotator; Wyman et al. 2004). These contigs were then compared with previously published reference sequences from the complete plastid genomes of Daucus carota subsp. sativus (GenBank accession number DQ898156) and Panax schin-seng (AY582139) and assembled using Sequencher 4.7 (Gene Codes Corp., Ann Arbor, Michigan). The lack of genome rearrangements, low frequency of repeated sequences, and availability of gene mapping data for three of the five Apiales plastid genomes (Plunkett and Downie 1999; Lee and Downie 2000) greatly facilitated assemblies. Gaps between contigs were filled in (or the contigs determined to abut, because there was no missing sequence between adjoining contigs) by identifying sequences at both ends of adjacent contigs and using these regions to pull out high-quality individual sequence reads from the 454 FASTA-formatted output. A minimum of 20 overlapping sequence reads in each direction was aligned using ClustalX (Jeanmougin et al. 1998) and their consensus sequence used to link the contigs. This same method was used to identify IR-single copy junctions. Sequences between a dozen pairs of contigs were confirmed by designing PCR primer pairs (18–20 bp in size) flanking gaps and Sanger sequencing using an ABI Prism BigDye Terminator v3.1 ready reaction cycle sequencing kit (Applied Biosystems, Foster City, California) and an ABI 3730x1 DNA Analyzer. Because breaks between contigs generally occurred at or near IR-single copy junctions and in noncoding regions characterized by homopolymer runs, this Sanger sequencing verified IR boundaries and confirmed homopolymeric repeats in these regions. Genomic data for three taxa (Anthriscus, Petroselinum, Crithmum) were also assembled into contigs using MIRA (Mimicking Intelligent Read Assembly), under slightly less stringent alignment settings (Chevreux et al. 2000). DOGMA was implemented to assist in annotating all genes and to identify coding sequences, tRNAs, and rRNAs using the plastid/bacterial genetic code; start and stop codons were added manually, as were small exons. Exonintron boundaries were determined by detailed comparison with other annotated genomes and individual gene sequences, especially that of N. tabacum (Shinozaki et al. 1986). The identification of tRNAs was confirmed using tRNAscan-SE (Lowe and Eddy 1997; Schattner et al. 2005). Following convention and for annotation purposes, numbering begins with the first base of the LSC region following I_LA and proceeds counterclockwise around the circular genome.

Plastome and Sequence Comparisons —Plastid genome sequences (minus one copy of the IR) of Anthriscus, Crithmum, Petroselinum, and Tiedemannia were compared to the reference Daucus using MultiPipMaker (Schwartz et al. 2003). The options “Search both strands” and “Show all matches” were chosen. REPuter (Kurtz and Schleiermacher 1999) was used to identify and locate forward (direct) and inverted (palindromic) repeat sequences. Searches of plastomes (upon the removal of one copy of the IR) were carried out with minimal repeat size ≥30 bp and a Hamming distance of 3. Because REPuter overestimates the number of repetitive elements in a given sequence by recognizing nested or overlapping repeats within a given region containing multiple repeats (Chumley et al. 2006; Timme et al. 2007), all redundant repeats were identified manually and excluded. Circular genomic maps were constructed using OGDraw (Organellar Genome Draw; Lohse et al. 2007).

The “Extract Sequences” option of DOGMA was used to extract intergenic spacers and introns from the annotated sequences in order to construct data matrices to identify highly divergent regions. The more divergent Araliaceae representatives, Panax and Hydrocotyle, were excluded from these analyses. A total of 70 noncoding ptDNA regions were compared, representing all intergenic spacers >150 bp in Daucus and introns from both LSC and SSC regions. Coding regions were excluded because they tend to provide fewer variable characters than noncoding regions; similarly, all regions contained in the IR were also excluded because many of them diverge at slower rates compared to sequences located in their adjacent single-copy regions (Wolfe et al. 1987; Perry and Wolfe 2002; Kim and Lee 2004). Four regions of combined noncoding and coding loci were also compared, as these loci are small enough to be PCR-amplified and sequenced together (i.e. trnD-trnY-trnE-trnT; trnS-psbZ-trnG-trnfM; psbB-psbT-psbN-psbH; rpl36-infA-rps8-rpll4). Sequences from the five Apiaceae genomes were aligned using ClustalX and, if necessary, manually adjusted using PAUP^* 4.0bl0 (Swofford 2002) to produce an alignment with the fewest number of changes (indels or nucleotide substitutions). The comparison of these 74 regions included numbers of constant, uninformative, parsimony-informative, and variable alignment positions, range of sequence divergence in pairwise comparisons (uncorrected “p” distance), and numbers of total and parsimony-informative indels. The number of nucleotide substitutions (uninformative + informative alignment positions) and indels for each ptDNA region was tallied; this number divided by the total aligned length resulted in the proportion of observed mutational events (or percent variability) for each region compared (Shaw et al. 2005).

Results

Plastome Assemblies —Newbler assemblies of the nucleotide sequence reads resulted in 54 (Anthriscus) to 333 (Hydrocotyle) non-redundant contigs, a portion of each of which blasted to published plastid genomes. Average read lengths varied from about 225 bp (Anthriscus, Crithmum, Petroselinum) to 336 and 381 bp (in Hydrocotyle and Tiedemannia, respectively). Average sequencing depth ranged from about 90 × (Anthriscus, Tiedemannia) to 350 × (Hydrocotyle). A comparison of assemblies of the five Apiales plastid genomes is presented in Table 1.

Plastid contigs assembled using MIRA were substantially longer than the Newbler assemblies and, in most cases, confirmed the consensus sequences that were assembled manually from the individual 454 reads. As one example, a MIRA re-assembly of Crithmum under the “accurate setting” resulted in one 68-kb contig spanning many AT-rich intergenic spacers, whereas the Newbler-generated assembly resulted in 18 smaller contigs covering the same ptDNA region. In general, Newbler generated more contigs by breaking at homopolymer runs, even though high-quality reads existed at the ends of all contigs. In only a few instances did the MIRA assemblies differ from those constructed using Newbler and these were all single nucleotide differences in regions of poly-A/T uncertainties in intergenic spacers.

Plastome Structural Characteristics —With the exceptions of Crithmum and Petroselinum, where J_LB and J_LA differ and unique DNA has been incorporated into the LSC region adjacent to the J_LA boundary, all plastomes possess the ancestral angiosperm plastid genome organization. These genomes ranged in size from 152,890 bp (Petroselinum) to 158,355 bp (Crithmum), with size extremes caused by contraction or expansion of the IR at J_LB (Table 2). All plastomes except Crithmum and Petroselinum share identical complements of coding genes, each containing 30 unique tRNA genes, four unique rRNA genes, and 79 unique protein-coding genes. Allowing for duplication of genes in the IR, there are a total of 130 complete, predicted coding regions in these plastomes. Eighteen genes have introns, 16 of which contain a single intron, while two (ycf3, clpP) have two introns. All maintain conserved intron boundaries, and with the exception of the sole group I intron found in trnL-UAA, all are group II introns (Michel et al. 1989). The intron-containing gene rps12 is trans-spliced; its 5′ end exon (5′rpsl2) occurs in the LSC region and its second and third exons (3′rpsl2) reside some 27 kb away in IR_B separated by an intron (an additional copy of 3′rpsl2 with intron exists in IR_A). Similar to other angiosperm plastid genomes, Apiales plastomes are AT-rich.

Table 1.

Comparison of assemblies of five Apiales plastid genomes.

Table 2.

Comparison of major structural features of seven Apiales plastid genomes. Features of Panax schin-seng are provided by Kim and Lee (2004) and those of Daucus carota are provided by Ruhlman et al. (2006).

A comparison of the major structural features of each of these five plastid genomes, along with comparable data for the previously published Daucus and Panax plastomes, is provided (Table 2, Fig. 1). In the Anthriscus cerefolium plastome J_lb occurs within gene rpsl9, resulting in the duplication of a portion of this gene (99 bp) in IR_a at J_LA. There are 2 bp of noncoding sequence between J_la and the 3′ end of gene trnH-GUG in the LSC region. The Tiedemannia filiformis subsp. greenmanii plastid genome contains the same complement of uniquely occurring and duplicated genes as in Anthriscus. In Tiedemannia, J_lb occurs within rpsl9 resulting in the duplication of a portion of this gene (56 bp) in IR_A, and there are 2 bp of noncoding sequence between J_la and trnH-GUG. The Hydrocotyle verticillata plastome contains the same complement of uniquely occurring and duplicated genes as in Anthriscus and Tiedemannia. Here, J_lb extends into rpsl9 resulting in the duplication of a portion of this gene (62 bp) in IR_A, and there are 6 bp of noncoding sequence between J_la and trnH-GUG. Relative to the previous three plastomes, the IR of Petroselinum crispum has contracted ∼1.5 kb such that all of rpsl9 and most of rpl2 are now single-copy. J_lb is contained within rpl2 and a small portion of this gene (37 bp) is duplicated at the end of IR_A. In Petroselinum, there are now only four intron-containing genes duplicated in the IR, not five. Between J_la and the 3′ end of trnH-GUG in the LSC region there are 345 bp of noncoding sequence. Unlike the previous four genomes, the IR of Crithmum maritimum has expanded outwards ∼1.5 kb such that all of rps19 and previously single-copy genes rpl22 and rps3 are now contained within the IR. Nine protein-coding genes are contained entirely within the IR, not six as in the other plastomes, and between J_la and trnH-GUG there are 1,463 bp of noncoding sequence. The circular plastome gene map of Crithmum can be viewed as online supplementary data (supplementary Fig. S1).

Fig. 1.

Comparison of IR-single copy border positions across seven Apiales plastid genomes (vertical lines), with sizes of each of the four major plastid genome components indicated (LSC, IR_B, SSC, IR_A). The various lengths of complete genes (rpsl9, rpl2, ycfl), pseudogenes denoted by asterisks (ycfl^*, rpsl9^*, rpl2^*), and intergenic spacer regions (J_SB-ndhF,J_LA-trnH) adjacent to Jlb, Jsb, Jsa and Jla are also indicated. In Crithmum, Jlb occurs between genes rpl16 and rps3, whereas in Petroselinum, Jlb occurs in the rpl2 5′exon. In all other Apiales examined, Jlb occurs at different positions within rpsl9. In Hydrocotyle and Crithmum, the final 32 and 7 nucleotides, respectively, of Irb adjacent to Jsb are shared by the genes ycfl^* and ndhF, in opposite transcriptional orientations. Crithmum and Petroselinum have incorporated large regions of novel DNA (of 1,463 and 345 bp, respectively) into Jla-trnH.

Junctions J_sb and J_SA are in the same relative gene positions in all seven Apiales plastomes (Fig. 1). Jsb occurs within the 3′ end of ndhF (Hydrocotyle, Crithmum) or between genes ndhF and ycfl^*, the latter existing as a pseudogene in IR_B (Panax, Daucus, Anthriscus, Tiedemannia, Petroselinum). In Hydrocotyle and Crithmum, the final 32 and 7 nucleotides, respectively, of IR_b are shared by genes ycfl^* and ndhF, albeit transcribed in opposite directions (Fig. 1). Ycfl varies in size from 5,382– 5,760 bp, and possesses numerous indels. In all seven plastomes, J_sa occurs in ycfl, but at different positions within this large and variably sized gene. As examples, in the Tiedemannia plastome, 2,069 bp of ycfl are duplicated in the IR, whereas in Hydrocotyle only 768 bp of ycfl are duplicated. Gene rpsl9, at 279 bp in size in all plastomes, overlaps J_lb in Panax, Hydrocotyle, Daucus, Anthriscus, and Tiedemannia, and as a result various lengths (51–99 bp) of rpsl9 pseudogenes are located adjacent to J_la in IR_A (Fig. 1).

An examination of DNA sequences flanking J_sb and J_sa across all seven species reveals that, other than Daucus and Anthriscus (both Apiaceae tribe Scandiceae) whose IR-SSC junctions are in identical positions, there are no significant sequence motifs or repetitive elements in common at or near these junctions which may have instigated IR boundary changes. Alignment of DNA sequences flanking J_la across Anthriscus, Daucus, Hydrocotyle, Tiedemannia, Panax, and Araliaceae species Tetrapanax papyrifer (Wang et al. 2008) and Eleutherococcus senticosus (Yi et al. 2012), also revealed no common sequence on either side of J_la. Differences among these seven species are due to length rather than point mutations, and variable poly(A) tracts are apparent between genes rpl2 and rps19^*, ranging in length from A₈ in Eleutherococcus, Panax, and Tetrapanax, to A₂₁ in Anthriscus. The sequences flanking J_la in Crithmum and Petroselinum plastomes also each appear to be unique, although in these species novel, noncoding DNA is apparent between J_la and 3′trnH-GUG.

The positions and percent identities of gap-free segments of pairwise alignments of Anthriscus, Crithmum, Tiedemannia, and Petroselinum with Daucus are shown as a MultiPip (Fig. 2). The positions of the genes shown on the MultiPip are relative to the Daucus plastome coordinates, the reference sequence. Most of the Daucus sequence aligns with that of the other four species of Apiaceae. Coding regions are represented as mostly unbroken aligning segments of relatively high percent identity (95–98%); noncoding regions, on the other hand, show more mismatches and indels, with percent identity values much lower. Highly variable noncoding regions are scattered throughout the single-copy regions and include intergenic spacers ndhF-rpl32 and rpl32-trnL in the SSC region and 5′trnK-3′rpsl6, 5′rpsl6-trnQ, and trnD-trnY-trnE-trnT in the LSC region. A large segment of the Daucus plastome is not present in any of the secondary sequences. This region occurs in the IR between genes 3′rpsl2 and tmV-GAC at Daucus coordinates ∼99,300–100,800 and has been identified previously as a putative mitochondrion to plastid transfer region (Goremykin et al. 2009; Iorizzo et al. 2012a,b).

Repeat Analysis —Repeat analysis identified 9–29 direct and inverted repeats of 30 bp or longer with a sequence identity of ≥90% among the five newly sequenced plastomes (Table 2). These repeated sequences were detected in coding regions (primarily ycf2 and the three serine transfer-RNA [trnS] genes that recognize different codons, but also in psaA and psaB), introns (ycf3, petB, petD, ndhB, and ndhA), and most commonly, within 21 intergenic spacer regions scattered throughout the entire plastomes. The majority of these repetitive elements were direct repeats. Maximum repeat sizes ranged from 39 bp (Tiedemannia) to 71 bp (Crithmum). Several of these repeats occurred in the same locations and were shared by all five Apiaceae plastomes, either within genes, introns, or intergenic spacers.

Repeat analysis of the 345-bp region between J_la and trnH of Petroselinum failed to detect any direct or inverted repeats ≥30 bp. In contrast, a series of direct repeats characterize the 1,463-bp Crithmum J _la-trnH region, with the four largest of these being perfect or nearly perfect repeats of 71, 67, 49, and 39-bp. Many smaller repeats were also identified through manual examination, with the largest of these being 27- and 21-bp in size, plus an 18-bp repeating unit occurring three times in tandem. All of these repeats were restricted to positions 4-943 of the Crithmum plastome.

The number of homopolymers ≥7 bp in each of the five newly sequenced plastid genomes ranged between 283 and 326, with A or T polymers outnumbering G or C polymers by ratios of 9.0–11.5 to 1 (Table 2). Generally, the largest homopolymers were composed of A′s or T′s, the largest being 21-bp in size (Table 2). In each genome, the majority of homopolymers was 7 bp long (173–186), followed by 8 bp long (42–65), then 9 bp long (28–40). As the length of the homopolymer increased, fewer numbers were identified. Homopolymer sequences of >13 bp occurred infrequently in each genome (0–6). Among the largest homopolymers, the single 21-bp poly-A homopolymer (Anthriscus) occurred between IR genes rpl2 and rps19^*; the single 20-bp poly-A homopolymer (Crithmum) occurred between LSC genes trnT-GGU and psbD; the single 19-bp poly-T homopolymer (Hydrocotyle) occurred between LSC genes rps8 and rpll4; and the single 18-bp poly-A homopolymer (Hydrocotyle) occurred in SSC gene ccsA. The vast majority of homopolymers, and all but one of the largest ones (>13 bp), occurred in intergenic spacers; several others occurred in the variable and large hypothetical coding frames ycfl and ycf2.

Novel DNA —Both Petroselinum and Crithmum contain novel DNA sequences between J_la and LSC gene trnH-GUG of 345 and 1,463 bp, respectively (Fig. 1). In contrast, the other three newly sequenced plastomes are similar to Daucus and Panax in having only 2–6 bp of noncoding sequence in this region. These two large insertions show no sequence similarity to any other portions of their respective plastid genomes (by comparing these insert sequences to their truncated plastid genomes under varying degrees of stringency), nor do they match any ptDNA data currently available in GenBank. A BLAST search querying the 345-bp Petroselinum sequence resulted in multiple hits to several angiosperm mitochondrial genome sequences. The best alignment score, with 35% query coverage (positions 224–345, numbering beginning the first bp in the LSC region following J_la and proceeding counter-clockwise), showed 91% sequence similarity to a mitochondrial DNA (mtDNA) intergenic spacer between genes cytochrome b (cob) and ORF25 in several Daucus carota subsp. sativus cultivars and breeding lines (Bach et al. 2002). The same high match was obtained against the complete Daucus carota subsp. sativus mitochondrial genome sequence, at positions 193,593–193,712 (Iorizzo et al. 2012b) representing an intergenic spacer also bounded by genes cob and atp4 (ORF25). A BLAST search querying the 1,463-bp novel sequence of Crithmum resulted in no significant alignments; however, 39 bp of this sequence immediately adjacent to trnH-GUG showed 87% sequence similarity to a portion of the 345-bp Petroselinum insert and to the Daucus carota mitochondrial genome fragment. The Crithmum trnH-psbA region also contains novel DNA and this region was unalignable with the other Apiales species examined; a BLAST search of this sequence resulted in a single hit to the same spacer region in another accession of Crithmum maritimum (Degtjareva et al. 2009).

Fig. 2.

Multiple percent identity plots (MultiPip) comparing the Daucus plastid genome (as the reference sequence) to Anthriscus, Tiedemannia, Petroselinum, and Crithmum plastid genomes using MultipPipMaker (Schwartz et al. 2003). IR_a is excluded. The top line is a gene map of Daucus (Ruhlman et al. 2006); introns are depicted by white boxes. Sequence similarity of aligned regions in Anthriscus, Tiedemannia, Petroselinum, and Crithmum is shown as horizontal bars in each plot indicating average percent identity between 50–100% (vertical axis). The horizontal axis represents the coordinates in the Daucus plastid genome. Parallel lines, as in psaA and psaB, indicate repeated sequence.

Fig. 2. (cont)

Multiple percent identity plots (MultiPip) comparing the Daucus plastid genome (as the reference sequence) to Anthriscus, Tiedemannia, Petroselinum, and Crithmum plastid genomes using MultipPipMaker (Schwartz et al. 2003). IR_a is excluded. The top line is a gene map of Daucus (Ruhlman et al. 2006); introns are depicted by white boxes. Sequence similarity of aligned regions in Anthriscus, Tiedemannia, Petroselinum, and Crithmum is shown as horizontal bars in each plot indicating average percent identity between 50–100% (vertical axis). The horizontal axis represents the coordinates in the Daucus plastid genome. Parallel lines, as in psaA and psaB, indicate repeated sequence.

Identification of Divergent Regions —A comparison of sequence features across five Apiaceae plastid genomes and 70 noncoding ptDNA loci (>150 bp in Daucus) is presented in Table 3. Aligned lengths for each locus ranged from 172 (trnfM-rpsl4) to 1,545 bp (5′rpsl6-trnQ) and the number of parsimony-informative positions ranged from two (for several loci) to 73 (trnE-trnT). The number of variable alignment positions ranged from seven (trnfM-rpsl4) to 337 (rpl32-trnL). There was no obvious relationship between the length of the aligned sequence and the number of variable or parsimonyinformative positions, a relationship consistent with other plastome sequence comparisons (Shaw et al. 2007). Pairwise sequence divergence values for each locus ranged from 0–3.1% for trnfM-rpsl4 to 8.5–26.2% for rpl32-trnL. The number of indels in each region ranged from two (rpoC2-rpoCl) to 55 (ndhF-rpl32) and the number of parsimony-informative indels ranged from zero (for multiple loci) to 10 (5′rpsl6-trnQ). Percent variability ranged from 6.4 (trnfM-rpsl4) to 38.7% (rpl32-trnL). Of the four combined regions examined, trnD-trnY-trnE-trnT was the most variable, with 30.5% variability, 7.6–19.3% sequence divergence, 321 variable alignment positions, and 51 indels.

Table 3.

Comparison of 74 noncoding loci from five plastid genomes of Apiaceae (Daucus, Anthriscus, Tiedemannia, Petroselinum, and Crithmum). For the individual intergenic spacers and introns, the loci are listed as they occur on the Crithmum ptDNA map (Fig. S1), starting at J_la, proceeding counterclockwise, and excluding the two IR regions. These noncoding loci represent all intergenic spacers >150 bp and introns from both LSC and SSC regions. The trnH-psbA locus excludes Crithmum because of alignment ambiguities.

Continued.

A comparison of the 25 most divergent noncoding plastome regions, as represented by numbers of variable and parsimony-informative alignment positions, is presented in Fig. 3. Among the most variable loci include rpl32-trnL, trnD-trnY-trnE-trnT, trnE-trnT, ndhF-rpl32, 5′rpsl6-trnQ, trnT-psbD, ndhC-3′trnV, rpoB-trnC, petN-psbM, atpH-atpl, 5-trnK-3′rpsl6, petA-psbJ, and psbE-petL. These 13 regions had among the highest percent variability of all loci examined (14.3– 38.7%). They also had among the highest pairwise sequence divergence values (19.3% for trnD-trnY-trnE-trnT, 23.9% for tmE-trnT, and 26.2% for rpl32-trnL). No intron was included in these top 13 most divergent loci, but of those included in the analysis four (ndhA, rpl16, rps16, and atpF) had the greatest numbers of variable alignment positions (88–137) and among the highest sequence divergence and % variability values (Table 3).

Fig. 3.

Histogram showing the numbers of variable (dark boxes) and parsimony-informative (white boxes) positions in 25 of the most variable noncoding loci (intergenic spacers and introns) identified among the 74 regions compared from five Apiaceae plastid genomes.

Discussion

Plastome Comparisons —The plastomes of Anthriscus, Crithmum, Hydrocotyle, Petroselinum, and Tiedemannia are highly conserved in size, structure, and gene order and content compared to the previously published plastomes of Panax (Kim and Lee 2004) and Daucus (Ruhlman et al. 2006). They are also similar to the plastome of Eleutherococcus senticosus (Rupr. & Maxim.) Maxim. (Siberian ginseng), the most recent member of Araliaceae to have its plastid genome sequenced (Yi et al. 2012). Each of these Apiales plastomes shares characteristics of gene content and organization with the ancestral angiosperm plastid genome, as most frequently represented by N. tabacum (Raubeson et al. 2007). The ptDNAs of Anthriscus and Daucus, both Apiaceae tribe Scandiceae, are colinear, although the plastid genome map presented in Ruhlman et al. (2006) for Daucus indicates an additional open reading frame (ORF80) in the IR. Additionally, Daucus contains two copies of trnG-GCC in the LSC region, whereas in all other Apiales plastomes considered herein, one of these is replaced by trnG-UCC (between trnS-GCU and trnR-UCU). The sizes of these Apiales plastomes are within the known range for most angiosperms possessing an IR; however, the size of the Crithmum IR, at 27,993 bp, is at the higher end of its 20–30 kb average size (Plunkett and Downie 2000; Cai et al. 2006; Chumley et al. 2006; Ruhlman et al. 2006; Haberle et al. 2008; Wicke et al. 2011).

The IR boundaries of Anthriscus, Daucus, Hydrocotyle, Tiedemannia, Panax, and Eleutherococcus are also in the same relative positions as in many other unrearranged non-monocot angiosperm plastid genomes. In these Apiales taxa, J_LB occurs within rpsl9, resulting in the duplication of a portion of this gene at the end of IR_a adjacent to J_la. J_la occurs downstream from gene trnH-GUG in the adjacent LSC region. In Petroselinum, however, the IR has contracted 7sim;1.5 kb at the IR_B-LSC boundary, such that J_LB is contained within the rpl2 5′exon, and in Crithmum, the IR has expanded ∼1.5 kb at the IR_B-LSC boundary to duplicate genes rpl22 and rps3. These large IR junction shifts (>1,000 bp) are unusual, for they occur less frequently in angiosperm plastomes than do IR contractions and expansions of <100 bp (Goulding et al. 1996). Another major difference between the plastomes of Petroselinum and Crithmum and those of most other angiosperms is the presence of a large insertion of novel, noncoding DNA in the region between J_LA and the 3′end of trnH-GuG, with portions of these inserts showing high sequence similarity to mitochondrial DNA.

Complete sequence data for 199 angiosperm plastid genomes (representing 134 genera) are publicly available in NCBI′s Organelle Genome Resources database ( http://www.ncbi.nlm.nih.gov/genomes/GenomesGroup.cgi?taxid=2759&opt=plastid; accessed February 1, 2013). Excluding most monocots and other angiosperm species whose plastid genomes are rearranged relative to the ancestral angiosperm plastome gene order (as exemplified by N. tabacum) and those species whose plastomes have lost one copy of the IR (Downie and Palmer 1992; Raubeson and Jansen 2005; Guisinger et al. 2011; Jansen and Ruhlman 2012), the Jla boundaries for the remaining 108 plastomes (73 genera) typically occur between rpl2 and trnH-GuG, with the latter present in its entirety only in the LSC region (Wang et al. 2008). While the N. tabacum plastome is frequently inferred to represent the ancestral plastid genome organization for angiosperms, it may not actually contain the exact ancestral IR boundaries (Raubeson et al. 2007). In many instances (like those of some Apiales), an rpsl9 pseudogene of varying length remains in IR_A immediately adjacent to the LSC region; in three other instances, rpsl9 is full-length (and therefore exists as two complete copies in the IR) and J_la occurs between rpsl9 and trnH-GuG. The majority of these plastomes (94) have between 0 and 30 nucleotides of noncoding sequence between J_la and the 3′end of trnH-GuG. Twelve of the remaining 14 plastomes have between 34 and 71 bp of noncoding sequence in this region, one species of Gossypium L. has 121 bp of noncoding sequence, and Hevea brasiliensis (Willd. ex A. Juss.) Müll. Arg. has 198 bp of noncoding sequence in this region. Therefore, the large insertions of novel, noncoding sequence between J_la and trnH-GuG of 345 nucleotides (Petroselinum) and 1,463 nucleotides (Crithmum) are highly unusual among angiosperms. Similar large insertions of novel DNA (of between 214 and 393 nucleotides) have been reported within this same region in other Apiaceae species and their pattern of distribution suggests that these insertions may be restricted to the “apioid superclade” of Apiaceae subfamily Apioideae (Peery et al. 2007). In most monocot plastomes, the gene trnH is located between rpl2 and rpsl9 in the IR (Mardanov et al. 2008; Wang et al. 2008; Wang and Messing 2011), and between J_la and the first gene in the LSC region, psbA, there are 82–148 nucleotides of noncoding sequence. These sequences are also unique within the plastome (Hansen et al. 2007), but to the best of our knowledge these regions do not show any similarity to angiosperm mitochondrial DNA.

Mitochondrial DNA Transfer into the Plastid Genome — Angiosperm mitochondrial genomes readily accept DNA of plastid and nuclear origin (Knoop 2004; Kubo and Mikami 2007; Richardson and Palmer 2007; Goremykin et al. 2009; Iorizzo et al. 2012b), as well as sequences of horizontal origin from foreign genomes (Bergthorsson et al. 2003, 2004; Mower et al. 2010). In contrast, intracellular or horizontal gene transfer to the plastid genome is extremely rare, with the few reported cases restricted to some red and green algal plastids (Sheveleva and Hallick 2004; Rice and Palmer 2006; Kleine et al. 2009). Until recently, no convincing evidence of DNA transfer from either the nucleus or the mitochondrion into the angiosperm plastome has been documented (Palmer 1985b, 1990; Rice and Palmer 2006; Richardson and Palmer 2007; Bock and Timmis 2008; Bock 2010); in fact, such a transfer has generally been considered highly unlikely, possibly reflecting the lack of a DNA uptake system in plastids (Richardson and Palmer 2007; Kleine et al. 2009; Smith 2011).

Several recent papers, however, have identified DNA sequences of putative mitochondrial origin in angiosperm plastomes or invoked intracellular or horizontal transfer to explain the origin of novel sequences within plastid genomes. In the Pelargonium × hortorum L. H. Bailey (Geraniaceae) plastome, for example, ORF56 and ORF42 in the trnA-UGC intron showed high sequence similarities to the mitochondrial ACR-toxin sensitivity (ACRS) gene of Citrus jambhiri Lush. (Rutaceae) and to the mitochondrial pvs-trnA gene of Phaseolus L. (Fabaceae), respectively (Chumley et al. 2006). ORF56 is also shared by both plastid and mitochondrial genomes in other taxa, suggesting transfer from the mitochondrion into the plastid (Chumley et al. 2006). The plastid genome of Trifolium subterraneum L. (Fabaceae) contains 19.5 kb of unique DNA distributed among 160 fragments ranging in size from 30–494 bp, some of which was hypothesized to represent instances of horizontal transfer from bacterial genomes (Cai et al. 2008). Other legume plastid genomes contain much repetitive sequences and unique DNA too, and speculation on their origins included intracellular transfers from the mitochondrion or nucleus, or horizontal transfers from other genomes, possibly pathogenic bacteria (Cai et al. 2008). Another example of transfer of mitochondrial DNA to the plastid genome has been identified in Apocynaceae tribe Asclepiadeae (Ku et al. 2013; Straub et al. 2013). In this case, it was suggested that homologous recombination is responsible for the transfer of the mitochondrial sequence into the plastome. Thus, it appears the transfer of mitochondrial DNA to the plastome is more common than previously thought.

Goremykin et al. (2009) uncovered two small sequences (of 74 and 126 bp) of the Vitis vinifera L. (Vitaceae) mitochondrial genome that were highly similar to regions of the Daucus carota plastome. BLAST analysis of the larger sequence revealed a high similarity to the coding region of the mitochondrial cytochrome c oxidase subunit 1 gene (cox1l), prompting the authors to suggest that its presence in the Daucus plastome might possibly represent a rare transfer of DNA from the mitochondrion into the plastid. These two sequences are contained within a large 1,439-bp fragment of the D. carota IR (positions 99,309–100,747 and 139,407–140,845 in Ruhlman et al. 2006) that is a part of the 3′rpsl2-trnV-GAC intergenic spacer region; this fragment, however, has no similarity to any other published plastid nucleotide region (Goremykin et al. 2009). The MultiPip (Fig. 2) shows that this large fragment of the Daucus plastome is not present in any other Apiaceae plastome sequenced herein, including its closest ally Anthriscus (both Apiaceae tribe Scandiceae). Subsequently, Iorizzo et al. (2012a,b) confirmed the presence of this unique 1,439-bp fragment (identified as 1,452 bp in their study) in the Daucus carota plastome, while also discovering that it is present as three noncontiguous sequences in the D. carota mitochondrial genome. Further, the authors documented consistent conservation of a large portion of this region across all mitochondrial genomes of the diverse Apiaceae species they examined (representing seven additional species of Daucus and six other genera of Apiaceae, the latter including five species from the “apioid superclade” and Eryngium planum L. of subfamily Saniculoideae). In the plastid genome, this fragment, or a large portion thereof, was present only in Daucus and Cuminum L., both of tribe Scandiceae subtribe Daucinae (Lee et al. 2001). Iorizzo et al. (2012a) concluded that their results provided strong evidence of a mitochondrial to plastid transfer of DNA, and the presence of this putative mitochondrial to plastid transfer region in Scandiceae subtribe Daucinae but not in subtribe Scandicinae (Anthriscus) suggests this transfer event occurred sometime during the early evolution of the Daucinae clade. The authors also suggested that a retrotransposon-like event might have been responsible for transferring this sequence into the plastid genome (Iorizzo et al. 2012b).

Crithmum and Petroselinum also incorporate novel DNA into their plastomes. BLAST searches of the 345-bp Petroselinum unique insert resulted in 122 bp of this region having multiple hits to angiosperm mitochondrial DNAs from a variety of angiosperms, with the best alignment score matching a portion of the cob-atp4(ORF25) intergenic spacer in the Daucus carota mitochondrial genome. This high sequence similarity is suggestive of a possible intracellular transfer of DNA into the plastid from the Petroselinum mitochondrial genome. BLAST searches of the 1,463-bp Crithmum insert resulted in no significant hits, except for a small (39-bp) region that also matched a portion of the Petroselinum insert and the D. carota mitochondrial genome fragment. The origin of the large number of repetitive elements within the Crithmum insert remains unclear. We interpret these results cautiously, for we have not firmly demonstrated the mitochondrial provenance of this novel DNA, and suggest that the sequencing of whole mitochondrial genomes from the Apiaceae, specifically from Petroselinum, may eventually shed light on its origin.

Repeated Sequences —For the five newly sequenced plastomes, repeat analysis identified 9–29 dispersed repeats of 30 bp or longer with a sequence identity of 90% or greater. In the previously published Daucus plastome, repeat analysis identified 12 direct and two inverted repeats of ≥30 bp, four of which were direct repeats (of 30–70 bp) occurring in ycf2 (Ruhlman et al. 2006). In Panax, nine direct and two inverted repeats ≥30 bp were identified, the largest of these being four tandem repeats of 57-bp in ycf1 (Kim and Lee 2004). In the Eleutherococcus plastome, a total of 23 repeats ≥30 bp were located, the largest being 79 bp in size and also occurring in ycf2 (Yi et al. 2012). Repeated sequences, other than the IR, are considered to be uncommon in plastid genomes and are more prevalent in those plastomes that have undergone major changes in genome organization (Palmer 1985a, 1991; Chumley et al. 2006; Haberle et al. 2008; Weng et al. 2014). Considering all Apiales plastomes investigated to date, repeats occurred most often in five introns and numerous intergenic spacer regions scattered throughout the genome; within coding regions, repeats were prevalent in genes ycf2, trnS (3 genes), psaA, and psaB. The repeated elements identified in Apiales are similar in number, size, and location to those of other angiosperms whose plastomes are structurally unrearranged (Ruhlman et al. 2006; Yi et al. 2012). The origin of these repeats is not known, although replication slippage could be responsible for generating direct repeats (Palmer 1991).

Similarly, homopolymer sequences of >13 bp occurred infrequently in each plastid genome. In Panax, 18 homopolymers ≥10 bp were identified previously (the largest being 13-bp in size), the majority being composed of multiple A's or T's (Kim and Lee 2004). In Eleutherococcus, 27 homopolymers ≥10 bp were identified, of which 22 were composed of multiple A or T bases and the largest were also of 13 bp in size (Yi et al. 2012). These repetitive sequences, of which there are many in the Apiales, may be potentially useful for population-level genetic studies, as sequence length polymorphisms can be surveyed among multiple accessions of a species revealing haplotype variation (Yi et al. 2012, and references cited therein). Homopolymer runs are difficult to resolve using the Roche/454 platform (and other sequencing methods) and contribute to sequencing errors (Moore et al. 2006). In our study, Sanger sequencing was used for independent confirmation of only a few homopolymeric repeats at or near IR-single copy junctions. However, with the low depth of coverage at some contig ends and the fact that very few homopolymer runs were confirmed by Sanger sequencing, we caution that sequencing errors are likely present, especially within the largest homopolymeric regions.

All major shifts in J_LB position documented in Apiaceae are restricted to members of the “apioid superclade” (Plunkett and Downie 1999, 2000). These include both major IR expansions and contractions, such as those presented herein for Crithmum and Petroselinum. The prevalence of IR junction shifts in this one (albeit large) clade suggests that some common mechanism is responsible for flux in position of J_LB and that this mechanism likely originated in the immediate common ancestor of the group (Plunkett and Downie 2000). The underlying mechanisms responsible for causing such large IR junction shifts are not yet known, but may involve homologous recombination between small, dispersed repeated sequences, adjacency to tRNA genes, or double-stranded breakage followed by recombination between poly(A) tracts (Palmer 1985a; Goulding et al. 1996; Plunkett and Downie 2000; Wang et al. 2008). Small IR expansions (<100 bp), such as those found in many Nicotiana species, are likely caused by gene conversion (Goulding et al. 1996; Wang et al. 2008). Comparisons of DNA sequences flanking IR-single copy regions in Apiales reveal that there are no dispersed repeats at or near these junctions, although this doesn't preclude such repetitive elements from instigating IR boundary changes in other apioid taxa. However, the absence of any significant sequence motifs or repetitive elements at or near IR-single copy junctions in other angiosperms exhibiting IR boundary fluxes suggests that these types of sequences may not be the typical cause of IR spreading (Palmer 1985b; Goulding et al. 1996), although Wang et al. (2008) observed that IR-LSC junctions are indeed found at either poly(A) tracts or in A-rich regions. In the “apioid superclade,” the presence of both large-scale expansions and contractions of the IR at J_LB and large, novel insertions in the LSC region adjacent to J_LA may be coincidental, or the latter may be mediating, in some way, these IR boundary changes, or vice versa.

Expansion/Contraction of the IR as a Phylogenetic Marker— Because of their infrequent occurrence, major structural rearrangements of the plastid genome usually can provide strong evidence of monophyly (Downie and Palmer 1992; Raubeson and Jansen 2005; but see Downie et al. 1991, for an exception). These rare, large-scale mutational changes provide complementary markers that can be used alongside nucleotide substitutions in molecular systematic studies. Variation in overall size of the IR is common in angiosperms, but such a character cannot readily be used in a phylogenetic analysis because length mutations can occur anywhere within the IR, making homologous size variants difficult to assess. Specific expansions and contractions of border positions of the IR, however, can be used to demarcate monophyletic groups, and these markers have been used previously to delimit major clades in Ranunculaceae (Johansson and Jansen 1993; Hoot and Palmer 1994; Johansson 1998), Berberidaceae (Kim and Jansen 1994), Nicotiana (Goulding et al. 1996), Campanulales (Knox and Palmer 1999), Poaceae (Guisinger et al. 2010), and monocots (Wang et al. 2008). Minor changes in IR junction positions (<100 bp) are more common (Goulding et al. 1996), but this increased frequency suggests that they would not make very reliable phylogenetic markers.

In Apiaceae, the IR has both expanded and contracted significantly, with successive contractions identified within the “apioid superclade,” five of which are parsimonyinformative (Plunkett and Downie 2000). Members of Apiaceae tribes Careae and Pyramidoptereae, both previously attributable to the “Aegopodium group,” share a large expansion of J_LB relative to the Nicotiana plastome and are monophyletic sister groups based on a nuclear ribosomal DNA internal transcribed spacer (nrDNA ITS) sequence phylogeny (Plunkett and Downie 2000; Banasiak et al. 2013). Members of Apiaceae tribe Apieae (the “Apium group”), as represented herein by Petroselinum, but also by Ammi L., Anethum L., Apium L., Foeniculum Hill., and Ridolfia Moris, in Plunkett and Downie (2000), are all characterized by a ~1.5 kb J_LB contraction. Conium L. (Conium clade) and Pimpinella L. (tribe Pimpinelleae) also possess a similar-sized IR contraction as that found in tribe Apieae (Plunkett and Downie 2000), but whether this contraction serves as a synapomorphy uniting these three major clades of taxa remains to be seen. In the absence of sequence data flanking their IR junctions, ambiguity remains in assessing the precise location of J_LB.

To date, all major IR expansions and contractions uncovered in the Apiales are restricted to the “apioid superclade,” a large, distally branching group within Apiaceae subfamily Apioideae whose higher-level relationships are unclear (Plunkett and Downie 2000). Resolution of relationships among its constituent 14 tribes and major clades (Downie et al. 2010), as well as the formal recognition of these major clades, must await supporting data, such as that provided by the distribution of well-characterized, coincident IR junction locations. Such rearrangements are rare enough that they have the potential to resolve with confidence a particular branching point in a phylogeny and define monophyletic groups (Downie and Palmer 1992). We are currently surveying for, characterizing, and circumscribing the distribution of major IR junction shifts in the “apioid superclade” to more rigorously evaluate the extent of IR junction flux in the group (R. Peery et al. unpubl. data). Even though these rearrangements alone are unlikely to provide a comprehensive framework of relationships, when used in conjunction with other molecular and traditional approaches, they have the power to illuminate phylogeny or to provide additional support for otherwise weakly-supported clades. Of course, there is also the possibility that IR fluxes are completely random, as they appear to be among closely related Nicotiana species (Goulding et al. 1996), or successive expansion-contraction events have occurred, as has been suggested in adzuki bean (Perry et al. 2002) and Pelargonium × hortorum (Chumley et al. 2006), which would confound interpretation of phylogeny.

Noncoding Regions as a Source of Phylogenetic Information in Apiaceae— At present, nrDNA ITS sequences comprise the most comprehensive database for Apiaceae phylogenetic study (Downie et al. 2010). These data have often been supplemented with sequences from a variety of ptDNA introns and intergenic spacers, yet for those taxa where radiations have been relatively recent it has been difficult to generate sufficient phylogenetic signal because of the relatively slow rate of evolution of the various plastid loci being compared. Within the “apioid superclade,” phylogenetic resolution is particularly poor within several groups of genera, such as the clade of perennial, endemic Apiaceae subfamily Apioideae of western North America (Sun et al. 2004; Sun and Downie 2010), Bunium (Degtjareva et al. 2009), Conium (Cordes 2009), Pimpinella (Magee et al. 2010), the Arracada clade (Danderson 2011), Heracleum (Yu et al. 2011), and Angelica (Liao et al. 2013). Each of these groups contains unresolved terminal clades of presumably closely related species, which would benefit considerably by the inclusion of additional phylogenetically informative characters.

The development of genetic markers to infer evolutionary relationships among recently diverged taxa remains an important challenge in molecular systematic studies and, to this end, considerable advances have been made in the identification and development of markers from the angiosperm plastid genome (Shaw et al. 2005,2007). Despite these efforts, a plastid region identified as highly variable in one group may not actually be phylogenetically informative in another. As an example, the Crithmum trnH-psbA intergenic spacer region shows extensive sequence divergence when compared with the other four Apiaceae species and could not be aligned unambiguously with these sequences. This region is substantially shorter in Crithmum ptDNA (127 bp) than that of the other species (180–191 bp); it also contains novel DNA, many small, repeated sequences in direct orientation, and is AT-rich (73%). This spacer region is highly variable in angiosperms (Aldrich et al. 1988; Shaw et al. 2007) and has been proposed for DNA barcoding (Kress et al. 2005), yet in our study this locus ranked 59th out of 74 in its number of variable alignment positions. Furthermore, in Heracleum and its allies (Apiaceae tribe Tordylieae), the locus is quite conserved, with one small internal region characterized by an inversion and duplication (Logacheva et al. 2008). Repeats, large deletions, and small inversions in trnH-psbA have been reported for other species of Apiales, severely confounding the use of this locus in phylogenetic analyses (Degtjareva et al. 2012).

Practical information to emerge from our study is the identification of divergent plastid regions for enhancing resolution of low-level phylogenetic relationships in Apiaceae. A comparison of levels of DNA sequence divergence across five Apiaceae plastomes, representing species from four disparate lineages within the family (tribes Apieae, Oenantheae, Pyramidoptereae, and Scandiceae), revealed that intergenic spacer regions rpl32-trnL, trnE-trnT, ndhF-rpl32, 5′rpsl6-trnQ, and trnT-psbD are among the fastest-evolving loci, as measured by numbers of variable and parsimony-informative characters, and also by total indels, pairwise sequence divergence estimates, and percent variability. Aligned lengths of each of these regions ranged from 850–1,545 positions, and these sizes are of sufficient length for ease of PCR-amplification and DNA sequencing. These five spacers appear to be good candidates for resolving recent divergences. Consideration of two spacers simultaneously because of gene adjacency and co-amplification as a single contiguous unit (i.e. ndhF-rpl32 and rpl32-trnL, and trnE-trnT and trnT-psbD) offers even better regions for phylogenetic analysis. The trnD-trnY-tmE-trnT combined region is substantially longer and presents the greatest number of parsimony-informative characters overall. Shaw et al. (2007) identified nine noncoding regions as the best possible choices for low-level phylogenetic studies of angiosperms. All nine of these regions are within our top 13 most divergent loci. To these nine, we also include tmE-trnT, rpoB-trnC, petN(ycf6)-psbM, and the trnD-trnY-trnE-trnT combined region; these regions, too, were also noted as highly variable or in the “Tier 1” category of Shaw et al. (2005, 2007). Although we acknowledge that a plastid region identified as highly variable in one group may not be as useful phylogenetically in another group, it is interesting that the Shaw et al. (2005, 2007) studies identified the same highly variable loci for several disparate angiosperm lineages as we did for the Apiaceae. The resolution of recent divergences in Apiaceae would benefit considerably by the inclusion of any or all of these highly variable loci.

No intron was included in the top 13 most divergent noncoding loci, although four (ndhA, rpl16, rps16 and atpF) made the top 25. Of these, the rpl16 and rps16 introns have been used extensively in Apiaceae phylogenetic studies, whereas we are unaware of any studies of the family incorporating ndhA or atpF intron sequences. Group II introns of the plastid genomes of land plants, which includes the aforementioned four introns, show a strong relationship between the functional importance of its secondary structural features and the likelihood of mutational change, with those domains and subdomains essential for intron-associated functions most conserved evolutionarily (Downie et al. 1996, 2000; Kelchner 2002). Intergenic spacers are less functionally constrained than introns and therefore more likely to vary, thereby providing more variable and phylogenetically informative characters for phylogenetic analyses. The greater variability of intergenic spacer regions over introns in plastome comparisons has been reported previously (Shaw et al. 2007). In Apiaceae, intron sequences do not provide enough informative characters to resolve relationships below generic levels (Downie et al. 1996, 2000, 2008; Downie and Katz-Downie 1999).

In conclusion, a comparative analysis of seven complete plastid genomes from the Apiales, five of which were sequenced during the course of this investigation, has provided a plethora of information on their structural organization and sequence evolution. The features of five of these genomes are consistent structurally with typical plastomes of other non-monocot angiosperm species, while two, Crithmum and Petroselinum, each demonstrate a ~1.5 kb expansion or contraction of its IR. These two species also incorporate novel DNA adjacent to J_LA, indicative of a mitochondrial to plastid transfer of DNA. A comparison of 74 loci, representing all noncoding regions >150 bp from throughout the single-copy portions of five Apiaceae plastomes, demonstrated a wide range of sequence divergence in different regions. Several regions were identified that yielded greater numbers of variable and parsimony-informative characters than many of those regions heretofore commonly employed in Apiaceae molecular systematic studies, providing new information on the potential phylogenetic utility of plastid noncoding loci. Therefore, these highly variable loci should be more useful in future interspecific phylogenetic, intraspecific phylogeographic, and population-level genetic studies of these plants. The distribution of specific IR junction shifts, on the other hand, has the potential to demarcate major clades within the family, particularly those within the “apioid superclade” whose higher-level relationships remain elusive.

The frequency and large size of IR junction shifts in Apiaceae plastomes appear to be unusual in angiosperm families whose plastid genomes are otherwise unrearranged, indicating that the group represents a model system in which to study the mechanisms leading to large-scale expansions and contractions of the IR. Concomitant with these IR shifts is the incorporation of novel DNA in the LSC region adjacent to J_LA, some of which represents transfer of DNA from the mitochondrial genome into the plastid. With other reports of a putative mitochondrial to plastid transfer region elsewhere in the plastomes of Daucus and Cuminum (Goremykin et al. 2009; Iorizzo et al. 2012a) it will be well worthwhile to examine additional Apiaceae plastomes for instances of intracellular transfer of DNA. To improve understanding of the mechanisms causing IR junction shifts and to further characterize and circumscribe the putative mitochondrial DNA insert, the sequences spanning the four IR-single-copy junctions are being compared in additional species of Apiaceae (R. Peery et al. unpubl. data). By combining these data and analyzing them within the context of a ptDNA-derived phylogeny of the group, it will be possible to determine when the first IR junction shift occurred and the extent of IR junction mobility during evolution of the group. Complete plastome sequences of other members of the “apioid superclade” will further advance understanding of the evolutionary events in the plastome that accompanied umbellifer diversification.