Plant materials from ten Aconitum spp. were collected from the Weixi County, Lijiang County, and Zhongdian County, Yunnan Province, in September 2020 (Table 1; Figure 1). Detailed information about the materials used and sequences obtained in the study is provided in Table 1. Fresh leaf material without lesions was collected and stored dry with discoloured silica gel. Professor Conglong Xia of Dali University identified the species according to the morphological characteristics recorded in the Flora of China. The voucher specimens were deposited in the Plant and Medicinal Herbology, College of Pharmacy, Dali University. The voucher specimens number was recorded as WT001∼WT0010.

Total genomic DNA was extracted from clean leaves from samples frozen at −80°C using the E.Z.N.A® Plant DNA kit (OMEGA, Beijing). Shanghai Origingene Co., Ltd. did the work. The DNA quality was assessed by electrophoresis in a 1% (w/v) agarose gel. The OD₂₆₀/₂₈₀ values ranged from 1.8 to 2.0, and >2 μg of DNA was equally pooled from individuals of the ten species to generate shotgun libraries. DNA samples were randomly sheared, incubated with fragmentation buffer, and broken into 300–500 bp fragments in an M220 focused ultrasonicator (Covaris, Woburn, MA, United States). The A & B connectors at both ends of the DNA fragment were connected, the segments were screened, and the self-connecting segments of the connector were removed. Fragment screening by electrophoresis in agarose gel keeps the fragment with A connector at one end, and B connector at the other end; the subsequent additions of NaOH to denatures produces single-stranded DNA fragments. This DNA was further paired-end sequenced using the Illumina Novaseq 2500 (Shanghai Origingene Co., Ltd., China) platform.

