Introduction

DNA barcoding is an effective tool to identify many plant species rapidly and accurately1,2,3,4. However, there is no single universal barcode that can be successfully used to identify all plants to the species level5. Consequently, two alternative strategies have been proposed to distinguish among plant species: the first one is the use of complete chloroplast genomes6, named ‘super-barcoding’, and the second one is an approach that involves searching for mutational hotspots7, or using comparative plastid analyses to find loci with suitable species-level divergence8,9. Analyses of entire chloroplast genome sequences provide an effective way to develop both of these strategies.

In most angiosperms, the chloroplast genomes are inherited maternally and have a consistent structure, including two inverted repeats (IR), one large (LSC) and one small (SSC) single copy region. The chloroplast genome always contains 110–130 genes that exhibit a range of levels of polymorphism4,5. Thus, chloroplast genome sequence data are extremely valuable for studies of plant population genetics, phylogeny reconstruction, species identification, and genome evolution2,4,5,6.

The Ranunculaceae is a large family, which includes approximately 59 genera and 2500 species. Many plants of Ranunculaceae are pharmaceutically important10. The genus Pulsatilla Adans. consists of about 40 species which are distributed in temperate subarctic and mountainous areas of the Northern Hemisphere10. There are always long, soft hairs covering plants of Pulsatilla species. Most of the flowers of Pulsatilla are large and showy, and therefore the genus has horticultural importance11,12. The flowers are solitary and bisexual. In one flower, there are always six tepals, numerous stamens and carpels, with the outermost stamens resembling degenerated petals, excluding P. kostyczewii11,13,14,15.

In China, there are eleven species of Pulsatilla. Some species of Pulsatilla have been used in traditional Chinese medicine for many years, such as for “detoxification” or “blood-cooling”, because Pulsatilla species contain numerous secondary metabolites, including phytosterols, triterpenoid saponins and anthocyanins16. At the same time, all members of Pulsatilla produced the lactone protoanemonin17,18,19.

In Europe, some species of Pulsatilla are rare, endangered and endemic. Those taxa are protected due to their small populations and disappearing localities, and those species have been placed on the Red Lists of Endangered Species20.

Taxonomically, Pulsatilla is an especially complex and challenging group. In all treatments published before, three subgenera have been recognized: subgenus Kostyczewianae (only one species), subgenus Preonanthus, and the largest subgenus Pulsatilla. However, the intragenic morphological variability of Pulsatilla was especially complicated12. The recognition and identification of wild Pulsatilla species is particularly difficult based on traditional approaches.

Molecular markers are significant to explore the phylogenetic relationships of the genus Pulsatilla. Phylogenetic relationships between Pulsatilla and closely related genera have been dedicated during the past years21,22,23,24,25. Previous studies have attempted to identify these species among Pulsatilla with universal molecular markers, but the species resolution was relatively low15.

In this study, we present seven complete cp genomes from two subgenera of Pulsatilla obtained through next-generation sequencing (NGS) and genomic comparative analyses with four previously published cp genome sequences of Pulsatilla from NCBI, with Anemoclema glaucifolium as the outgroup. We identify microsatellites (SSRs), larger repeat sequences, and highly variable regions, with the aim of developing DNA barcodes and testing the feasibility of phylogenetic analyses of Pulsatilla using the chloroplast genome.

Results and discussion

Chloroplast genome features

We have obtained 1.95 Gb of average NGS clean data for each species, with minimum and maximum values of 1.14 Gb (P. dahurica), and 3.56 Gb (Pulsatilla alpina), respectively. The read number for each species ranged from 6,468,944 (P. dahurica) to 15,816,765 (P. alpina). The average length of the reads was 150 bp on the Illumina Sequencing System. The seven new Pulsatilla cp genomes ranged from 161,501 bp (P. grandis) to 162,669 bp (P. alpina) in length and 151.5 × to 503.4 × coverage. These seven novel Pulsatilla cp genome sequences were submitted to GenBank (Tables 1, 2). Their quadripartite structure is similar to the majority of cp genomes of land plants, which are composed of a pair of IRs (31,184–31,416 bp), separated by the LSC (81,615–82,149 bp) and SSC (of 17,431–17,908 bp) regions26,27 (Fig. 1; Table 2). Previous studies of other angiosperm groups have found that chloroplast genomes are conserved28 or highly polymorphic29,30. These genomes which we reported are highly conserved in gene order, gene content and intron number, which is in accordance with the results from many other taxa26,27. However, in some taxa, e.g. Amorphophallus of Araceae, some genes were lost31, and in others, e.g. Pelargonium32, the structure and gene order diverges from what is reported here and in most other angiosperms. The cp genomes of P. alpina, P. grandis, P. hirsutissima, P. ludoviciana, P. multifida and P. occidentalis had the same GC content of 37.6%, while P. dahurica had a subtle difference (37.5%) compared with the others.

