1. Introduction

The genus Panax L. (Araliaceae) is one of the most medicinally important plant in East Asia. It includes seven well-recognized species and one species complex that is disjunctly distributed in East Asia and eastern North America (Wen and Zimmer, 1996; Wen, 1999; Lee and Wen, 2004). Almost every species within this genus has been used as a medicinal herb in East Asia, especially in China (Yang et al., 1988). Because of its considerable medicinal benefits, Panax has been involved in many molecular-based phylogenetic analyses in the past few decades. However, high-resolution and well-supported phylogeny within this genus remain elusive. Phylogenetic analysis relying on ITS, 18S rRNA, and plastid fragments have shown that the genus Panax is a monophyletic group; however, the same approaches have not resolved infra-generic relationships within the clade (Wen and Zimmer, 1996; Wen, 1999; Zhu et al., 2003; Lee and Wen, 2004). Although recent studies have clarified the phylogenetic relationships between Panax trifolius, Panax stipuleanatus, Panax pseudoginseng, Panax notoginseng, and the Panax bipinnatifidus species complex, the phylogenetic relationships between Panax ginseng, Panax japonicus, and Panax quinquefolius remain unresolved (Shi et al., 2015; Zuo et al., 2017).

The plastomes of angiosperms are typically circular DNA, consisting of two copies of a large inverted repeat (IR) region separated by a large single-copy (LSC) region and a small single-copy (SSC) region (Raubeson and Jansen, 2005; Wicke et al., 2011). Because plastome sequences show a high level of divergence between species and even populations, these sequences provide valuable information for resolving complex relationships in plants (Moore et al., 2007, 2010; Jansen et al., 2007; Parks et al., 2009). With the advent of second-generation DNA sequencing technologies, plastomes have been widely used in recent years to reconstruct robust phylogenies and to identify species (Nock et al., 2011; Yang et al., 2013; Ruhsam et al., 2015; Huang et al., 2016). To date, a total of 27 complete Panax plastomes (Supplementary Table S1) have been sequenced (e.g., Kim and Lee, 2004; Choi et al., 2014; Zhao et al., 2015; Han et al., 2016; Zhang et al., 2016; Nguyen et al., 2017). These plastid genomic resources provide an opportunity to reconstruct a well-supported phylogeny of this medicinally important genus.