Raw image data obtained by Illumina sequencing was converted into FASTQ format sequencing data files by Base Calling (Paul, 2015). The FASTQ data file used FastQC v0.11.4 (www.bioinformatics.babraham.ac.uk/bugzilla/) for quality control. The sequenced reads were spliced with several iterations using NOVOPlasty (https://github.com/ndierckx/NOVOPlasty) stitching software (Dierckxsens et al., 2017). The seed sequence required for the assembly was MN556604.1 (complete chloroplast genome of Delphinium grandiflorum), and the optimal assembly result was obtained.

PGA software (https://github.com/quxiajian/PGA) was used for gene annotation (Qu et al., 2019), and manual corrections were made for all sample annotation results. The annotated chloroplast genome GB file was submitted to the online tool (https://chlorobox.mpimp-golm.mpg.de/OGDraw.html) to generate a physical map (Lohse et al., 2007). Finally, these TBL files were generated to submit the sequences to NCBI. The complete and correct CP genome sequences of the ten Aconitum species were deposited in the GenBank database (https://www.ncbi.nlm.nih.gov/genbank/) of NCBI, and accession numbers are OM289057.1 (A. ouvrardianum), OM289058.1 (A. delavayi), OM289059.1 (A. ramulosum), OM328065.1 (A. vilmoinianum), OM328066.1 (A. episcopale), OM328067.1 (A. stapfianum), OM328068.1 (A. nagarum), OM328069.1 (A. weixiense), OM328070.1 (A. duclouxii), and OM328071.1 (A. stylosum), respectively.

The CodonW software (University of Texas, Houston, TX, United States) was used to investigate the distribution of codons based on the relative synonymous codon usage (RSCU) ratio and the effective number of codon usage (ENC) to analyze the codon usage preference (Sharp et al., 1986).

The FASTA files of whole chloroplast genome sequences of 42 Aconitum species were sorted out. Simple sequence repeats (SSRs) were detected using the MISA software (http://pgrc.ipk-gatersleben.de/misa/) (Beier et al., 2017). The minimum values for the number of repeats of mono -, di -, tri -, tetra -, penta -, and hexanucleotide repeats were set to 10, 5, 4, 3, 3, and 3 respectively. The sizes and locations of repeat sequences in the CP genomes of the 42 Aconitum spp. were identified using REPuter (http://bibiserv.techfak.uni-bielefeld.de/reputer/) (Kurtz et al., 2001) with the parameters set to a similarity percentage of scattered repeat copies ≥90%, a minimal repeat size of 30 bp, the hamming distance of 3 and maximum computed repeats of 5,000.

The mVISTA online program (https://genome.lbl.gov/vista/mvista/submit.shtml) (Frazer et al., 2004) was used in Shuffle-LAGAN model to compare 42 Aconitum species CP genome using A. vilmorinianum CP genome as a reference. The IRscope (https://irscope.shinyapps.io/irapp/) online program (Amiryousefi et al., 2018) was used to obtain IR regions comparative analysis of the chloroplast genome.

The DnaSP software was used for sliding window analysis (Rozas et al., 2017), the nucleotide diversity values (Pi) were calculated, and interspecific high variation sequences (hotspots) were screened according to the analysis results. The size of the windows length was set to 600 bp, and the step size was set to 200 bp.

To determine the phylogenetic positions of the ten Aconitum species within Ranunculaceae, we analyzed the CP genomes of 44 species, encompassing 34 additional taxa within this lineage. The CP genome sequences of 34 species were downloaded from Genebank database. The complete CP genome sequences and PCGs were used to reconstruct the Aconitum species phylogenetic tree. The chloroplast genome sequences were aligned using MAFFT online tool (Katoh and Standley, 2013). We used the CP genomes of Delphinium anthriscifolium Hance (MK253461.1) and Delphinium grandiflorum L. (NC_049872.1) as outgroups. The MEGA X software (Kumar et al., 2018) was used to construct phylogenetic trees employing 44 CP genomes sequences based on neighbour-joining (NJ) with 1,000 bootstrap replicates, and the model was Kimura 2-parameter. The IQ-tree software (Nguyen et al., 2015) was used to construct phylogenetic trees employing 44 CP genomes sequences based on maximum likelihood (ML) with 1,000 bootstrap replicates.

The genome sequences assembled using the reads obtained from the Illumina sequencing platform ranged from 155,475 bp for A. stylosum to 155,921 bp for A. vilmorinianum (Table 2). The genome exhibited a typical cyclic tetramer structure, including the SSC and LSC regions separated by two IR regions (Figure 2). The sequence length of LSC regions ranged from 86,098 bp for A. stylosum to 86,524 bp for A. vilmorinianum, the sequence length of SSC regions ranged from 16,903 bp for A. stylosum to 16,983 bp for A. nagarum, and the sequence length of IRs regions ranged from 52,218 bp for A. duclouxii to 52,476 bp for A. ouvrardianum (Table 2).

The chloroplast genomes of ten Aconitum species can be divided into the coding and non-coding regions; the coding region sequence lengths ranged from 79,014 to 78,996 bp, and the non-coding region sequence lengths ranged from 76,569 to 77,475 bp. The ten Aconitum chloroplast genome sequences all encoded 129–132 genes (including IR region duplicates), and if 18 duplicates in the IR regions were excluded, a total of 111–114 genes were encoded, including 83 to 85 PCGs (seven duplicated genes), 37 tRNAs (seven duplicated genes), eight rRNAs (four duplicated genes). It is noteworthy that two pseudogenes (ycf1 and rps19) were found in A. ouvrardianum, A. delavayi, and A. ramulosum, whereas only the pseudogene ycf1 was found in the remaining species. Moreover, the GC contents of the ten Aconitum CP genomes were 38.1%, and the IR regions (43%) had higher GC contents than the single-copy regions (LSC: 36.2% and SSC: 32.6%). Furthermore, we also analyzed the length and frequency of four bases that included adenine (A), thymine (T), guanine (G), and cytosine (C) in the chloroplast genome of ten Aconitum species. The results showed that the lengths of adenine ranged from 47,585 to 47,840 bp, the lengths of cytosine ranged from 29,781 to 30,212 bp, and the lengths of guanine ranged from 29,058 to 29,524 bp, and the lengths of thymine ranged from 48,391 to 48,825 bp, respectively. A and T were used significantly more frequently than G and C.

The base composition of the entire coding region and different position codons in the chloroplast genomes of the ten Aconitum species were analyzed (Table 3). The coding region length ranged from 63,372 bp for A. duclouxii to 79,014 bp for A. delavayi. The total AT content in the coding region was 60.6–61.7%, and the total GC content was 38.3–39.4%. The length of the first, second, and third codons ranged from 21,124 bp to 26,338 bp, and GC contents ranged were 45.9–47.3% (first codon), 38.5–39.4% (second codon), and 30.7–31.5% (third codon), respectively. Moreover, the third codon site has AT preference.

Furthermore, we also analyzed the base composition in different regions of the ten Aconitum species chloroplast genomes (Supplementary Table S1). The AT content of SSC was the highest, followed by the LSC region and IRs region. The AT content of the protein-coding gene region was the highest, followed by the tRNAs region and rRNAs region. The results indicated that the AT content of the chloroplast genome was higher than the GC content and also demonstrated the codon preference for using bases A and T(U).

Chloroplast genomes of ten Aconitum species were annotated, and their genes were functionally classified (Table 4). The number of chloroplast genome coding genes in the ten Aconitum species ranged from 129 to 132, including 83–84 protein-coding genes (PCGs), 37 tRNA genes, and eight rRNA genes. The first group of genes involved in transcription and translation consists of 69–70 genes, including small subunit of ribosome genes (rps2, rps3, rps4, rps7 × 2, rps8, rps11, rps12 × 2, rps14, rps15, rps16, rps18, and rps19), ten large subunit of ribosome genes (rpl2 × 2, rpl14, rpl16, rpl20, rpl22, rpl23 × 2, rpl33, and rpl36), one transcription initiation factor genes (infA), four RNA polymerase genes (rpoA, rpoB, rpoC1, rpoC2), eight rRNA genes (rrn16 × 2, rrn23 × 2, rrn4.5 × 2, and rrn5 × 2), and 37 tRNA genes (trnA-UGC × 2, trnC-GCA, trnD-GUC, trnE-UUC, trnF-GAA, trnfM-CAU, trnG-GCC, trnG-UCC, trnH-GUG, trnI-CAU × 2, trnI-GAU × 2, trnK-UUU, trnL-CAA × 2, trnL-UAA, trnL-UAG, trnM-CAU, trnN-GUU × 2, trnP-UGG, trnQ-UUG, trnR-ACG × 2, trnR-UCU, trnS-GCU, trnS-GGA, trnS-UGA, trnT-GGU, trnT-UGU, trnV-GAC × 2, trnV-UAC, trnW-CCA, and trnY-GUA). The second group of 45 genes involved in the photosynthesis, including 20 Photosystem I and Photosystem II genes (psaA, psaB, psaC, psaI, psaJ, psbA, psbB, psbC, psbD, psbE, psbF, psbH, psbI, psbJ, psbK, psbL, psbM, psbN, psbT, and psbZ), six Cytochrome b/f complex genes (petA, petB, petD, petG, petL, petN), six ATP synthase genes (atpA, atpB, atpE, atpF, atpH, atpI), one ATP-dependent protease subunit gene (clpP), one RubiscoCO large subunit gene (rbcL), and 12 NADH oxidoreductase genes (ndhA, ndhB × 2, ndhC, ndhD, ndhE, ndhF, ndhG, ndhH, ndhI, ndhJ, ndhK). The third group of four genes involved in the biosynthesis of amino acids, fatty acids, etc., included one Maturase gene (matK), one Envelop membrane protein gene (cemA), one Subunit Acetyl-CoA-Carboxylate gene (accD), and one c-type cytochrome synthesis gene (ccsA); another group of unknown function genes has 6–8. Among the ten Aconitum species, A. ouvrardianum, A. delavayi, and A. ramulosum contained ycf15 and rps16 genes, and the other seven species did not have ycf15 and rps16 genes.

In eukaryotic and semi-prokaryotic systems, gene expression occurs in the cytoplasm and organelles of the nucleus, respectively. Introns play a vital role in the regulation of gene expression. Previous studies have proved that introns can improve the expression level of extraneous genes in eukaryotic genomes (Cui, 2020). A total of 17 genes contained introns in the ten Aconitum CP genomes, clpP, ycf3, and rps12 genes contained two introns and three exons, whereas the other 14 genes all contained one intron and two exons (Supplementary Table S2). In addition, the introns and exons contained by the genes of ten Aconitum were essentially the same or had minor differences in length, as shown by the results in Supplementary Table S2.

The codon usage and codon recognition patterns of the ten Aconitum CP genomes are shown in Figure 3 and Supplementary Table S3. The CP protein-coding genes of these ten species contained 61 codons encoding 20 amino acids. Leucine (10.4%–10.5%) was the most frequent, and cysteine (1.18%–1.19%) was the least, which was consistent with the results of amino acids coding of chloroplast genomes in most plants (Somaratne et al., 2020; Wang Z. et al., 2021). The statistical results of the frequency of codon usage in the chloroplasts of the ten Aconitum species showed that the coding regions (CDS) comprised 26,086 codons in A. ouvrardianum to 26,338 codons in A. vilmorinianum (Supplementary Table S5). Among amino acids encoded by the codon, Leu was the most frequent amino acid, encoding 2,701 times in A. ouvrardianum to 2,732 times in A. stylosum and A. nagarum; The least frequently coded amino acid was Cys, only encoding 302 times in A. ouvrardianum, A. ramulosum, and A. delavayi to 307 times in A. vilmorinianum. The results are consistent with those in Figure 3.

Furthermore, we analyzed the relative synonymous codon usage (RSCU) and codon base composition of protein-coding gene sequences (CDS) of ten Aconitum species. There were 64 codons in CDS sequences of ten Aconitum species genes, among which 61 codons encoded 20 amino acids and the other three codons were stop codons (Figure 4B). An RSCU value <1.00 indicates that the codon usage frequency is lower than expected, whereas an RSCU value >1.00 indicates that the codon usage frequency is higher than expected (Sharp et al., 1986). In this study, the RSCU values of 31 codons were greater than 1, RSCU values of 31 codons were less than 1, and RSCU values of two codons were 1 (Figure 4B). Among the 64 codons, there were 16 codons, each ending with A, U, G, and C. Meanwhile, among the codons with RSCU value >1, 13 codons ended in A, 16 codons ended in U, one codon ended in G, and one codon ended in C (Figure 4A). These results indicated that the chloroplast genome codon of Aconitum species prefers to end in A/U.

SSRs are widely present in CP genomes, consisting of 1-6-nucleotide repeat units, and are valuable molecular markers of high variation within the same species (Jiao et al., 2012). In order to explore the distribution and differences of SSRs among Aconitum species, we detected SSRs in the chloroplast genomes of 42 Aconitum species using MISA. Ten of them were the species sequenced in this study, and the sequences of the remaining 32 species were downloaded from NCBI. This study analyzed the number, type, and regional distribution of SSRs in 42 Aconitum chloroplast genomes (Figure 5). The total SSR loci in chloroplast genomes of 42 Aconitum species ranged from 48 in A. volubile to 79 in A. coreanum and A. reclinatum.

Furthermore, we analyzed SSRs loci in different regions of the chloroplast genome (LSC/SSC/IRs). The results showed that the distribution of SSRs significantly varied in different regions (Figure 6). The number of SSRs loci in the LSC region ranged from 38 to 65, accounting for 66.67%–85.92% of the total number of SSRs; SSC regions ranged from six to fifteen, accounting for 9.38%–20.00%; Two IR regions ranged from one to ten, accounting for 1.32%–17.54%. Most of the SSR sequences in the chloroplast genomes of 42 Aconitum plants were composed of mononucleotide and dinucleotide repeat units. The number of mononucleotide repeats ranged from 19 in A. volubile to 52 in A. finetianum. These were followed by dinucleotide (10–17), trinucleotide (6–12), tetranucleotide (6–9), pentanucleotide (1–5), and hexanucleotide (1–3).

We analyzed the repeat units of different types of repeat sequences in ten Aconitum species and found that A/T (96.67%–100% of the total SSRs of each species) and AT/TA (100%) dominated the repeat sequence types of mononucleotide and dinucleotide SSRs, respectively. The repeated trinucleotide SSRs were mainly composed of A/T base combinations (AAT/ATT, 71.43%–85.71%) (Supplementary Table S5). These results were consistent with previous reports that CP SSRs usually consist of short poly-A or poly-T repeats and rarely contain tandem G or C repeats in many plants (Kuang et al., 2011).

In this study, the analysis results of long repeats in chloroplast genome sequences of 42 Aconitum species showed 1,467 long repeats in all species, and each species had 17–77 long repeats, respectively (Figure 7). Among them, A. piepunense has the most significant number of long repeats, including 18 forward, 19 palindromic, 35 reverse repeats, and five complement repeats; A. monanthum is the least, including six forward and 11 palindromic repeats. Among the 42 Aconitum species, A. duclouxii, A. barbatum var. hispidum, A. barbatum var. puberulum, A. chiisanense, A. finetianum, A. coreanum, A. piepunense, A. pseudolaeve, A. puchonroenicum, A. quelpaertense, A. reclinatum, A. scaposum var. vaginatum, and A. scaposum which contained 30–39 bp complement repeat sequence.

Of the 42 species, the length range of the repeats was divided into five types (30–39, 40–49, 50–59, 60–69, ≥70 bp), with the highest number of repeats of 30–39 bp in length. All Aconitum species contained 30–39 bp repeat sequences, accounting for 41.18%–85.71% of the total repeat sequences. Among them, 30–39 bp long repeats in A. piepunense, and A. monanthum are the most and least abundant, respectively.

The 60–69 bp repeat sequences were the least. Only in A. nagarum, A. duclouxii, A. ouvrardianum, A. austrokoreense, A. brachypodum, A. chiisanense, A. ciliare, A. contortum, A. flavum, A. hemsleyanum, A. jaluense subsp. jaluense, A. japonicum subsp. napiforme, A. longecassidatum, A. pendulum containing repeats of 60–69 bp in length, accounting for 3.57% , 3.70%, 2.22%, 4.00%, 3.70%, 3.45%, 5.26%, 4.17%, 4.55%, 3.13%, 5.00%, 5.00%, 3.45%, and 3.57% of the total repeat sequence.

During the adaptive evolution of plant species, the expansion and contraction of the IR boundary of the chloroplast genome resulted in different levels of sequence replication at the four boundaries (LSC/SSC/IRa/IRb) (Hu and Zhang, 2021). Analysis of IR boundaries of chloroplast genome sequences of 42 Aconitum species showed that the four boundaries, the structure of the chloroplast genome, and the connection between IR regions were very conserved; LSC/IR and SSC/IR boundary distribution genes included rps19, rpl22, rpl2, ycf1, ndhF, trnH, and psbA (Figure 8). The rps19 gene in the chloroplast genomes of 42 Aconitum species showed various degrees of contraction and expansion at the LSC/IRb boundary. And then, the rps19 gene in chloroplast genomes of most Aconitum species exhibited a 3 bp protrusion in the IRb region. As the ycf1 gene straddles SSC/IRa boundary, a pseudogene ψycf1 was generated in the IRb region. No ycf1 pseudogene was found in A. reclinatum, A. scaposum var. vaginatum, A. sinomontanum, A. finetianum, A. barbatum, and A. angustius at the IRb region.

The length of the ndhF gene in A. duclouxii was 2,237 bp, extending across the SSC/IRa boundary to 30 bp in the IRa region, and the other Aconitum species all had ndhF gene of different lengths in the SSC region. These results indicated that some variation occurred in A. duclouxii. The rpl22, trnH, and psbA genes in the chloroplast genomes of the 42 Aconitum species were all located in the LSC region, the rpl2 gene was all located in the IRa region, and these genes were essentially unchanged in length across species, which can be used as a common feature of Aconitum chloroplast genomes for species identification.

In this study, the chloroplast genome of A. vilmorinianum was used as the reference sequence to compare the chloroplast genome in pairs with that of the other 41 species. The Shuffle—LAGAN model of mVISTA was used to map the chloroplast genomes of 42 Aconitum species (Figure 9) and observe the approximate gene composition and sequence. The results showed that the chloroplast genome sequences of 42 Aconitum plants were similar, revealing that the chloroplast genome of Aconitum was highly conserved. The IR region of the Aconitum chloroplast genome was more conserved than LSC and SSC regions, the coding region was more conserved than the non-coding region, and the variation degree of the intergenic spacer region was greater than that of the gene region. The intergenic spacer region trnK-UUU-trnQ-UUG, trnS-GCU-trnG-UCC, atpH-atpI, petN-psbM, rps18-rpl20, rpl16-rps3, trnL-CAA-trnL-CAA, ndhB-trnL-CAA had the largest variation; The protein-coding gene region was the most conserved, and only ycf1 and ycf2 had the largest variation; Four rRNA genes (rrn4.5, rrn5, rrn16, and rrn23) were the most conserved. These highly variable genes and intergenic spacers can be used as molecular markers for 42 Aconitum species phylogenetic and population genetics studies.

The recognition of highly variable sites in the whole chloroplast genome can be used as molecular markers for species identification and phylogenetic studies (Li et al., 2019). This study screened the highly variable sites of 42 Aconitum species by sliding window analysis. Seven hotspots with high variation were screened from the chloroplast genome sequence, including four intergenic spacer regions and three protein-coding gene regions (Figure 10B). They were distributed in the LSC region (trnK-UUU-trnQ-UGG, psbD, ndhJ-ndhK, clpP, and psbH-petB) and SSC region (ycf1 and trnA-UGC-trnI-GAU). The LSC region had the most variation sites, followed by the SSC region. In contrast, the IRb and IRa regions had relatively few variation sites. The distribution of Pi value was significantly lower than in LSC and SSC regions, which also revealed that in the chloroplast genome, the IRs region variation was conservative compared with LSC and SSC regions.

In addition, we performed chloroplast genome sliding window analyses of the ten Aconitum species included in this study (Figure 10A). We found that the variable sites were more divergent when variation analysis was performed on 42 Aconitum species. Analysis of the variant sites in all the 42 Aconitum species that have been sequenced allows screening of more suitable molecular identification markers for Aconitum species.

This study used Maximum Likelihood (ML) to construct phylogenetic trees based on two data sets (complete chloroplast genome and 84–85 protein-coding genes) from 44 species to determine the phylogenetic position of ten Aconitum species (Figures 11A,B). In addition, NJ tree with complete chloroplast genome sequences of 44 species were constructed. The topological structure of the NJ tree is essentially consistent with that of the ML tree (Supplementary Figure S1). The information on 44 species involved in phylogenetic tree construction is shown in Supplementary Table S6. The ML tree topology constructed based on the two data sets was essentially consistent, and a high support rate was detected on most nodes, but the support rate of the ML trees constructed based on different data sets was different (Figures 11A,B). The clustering of the two phylogenetic trees was the same, which strongly supports that Aconitum mainly consists of Subgen. Aconitum and Subgen. ParAconitum, with support values ≥100%. Two outgroup species (Delphinium grandiflorum and Delphinium anthriscifolium) were independent. Forty-two Aconitum species were clustered into a single large clade with support values ≥60%. Of these, A. stapfianum, A. weixiense, A. vilmorinianum, and A. episcopale were clustered together. A. stapfianum and A. episcopale exhibited a sister relationship with 96% and 100% support values, indicating the close relationship between the four Aconitum species. A. delavayi, A. ramulosum, and A. ouvrardianum clustered together, with 98% and 100% support values. A. delavayi and A. ramulosum exhibited a sister relationship. However, A. stylosum, A. stapfianum, A. weixiense, A. vilmorinianum, and A. episcopale were clustered together in ML phylogenetic tree based on PCGs. Our results will provide genomic resources for the phylogeny of Aconitum species and facilitate future phylogenetic studies and other studies of Ranunculaceae plants.