Table 1 Voucher information and GenBank accession numbers for Pulsatilla and outgroups.
Table 2 Summary of complete chloroplast genomes of Pulsatilla.
Figure 1
figure 1

Gene map of the Pulsatilla chloroplast genome. Dashed area in the inner circle indicates the GC content of the chloroplast genome. LSC, SSC and IR mean large single copy, small single copy and inverted repeat, respectively. Genes belonging to different functional groups are color-coded as indicated by icons on the lower left corner. The red line on the outside of the gene is three inversions.

Chloroplast genome comparison

In most angiosperms, the IR regions of cp genomes of angiosperms are highly conserved, but the expansion and contraction of IR region boundaries are ever present33,34. At the same time, several lineages of land plant chloroplast genomes show great structural rearrangement, even loss of IR regions or some gene families35. The expansion and contraction in IRs are significant evolutionary events, because they can change gene content and chloroplast genome size30,36. Expansion of the IRs has been reported in Araceae36,37. Sometimes, the size of LSC increases and that of SSC decreases, becoming only 7000 bp in Pothos38. At the same time, a linear chloroplast genome was also reported in some groups, e.g. maize35,39. Expansion and contraction of the IR regions can also lead to duplication of certain genes or conversion of duplicate genes to single copy, respectively30,36. Changes in the size of the IRs can also cause rearrangement of the genes in the SSC as recently observed in Zantedeschia36.

The Pulsatilla chloroplast genomes were compared to previously published data and showed typical Anemoneae (Ranunculaceae) genome structure37,39. As reported for Anemoclema, Anemone, Clematis and Hepatica, the IR regions of genus Pulsatilla are roughly 4.4 kb longer than those of other genera of the family Ranunculaceae, such as Aconitum, Coptis, Thalictrum, Megaleranthis, Ranunculus, and Trollius37,39. The gene orders located within the IR-SSC and IR-LSC boundaries are similar among tribe Anemoneae but different from those of other genera of Ranunculaceae (Fig. 2, Fig. S1).

Figure 2
figure 2

Comparisons of LSC, SSC, and IR region borders among the seven Pulsatilla chloroplast genomes.

We compared the IR/SC boundary regions of Pulsatilla, and the junction positions are very similar and conserved within genus Pulsatilla. In the four boundary regions (LSC, IRa, SSC, IRb) of seven Pulsatilla cp genomes, the LSC/IRa and IRb/LSC border was in the intergenic region, and the adjacent genes is rps36, rps8 and rps4, respectively. The genes ycf1 andψycf1 have crossed the SSC/IRb and IRa/SSC boundary, respectively, which was also found in Monsteroideae (Araceae)28. The pseudogene ycf1 has been found in other groups30,36. The IR regions were highly conserved, with nucleotide diversity values in those regions less than 2%.

Sequence divergence

Multiple alignments of plastid genomes were performed to investigate levels of genome divergence. Based on MAFFT analysis, there are three inversions in LSC of Pulsatilla, same as the tribe Anemoneae (Fig. 1, Fig. S1)37. The mVISTA analysis has revealed high sequence similarity across the coding region and there exists more variability in non-coding regions. Sequence identity among the seven species was 96.68–98.66%. The number of nucleotide substitutions and sequence distance (Pi) were the highest (1214, 0.0063) between P. alpine and P. dahurica, with the lowest (112, 0.0005) between P. vernalis and P. patens (Fig. 3; Table 3).

Figure 3
figure 3

Visualized alignment of the eight Pulsatilla chloroplast genomes. The mVISTA-based identity plots show the sequence identity among the seven chloroplast genomes, with P. chinensis serving as a reference. Blue represents coding regions, and pink represents non-coding regions.

Table 3 Numbers of nucleotide substitutions and sequence distance (Pi) in eleven complete cp genomes.

Identification of highly variable regions