The species P. stipuleanatus (Fig. 1) is restricted to the montane evergreen broad-leaved forests along the border between Vietnam and China in southeast Yunnan province (Xiang and Lowry, 2007). In this region, P. stipuleanatus has been traditionally used as a substitute for P. notoginseng, and its roots are extensively collected and sold in the local markets (Yang et al., 1988). Because of overharvesting, natural populations of Panax stipuleanatas have been markedly reduced. As a result, the species was listed as an endangered species by the Ministry of Environmental Protection of P. R. China in 2013 (http://www.zhb.gov.cn/gkml/hbb/bgg/201309/t20130912_260061). Available genomic resources for P. stipuleanatus are therefore limited, and little is known about the plastome features of the species.

Here, we present the complete plastome of P. stipuleanatus obtained through Illumina sequencing and a reference-guided assembly of de novo contigs. In conjunction with the previously published plastomes, the features of the plastomes among Panax species were compared, including gene content and order. Highly divergent DNA regions offering potential use in further species identification and phylogenetic analysis were also identified. Phylogenetic relationships among Panax species were analyzed with the plastome-based dataset. These results may well broaden our understanding of the evolutionary history of Panax.

2. Materials and methods 2.1. Genome sequencing and assembly

Total genomic DNA was isolated from the silica-gel-dried leaves of P. stipuleanatus, collected in Maguan County, Yunnan Province, China, using a modified CTAB method (Doyle and Doyle, 1987; Yang et al., 2014). Voucher specimens (JYH-2016, 466) were deposited in the Herbarium of the Kunming Institute of Botany, CAS (KUN). Genomic DNA was randomly fragmented into 400e600 bp with an ultrasonicator. Short-insert (500 bp) paired-end libraries were constructed using the Genomic DNA Sample Prep Kit (Illumina), according to the manufacturer's protocol, and then sequenced on the Illumina HiSeq 2500 system at BGI (Shenzhen, Guangdong, China). The Illumina raw data were filtered with a NGS QC Toolkit (Patel and Jain, 2012). High-quality reads were assembled into contigs using SPAdes v3.10.1 (Nurk et al., 2013) with its default parameters. The representative plastome sequence contigs were then mapped onto the reference plastome sequence of P. japonicus (Genbank accession number: KX247146) in Bowtie v2.2.6 (Langmead and Salzberg, 2012) with its default-preset options. Assemblies were then assessed and connected using Bandage (Wick et al., 2015). The validated complete plastome sequences were deposited in GenBank (Table S1).

2.2. Genomic annotation and comparison

The plastome of P. stipuleanatus was annotated with the online software tool, DOGMA (Wyman et al., 2004), coupled with manual corrections for start and stop codons. All tRNAs were further confirmed by the tRNAscan-SE v1.21 (Schattner et al., 2005). Functional classification of the plastid genes was determined by referring to the online database CpBase (http://chloroplast.ocean.washington.edu/). A circular map of the plastome was drawn in OrganellarGenomeDRAW (Lohse et al., 2007).

Complete plastomes of the P. bipinnatifidus species complex, P. ginseng, P. japonicus, P. notoginseng, Panax quiquefolius, P. trifolius, and Panax vietnamensis were downloaded from NCBI. Multiplesequence alignments were performed with MAFFT software (Katoh et al., 2002), and manually edited where necessary. Geneious v7.0 (Kearse et al., 2012) was used to compare the boundaries of the LSC, IR, and SSC regions among the Panax plastomes. To compare sequence divergence among different Panax plastomes, the mVISTA tool was used (Frazer et al., 2004), with P. ginseng set as the reference. Single nucleotide polymorphisms (SNPs) occurring across Panax plastomes were identified with the Shuffle-LAGAN model in Geneious v7.0 (Kearse et al., 2012). Divergent percentages of SNPs among the homologous regions across these species were also calulated.

2.3. Synonymous (Ks) and non-synonymous (Ka) substitution rate analysis

Non-synonymous (Ka) and synonymous (Ks) substitutions, and their ratios (Ka/Ks), are important indictors that reflect the plastome evolution and natural selection (Yang and Nielsen, 2000). Protein-coding exons were extracted from the plastomes of the eight Panax taxa and Aralia undulata Handel-Mazzetti. These genes were translated into protein sequences and aligned separately using Geneious v7.0 (Kearse et al., 2012). The Ks and Ka substitution rates for each protein-coding exon were estimated in DnaSP v5.0 software (Librado and Rozas, 2009).

2.4. Phylogenetic analysis

All 28 available plastome sequences of Panax were included in the analysis. A. undulata was used as the outgroup to root the phylogenetic tree. Maximum-likelihood (ML) analyses were conducted based on the following data partitions: (1) the whole plastome; (2) the protein-coding exons (Table S2); (3) the LSC regions; (4) the SSC regions; (5) the IR regions; and (6) the introns and intergenic spacers. All sequences were aligned with the software MAFFT (Katoh et al., 2002). All gaps in the sequence alignment were excluded. For each dataset, the best-fitting partition scheme and nucleotide substitution models were screened in the program PartitionFinder v2.1.1 (Lanfear et al., 2012). For each analysis, the branch lengths were linked, and the models of nucleotide substitution were restricted to RAxML (Stamatakis, 2006); the "greedy" search algorithm was selected. All ML analyses were done in RAxML-HPC BlackBox v8.1.24 software (Stamatakis, 2006). The bootstrap (BS) value for each branch was computed with 1000 bootstrap replicates.

3. Results 3.1. Plastome features

In total, we obtained 4, 948, 544 paired-end clean reads for P. stipuleanatus. These reads were used to assemble the P. stipuleanatus plastome. The overall size of the P. stipuleanatus plastome is 156, 069 base pairs (bp), and it shows a typical quadripartite structure, including a pair of IRs (each 25, 887 bp) that divide the genome into LSC (86, 126 bp) and SSC (18, 169 bp) regions (Fig. 2, Table 1). The total length of the coding regions (i.e., proteincoding genes, tRNA genes, and rRNA genes) and non-coding regions (i.e., introns and intergenic spacers) are 91, 632 bp and 64, 437 bp, respectively (Table 1). The plastome of P. stipuleanatus possesses 114 unique genes, including 80 protein-coding genes, 30 tRNAs, and four rRNA genes. Of these, 12 protein-coding genes (atpF, ndhA, ndhB, petB, petD, rpl16, rpl2, rpoC1, rps12, rps16, clpP, and ycf3), and six tRNAs (trnA-UGC, trnG-GCC, trnI-GAU, trnK-UUU, trnL-UAA, and trnV-UAC) contain at least one intron (Table 2).

3.2. IRs expansion and contraction

The IRs/LSC and IRs/SSC boundaries among Panax plastomes were compared (Fig. 3). The extensions of the IRs into rps19 and the intergenic spacer between rpl2 and trnH-GUG occur respectively at the IRa/LSC and IRb/LSC boundaries. Although the expansion of the IRs into the ycf1 pseudogene at the IRs/SSC junctions occurred in all species, the overlap between the ycf1 pseudogene and ndhF was only detected in P. quiquefolius and P. vietnamensis.

3.3. Sequence divergence in the Panax plastomes

Plastome sequence divergences were identified among the homologous regions across the eight Panax taxa (Fig. 4). We identified 1130 SNPs in the matrix of these plastomes, with average variant frequency of 0.72%. SNP mutations include 743 SNP sites detected in the LSC region, 291 SNPs in the SSC region, and 48 SNPs in the IR region. Their corresponding average variant frequencies were 0.86%, 1.60%, and 0.19%. In addition, 498 SNPs (average variant frequency = 0.62%) were detected in the protein-coding exons, while 620 SNPs (average variant frequency = 0.96%) were detected in the non-coding regions (Table 3). The divergent frequencies of the coding regions ranged from 0.07% to 2.41% (Supplementary Table S3), and those of the non-coding regions ranged from 0.11% up to 5.66% (Supplementary Table S4). According to the sequence divergences, we screened nine protein coding regions-ccsA, matK, ndhF, petL, psaI, rpl22, rpoA, rps3, and ycf1-which showed the percentage of SNPs >1%. We also scanned 10 non-coding regions-atpA-atpF, ccsA-ndhD, infA-rps8, ndhI-ndhA, psbK-psbI, rpl14- rpl16, rpl2-trnH-GUG, rpl22-rps19, rps19-rpl2, and trnY-GUA-trnEUUC-which showed a divergence proportion >3%. These particular plastid DNA regions (19 in all) may be utilized as potential molecular markers to reconstruct the phylogeny of Panax and to identify different species of this genus.

3.4. Synonymous (Ks) and non-synonymous (Ka) substitution rate

The Ks values of the eight Panax plastomes ranged from 0.0166 to 0.0218, while the Ka values ranged from 0.0032 to 0.0034; the Ka/Ks ratio ranged from 0.1391 to 0.1729 (Table 4, Fig. 5). We identified five genes (i.e., cemA, matK, ndhA, ndhG, and ycf2) with Ka/Ks values greater than one. Three of them (ndhA, ndhG, and ycf2) were identified in P. quinquefolius; one gene (ycf2) was identified in P. japonicus, P. vietnamensis and P. notoginseng; two genes were identified in P. ginseng (ndhG and ycf2), P. stipuleanatus (ndhA and matK), the P. bipinnatifidus species complex (ndhA and matK) and P. trifolius (cemA and ycf2). However, their P-values were greater than 0.05 (Supplementary Table S5), which suggested that no single protein-coding gene in the eight Panax plastomes has yet to be positively selected in a statistically significant way (Yang and Nielsen, 2000).

3.5. Phylogenomic analysis

Six data partitions from the 29 plastomes were used to perform phylogenetic reconstruction. The tree topologies for the whole plastome (Fig. 6A), protein-coding exons (Fig. 6B), and SSC (Fig. 6D) were congruent with each other, only differing with regards to support values at the interior nodes. By comparison, the phylogenetic relationships based on the LSC regions (Fig. 6C), IR (Fig. 6E), and introns and intergenic spacers (Fig. 6F) showed many similarities with the results from the other three datasets, except for the positions of P. notoginseng, which clustered into a new branch with P. ginseng and P. quinquefolius. Among the six data partitions, the phylogenetic tree for the protein-coding regions received the highest support for the major branches. As the tree shows (Fig. 6B), the eight Panax species were fully resolved as four, well-supported monophyletic clades. The clade represented only by P. trifolius was placed at the basally branching position. The second clade (BS = 100%) includes P. stipuleanatus and the P. bipinnatifidus species complex. The third clade (BS = 100%) contains P. ginseng and P. quinquefolius. The fourth clade (BS = 100%) consists of P. japonicus, P. vietnamensis, and P. notoginseng, in which the subclade of P. japonicus + P. vietnamensis (BS = 100%) diverged from P. notoginseng. In addition, those taxa (P. ginseng, P. notoginseng, P. quinquefolius, and P. vietnamensis) with more than one plastome currently available, were recovered as well-supported monophyletic lineages in all the datasets.

4. Discussion 4.1. Comparison of the plastomes in Panax

We compared the plastome of P. stipuleanatus with that of seven published Panax species. All Panax plastomes shared 114 unique genes (80 protein-coding genes, 30 tRNAs, and four rRNAs) in the same order. Although several previous studies (e.g., Millen et al., 2001; Jansen et al., 2007) have revealed that several proteincoding genes (i.e., accD, ycf1, ycf2, rpl22, rps16, rpl23, infA, and ndhF) have been independently lost over the course of angiosperm evolution, these genes were all identified in the eight Panax plastomes. Additionally, the sequence length of the whole plastome, the LSC, SSC, IRs, coding sequences, and non-coding regions of the eight Panax species were quite similar (Table 1). These results suggest that plastome structure and gene content in the genus Panax are highly conserved.

IR expansions often lead to size variations in the plastomes of angiosperms (e.g., Cosner et al., 1997; Plunkett and Downie, 2000; Chumley et al., 2006; Liu et al., 2017; Zhou et al., 2018). However, the IRs/LSC junctions of the eight Panax plastomes were highly conserved: the IRb/LSC boundaries were located between the rpl2 and trnH-GUG genes and the IR regions expanded into rps19 at the IRa/LSC junction, which is similar to other genera within the Araliaceae family (Li et al., 2013). This type of IRs/LSC boundary has been detected in Cornales (Yang and Ji, 2017), but it has not yet been observed in any other orders within the asterids clade (Kim and Lee, 2004; Huang et al., 2014; Downie and Jansen, 2015; Stull et al., 2015; Yao et al., 2016). In contrast to the IRs/LSC junctions, the IRs/SSC boundaries among the eight Panax plastomes were slightly variable, which may have contributed to the overall size variations among Panax plastomes.

During the evolutionary history of a certain lineage, environmental change may impose selective pressures that result in adaptive evolution (Yang and Nielsen, 2000). However, when we examined the eight Panax plastomes, we did not observe any signs or evolutionary fingerprints of positive selection in the proteincoding genes at the generic level (Fig. 5, Supplementary Table S5). In addition, the average Ka values for these genomes were relatively low, and their Ka/Ks values were all less than one (Table 4). These results imply that the plastid protein-coding genes of Panax species may have undergone strong purifying selection during their evolution (Yang and Nielsen, 2000).

4.2. Phylogenetic inferences

We used six datasets (complete plastome, protein-coding exons, LSC, SSC, IRs, and non-coding regions) to reconstruct the phylogeny of Panax. Although earlier studies revealed that both the phylogenetic resolution and the support values of nodes may be considerably improved by more and longer DNA sequences (Rokas and Carroll, 2005; Philippe et al., 2011), our results indicate that the phylogeny based on protein-coding exons generated the highest branch support (BS = 100%, Fig. 6B), even though its sequence length is shorter than that of whole plastomes and LSC regions. Furthermore, our analysis found that non-coding regions across Panax plastomes possessed the highest sequence divergence among the six datasets tested (Table 3), while the phylogeny based on non-coding regions did not have the highest branch support (Fig. 6F). The relatively lower node supports observed in the trees using whole plastomes and LSC regions may be attributed to the faster mutation rates of non-coding regions in the plastome, which, as Wiens (2003) has suggested, produces more evolutionary homoplasy. Accordingly, the protein-coding genes of plastomes can provide accurate information with which to trace the relationships within Panax. Compared to previous single- or multi-locus DNA sequence analyses, our plastid phylogenomic analysis produced higher-resolution nodes, with much higher support values within Panax. Similar findings have been reported in analyses of the plastome-wide protein-coding genes of major flowering plant lineages (Jansen et al., 2007), basal angiosperms (Moore et al., 2007), early diverging eudicots (Moore et al., 2010), commelinid monocots (Barrett et al., 2013), and basal Lamiid orders (Stull et al., 2015).

Phylogenetic analysis using protein-coding exons across the plastome recovered four, well-supported, monophyletic clades within Panax, and the relationships within each clade were also robustly supported (BS = 100%, Fig. 6B), which enabled the phylogenetic backbone of this genus to be recovered. The tree shows that P. ginseng is sister to P. quiquefolius with a high support value (BS = 100%), which is consistent with analysis based on expressed sequence tags (Choi et al., 2013), but differs from analyses of nuclear ITS sequences (Wen and Zimmer, 1996; Choi and Wen, 2000) and plastid intergenic spacers (Shi et al., 2015; Zuo et al., 2017). The sister relationship between P. ginseng and P. quiquefolius was also supported by similar morphological traits (carrot-like main roots), habitats (temperate forests at high latitudes in East Asia and North America), and chromosome number (tetraploid, 2n = 48) (Yi et al., 2004; Choi et al., 2014; Shi et al., 2015).

P. notoginseng is a commonly cultivated medicinal herb in southwest China (Yang et al., 1988). However, its relationship to other Panax species has been disputed (e.g., Wen and Zimmer, 1996; Choi and Wen, 2000; Shi et al., 2015; Zuo et al., 2017). Our results indicate that P. notoginseng, P. japonicus, and P. vietnamensis form a highly supported clade (BS = 100%), which is sister to the clade comprising of P. ginseng and P. quiquefolius; P. trifolius was found to be the earliest diverged clade of Panax, which is identical with previous analyses (Wen and Zimmer, 1996; Lee and Wen, 2004; Shi et al., 2015; Zuo et al., 2017). In previous studies, P. stipuleanatus has repeatedly shown a sister relationship with P. pseudoginseng (see Wen and Zimmer, 1996; Wen, 1999; Lee and Wen, 2004; Zuo et al., 2011, 2015, 2017); however, our phylogenomic analysis resolved P. stipuleanatus and the P. bipinnatifidus species complex as a clade (BS = 100%). Because we did not obtain a sample from P. pseudoginseng in the present study, future studies are required to resolve the position of this species.

The impact of polyploidization on plant speciation during the evolution of Panax is a particularly interesting issue. This genus contains three tetraploid species (2n = 46), P. ginseng, P. japonicus, and P. quiquefolius; furthermore, the P. bipinnatifidus species complex possesses both tetraploid and diploid species/populations (Yi et al., 2004). Our tree topologies clearly indicate that these tetraploid species were scattered in three well-supported clades, suggesting that whole genome duplication events may have occurred independently during the evolutionary history of Panax. This interpretation is supported by the study of Shi et al. (2015).

4.3. Utility of the Panax plastomes

Two plastid protein coding genes, rbcL and matK, and the psbA-trnH intergenic spacer, have been recommended as universal plastid DNA barcodes for land plants (Kress et al., 2005; Hollingsworth et al., 2011). Zuo et al. (2011) then proposed that psbA-trnH and ITS were sufficient for identifying species of Panax. However, we found that variation in rbcL and psbA-trnH was rela-tively low among the eight Panax species (at less than 1%) (Supplementary Tables S3-4). Hence, the universal plastid DNA barcodes rbcL and psbA-trnH may have limited power to identify Panax species. Thus, novel DNA barcodes for this genus are urgently needed.

Based on the sequence variations, we found nine protein-coding regions (ccsA, matK, ndhF, petL, psaI, rpl22, rpoA, rps3, and ycf1) and 10 non-coding regions (atpA-atpF, ccsA-ndhD, infA-rps8, ndhI-ndhA, psbK-psbI, rpl14-rpl16, rpl2-trnH-GUG, rpl22-rps19, rps19-rpl2, and trnY-GUA-trnE-UUC) harboring a high proportion of SNPs. We propose that these plastid DNA regions are potentially useful for identifying Panax species. Among them, matK and ycf1 have been proposed elsewhere as promising DNA barcodes (Hollingsworth et al., 2011; Dong et al., 2015), and rpl14-rpl16 have been widely used for phylogenetic studies (Shaw et al., 2014). In future studies, we will investigate whether or not these plastid DNA sequences can serve as reliable and effective DNA barcodes for rapid species identification of plants in the genus Panax.

Notably, in this study species with more than one plastome available were recovered as well-supported monophyletic groups, and different individuals of the same species generated notable genetic differentiation (Fig. 6). These results suggest that the plastome would be a reliable and accurate barcode for improving the resolution of species identification in the genus Panax. Further studies based on sampling at the population scale are needed to evaluate the efficiency of the plastome as an organelle-scale barcode.

5. Conclusion

The plastome of P. stipuleanatus, an endangered and medicinally important plant, was sequenced and assembled. The genome is 156, 069 bp in length and has a typical quadripartite structure. Comparative analysis showed that the plastomes of Panax are relatively conserved. We investigated the substitution rate of protein-coding regions, which suggest that the plastomes of Panax may have undergone strong purifying selection. We generated a well-supported phylogeny for Panax, in which P. stipuleanatus is sister to the P. bipinnatifidus species complex. Moreover, molecular markers with high sequence divergence were identified, which may be useful for phylogenetic analysis and species identification. Overall, the plastome of P. stipuleanatus will provide valuable genetic information for identifying species, phylogenetic research, as well as resource conservation.

Acknowledgements

This research was financially supported by the Major Program of National Natural Science Foundation of China (No. 31590823) and the National Natural Science Foundation of China (No. 31070297). We are grateful to Jiahui Chen at the Kunming Institute of Botany, Chinese Academy of Sciences, for his help with the data analyses.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2018.11.001.