Chloroplast genome markers, especially several universal chloroplast regions, have been widely used in plant systematics and identification at multiple taxonomic levels. Highly suitable polymorphic chloroplast loci have been identified and designed as unique markers in different groups28,40. However, relationships within the genus Pulsatilla have not been well resolved because of the low polymorphism of these universal markers15. In order to facilitate identification of closely related species of Pulsatilla, we sought to identify highly variable regions of the chloroplast genome, as previously described9,27,41,42,43,44. As a result, we identified nine divergent hotspot regions, including six intergenic spacer regions (rps4-rps16, rps16-matK, ndhC-trnV, psbE-petL, ndhD-ccsA, ccsA-ndhF) and four protein-coding regions (ycf1, ndhF, ndhI) (Fig. 4; Table 4). Most commonly employed loci, e.g. trnL-trnF, trnH-psbA were not selected in our finding. The nine highly variable regions included 684 variable sites, including 181 indels. However, these indels are not suitable for the phylogenetic inference because Maximum likelihood model used only substitutions not indels28. Their nucleotide diversity values ranged from 0.00802 to 0.02212. The region of ccsA-ndhF showed the highest variability, the next most variable regions were rps4-rps16, ndhC-trnV, and psbE-petL. The diversity level of two protein-coding regions (ycf1, ndhF) was the lowest.

Figure 4
figure 4

Sliding window analysis of the whole chloroplast genomes of Pulsatilla taxa.

Table 4 Sequence characteristics of eight high variable regionsamong eleven complete cp genomes of Pulsatilla.

Among the nine divergent hotspot regions, the ndhI is difficult to align. There are large numbers of indels in ndhI and the intergenic spacer between ndhI and ndhG, these regions were not considered suitable for the phylogenetic inference of the Pulsatilla. Thus, we selected eight regions, four (rps4-rps16, rps16-matK, ndhC-trnV, psbE-petL) in the LSC and four (ndhD-ccsA, ccsA-ndhF, ycf1, ndhF) in the SSC, with relatively high variability as potential molecular markers for the study of species identification and phylogeny in Pulsatilla. Five hotpots were found in chloroplast genome of Veroniceae (Plantaginaceae), and two universal marker, trnH-psbA and matK were identified, respectively45,46. Ten highly variable regions were selected as potential molecular markers for Fritillaria, including ycf144,47, which was also selected in this study. Sequences of these variable regions founded in this study could be regarded as potential molecular markers for species identification and evolutionary studies and have been shown to be valuable for studies in other groups (e.g., Fritillaria)44.

SSRs and large repeat sequences

Oligonucleotide repeats play an important role for generating indels, inversion and substitutions29. Repeat sequences in the chloroplast genome could provide valuable information for understanding not only the sequence divergence but the evolutionary history of the plants48,49,50. We have detected five types of large repeats (forward, reverse, palindromic, complement and tandem repeats) in the seven Pulsatilla cp genomes. Among them, the most common repeat types are forward and palindromic repeats, followed by reverse repeats, and only little complement repeats were found in Pulsatilla cp genomes (Fig. 5A). Most of the repeats were short, ranging from 30–49 bp (Fig. 5B).

Figure 5
figure 5

Analyses of repeated sequences in seven newly sequenced chloroplast genomes. (A) Number of five repeat types; (B) frequency of four repeats by length; (C) frequency of microsatellites by base composition; (D) frequency of microsatellites by types; (E) frequency of microsatellites by length; (F) number of all repeats by location.

We also identified multiple microsatellite repeats, also known as simple sequence repeats (SSR) or short tandem repeats (STR)49. Due to their codominant inheritance and high variability, SSRs are robust and effective markers for species identification and population genetic analyses49,50,51,52,53. Most of the mononucleotide repeats were composed of A/T. The other microsatellites types were also dominated by AT/TA, with very little G/C (Fig. 5C). In this study, plentiful microsatellite loci were found through the comparative analysis of Pulsatilla cp genome sequences. In total, we detected six types of microsatellite (mononucleotide, dinucleotide, trinucleotide, tetranucleotide, pentanucleotide and hexanucleotide repeats) based on the comparison of seven Pulsatilla cp genomes (Fig. 5D). Each Pulsatilla cp genome had 69–87 microsatellites. The lengths of repeat motifs of these microsatellites ranged from 10 to 21 bp (Fig. 5E). Among the four structural regions in the cp genomes, most of the repeats and microsatellites were distributed in LSC, followed by SSC, and fewest in IRa/IRb (Fig. 5F), which were also reported in other studies in angiosperms29,30. These SSRs and repeat sequences are uncorrelated with genome size and phylogenetic position of the species36, but will provide important information for further studies of phylogenetic reconstruction and infra- and inter-specific genetic diversity54,55.

Phylogenetic analyses

Chloroplast genomes have been widely used and have made significant contributions to phylogeny reconstruction at different taxonomic levels in plants7,9,24,35. To better clarify the evolutionary relationships within Pulsatilla, we used each data set to construct phylogenetic trees using the ML analytical methods. We also construct phylogenetic trees with those eight highly variable regions using the ML, MP analytical methods. All tree topology structures were identical. Therefore, here we presented the phylogenetic studies using the ML tree with the support values from the MP analyses recorded at the corresponding nodes (Fig. 6). The phylogenetic tree based on all data sets (except the IR and SSC regions) from the complete plastid genome sequences yielded the same topology. The phylogenetic tree based on chloroplast genome differed from that of the DNA barcode combination rbcL + matK + trnH-psbA, but with higher support values. The phylogenetic trees based on data from complete plastid genome sequences showed that the species of Pulsatilla formed a monophyletic group which in turn includes two strongly supported (bootstrap = 100) clades. One clade comprised P. alpina and P. occidentalis, members of subg. Preonanthus. The other comprised two subclades: (1) members of P. hirsutissima, P. ludoviciana, P. multifidi, P. patens and P. vernalis, and (2) species of P. chinensis, P. dahurica, P. grandis and P. pratensis. All the species of the two subclades are members of the subg. Pulsatilla. These results were congruent with our former results based on universal markers15.

Figure 6
figure 6

Phylogenetic relationships of the eleven Pulsatilla species inferred from maximum likelihood (ML). Including whole chloroplast genome, rbcL + matK + trnH-psbA, LSC region, coding region, SSC region, IR region, and the concatenation of the eight highly variable regions mentioned in Table 4 (Numbers above nodes are support values with ML bootstrap values on the left, and MP bootstrap values on the right).

In phylogenetic analyses, compared to the combination of barcodes, the full chloroplast genome sequence data formed distinct clades with high bootstrap support, improving the inadequate resolution of barcodes combination. The LSC regions and coding regions have the same topology structures with robust support. However, sequencing of genomic DNA is still expensive. It is necessary to utilize variation within chloroplast regions for rapid species-specific assay5,9,31,33,42. Here we found that phylogenetic inference based on highly variable regions yielded a tree with the same topology as the one recovered based on complete chloroplast genome sequences, demonstrating the high utility of hotspots of variability for species identification and phylogenetic analysis. More samples and laboratory works are needed in the future to increase the number of these variable regions available for study.

Conclusions

In this study, we generated complete chloroplast genomes of seven species of Pulsatilla and compared them to four previously published cp genome sequences of Pulsatilla. Chloroplast genomes of Pulsatilla share many features with those of other angiosperms. Informative differences between cp genomes of Pulsatilla were most evident in inversions in the large single copy region and expansion of the inverted repeat region. We identified multiple potentially valuable genetic markers, including large repeat sequences, numerous SSRs, and eight highly variable regions. Genetic markers have provided a reference for the improvement of plants fingerprints and the identification of similar Pulsatilla. In addition, it will be better to construct the phylogenetic relationships of Pulsatilla species using these highly variable regions and genetic markers in the future. Totally, this study provides a basis for future studies of horticultural cultivation, conservation, population genetics, phyletic evolution, development of DNA barcodes, and diverse research in Pulsatilla.

Materials and methods

Plant material, DNA extraction and sequencing

Fresh leaves of P. dahurica were collected from Jilin province of China and dried with silica gel. Dry leaves of other six Pulsatilla species were taken from herbarium specimens. We extracted total genomic DNA with the DNeasy Plant mini kits (QIAGEN, Guangzhou, China). The genomic DNA was sequenced using the Illumina Miseq platform (Illumina, San Diego, CA, USA).

Chloroplast genome assembly and annotation

Whole chloroplast genome sequencing was done for the seven species of Pulsatilla. For each species, high-quality Illumina sequencing reads were assembled into scaffolds with de novo sequence assembly software Spades, SOAPdenovo and CLC Genomics Workbench v.6.5 (CLC Bio), respectively56. We checked the contigs against the reference genome of P. chinensis (MG001341), using BLAST (https://blast.ncbi.nlm.nih.gov/) and oriented aligned contigs according to the reference genome. We mapped all the raw reads back to assembled sequences to check the assembly and then constructed the complete cp genomes using Geneious v.9.057. We submitted all the newly sequence data in raw format (fastq) and obtained SRA accesions (Table 1).

Annotations of cp genome sequences were performed using Plastid Genome Annotator (https://github.com/quxiaojian/PGA)58 and adjusted in Geneious v.9.0. To verify the accuracy of the annotations, we also used GeSeq59 to annotate each chloroplast genome in this study. We checked every boundary of tRNAs using tRNAscan-SE v.2.060. The circular genome maps were generated in OGDRAW (https://ogdraw.mpimp-golm.mpg.de/)61.

Genome comparisons

We aligned the cp genomes of Pulsatilla using multiple alignment of MAFFT v762 and manually edited in Geneious v.9.0. The contraction and expansion of inverted repeat regions were also examined among the seven species (excluded P. chinensis, MG001341) of the genus Pulsatilla using Irscope63. Then, we performed multiple alignments of the eight genomes of Pulsatilla in the mVISTA program64 under Shuffle-LAGAN mode, with default parameters for other options, using the annotation genome of P. chinensis as a reference, with the aim of comparing and visualized the similarities and differences among different plastid genomes.

To analyse chloroplast genome organisation and gene arrangement, we perform the analyses of collinear blocks with Mauve v 2.3.165 plugin in Geneious v.9.0, including only one copy of the IR, assuming collinear genomes for the full alignment. Detailed gene inversions were identified by comparing the gene order of Pulsatilla samples and Anemoclema to Berberis.

To observe the plastid genome divergence and determine parsimony informative sites, we conducted sliding window analysis after alignment to determine the nucleotide diversity (Pi) of the cp genome using DnaSP v5, with 200 bp of step size and 600 bp window length66. We defined hotspots as those regions with a higher value of Pi. We computed the variable sites across the complete cp genomes and the sequence characteristics of hotspots by DnaSP v5.0.

Repeated sequences identification

We identified repeat sequences, including palindromic, reverse and forward repeats, using the online software REPuter (https://bibiserv.cebitec.uni-bielefeld.de/reputer), with the following settings: Hamming distance of 3 and minimum repeat size of 30 bp67. We used the online program Tandem Repeats Finder (https://tandem.bu.edu/trf/trf.html) to find the tandem repeat sequences, in which the similarity percentage of two repeat copies was at least 90% and the minimal repeat size was 10 bp. The alignment parameters for match, mismatch, and indels were set at 2, 7, and 7, respectively. We identified microsatellites (SSRs) by MISA68 with thresholds of 10, 5, 4, 3, 3, and 3 for mono-, di-, tri-, tetra-, penta-, and hexa-nucleotide, respectively.

Phylogenetic analyses

For the purpose of reconstructing the phylogenetic relationships, four published complete cp genome sequences from the genus Pulsatilla and Anemoclema were also included in our analyses. The monotypic genus Anemoclema (MH205609) was selected as the outgroup. Because molecular evolutionary rates among the different cp genome regions are diverse, analyses of phylogenetic relationships were performed based on the following seven datasets: (a) the complete cp genome sequences; (b) coding genes (CDS); (c) one inverted repeat (IR) region (IRb); (d) the large single copy region (LSC); (e) the small single copy region (SSC), (f) the consensus sequences of eight highly variable regions; and (g) the DNA barcodes combination (rbcL + matK + trnH-psbA). We applied Maximum Likelihood (ML) analysis for each of the seven datasets to construct tree-sets. Maximum Parsimony (MP) analyses were also applied for the consensus sequences of eight highly variable regions and the DNA barcodes combination.

We conducted ML analyses with RAxMLHPC2 v.8.0.969 on the Cyberinfrastructure for Phylogenetic Research (CIPRES) Science Gateway v.3.370. Then the analysis of 1000 rapid bootstrap replicates (-x) was followed by a search for the best-scoring ML tree in one program (-f a). The best-fit model for nucleotide and amino acid sequences were evaluated using jModelTest 271. We applied the GTR + G model to nucleotide data for both bootstrapping and best-tree searching phases, with other parameters as the default settings. We performed the maximum Parsimony (MP) analysis on PAUP* v.4.0b1072. All the characters were treated as unordered and equally weighted. The heuristic search specified 1000 random sequence addition replicates with TBR branch swapping, saving only 10 trees per replicate. We obtained the strict consensus tree from all the most-parsimonious trees (MPTs) detected during the search. We calculated bootstrap percentages (BP) from 10,000 rapid bootstrap replicates, each comprising 10 random sequence addition replicates and saving only one tree per replicate.