Introduction

The Indian traditional medicinal plant genus Aconitum encompasses ~ 300 mountainous species of significant economic importance, albeit some being poisonous due to the presence of toxic diterpene alkaloids1. These plants have application in Ayurveda as well as Chinese traditional medicine addressing several ailments like neuralgia, sciatica, arthritis, gout, rheumatism, treatment of colds, sore throat, and inflammation of the respiratory tract1,2. As of now, their unique biodiversity is on the verge of extinction due to illegal human intervention triggered habitat loss, over-harvesting, and unrestricted trading, and based on this, several Aconitum species have been identified as endangered, critically endangered or vulnerable by IUCN3.

Phylogenetic studies of Ranunculaceae family based on diverse characteristics such as morphology, restriction site mapping, nuclear sequence and chloroplast sequence have complicated the classification of genus Aconitum4,5. Several classification models of genus Aconitum have been proposed, based on different morphological characteristics such as inflorescence, branching of stem, shape of sepals and petals and structure of embryo sac among others6,7. According to the chromosomal study by Schafer & La Cour (1934), Aconitum genus is classified into two subgenus: Lycoctonum and Aconitum based on ploidy levels of the two groups8. Based on several non-molecular characteristics such as phytochemical, cytological, anatomical and palynological (study of plant pollen and spores), Aconitum gymnandrum has been removed from genus Aconitum and converted into separate genus Gymnaconitum, with only one species Gymnaconitum gymnandrum9. In this context, Gymnaconitum gymnandrum has been widely used as an outgroup for phylogenetic and various other analysis.

Taxonomic and phylogenetic studies in Aconitum have long been challenged by conflicting signals from morphology, hybridization, and limited availability of molecular markers. Several phylogenetic approaches have been explored, including analyses based on rbcL–matK and ITS sequences10,11,12,13; however, these markers failed to consistently cluster sequences from the same species, limiting their utility for species-level discrimination. Several classification systems have been proposed based on different morphological characteristics and later refined using molecular evidence, yet the evolutionary relationships among subgenera, sections, and series remain unresolved4. Consequently, there is still no consensus regarding the phylogenetic boundaries within Aconitum, and reliable molecular markers for species identification are lacking14. This raises key scientific questions: How variable are chloroplast genomes within and between Aconitum species? Can plastome features be used to clarify phylogenetic relationships and identify species-specific markers?

In this study, we performed a comparative chloroplast genome analysis of 73 Aconitum species, along with the outgroup Gymnaconitum gymnandrum. The dataset includes multiple accessions for several species, allowing us to assess not only inter-specific but also intra-specific plastome variation—an aspect largely overlooked in previous studies. The present work aims to: (i) examine genome structure, SSR distribution, codon usage, and nucleotide diversity across Aconitum; (ii) construct a pan-plastome to identify core and variable regions; and (iii) evaluate phylogenetic relationships using comprehensive chloroplast datasets without enforcing monophyly on conspecific samples.

Methodology

Data collection and summary

Chloroplast genome data of genus Aconitum members (including genome nucleotide sequences, translated CDS sequences, and GenBank files) was retrieved from the National Centre for Biotechnology Information (NCBI) via Entrez using command line interface. A total of 105 chloroplast genomes were obtained for which a summary file (Supplementary Table 1) was generated using a custom script to extract relevant metadata from the GenBank files, including species name, accession number, genome size, GC content, and the counts of protein-coding genes, rRNA, tRNA, genes, pseudogenes, and voucher information. To eliminate redundancy, genomes were screened based on species names and genome sizes resulting in a refined dataset comprising 73 unique Aconitum chloroplast genomes. Gymnaconitum gymnandrum, a species closely related to Aconitum but belonging to a different subgenus, was included as an outgroup for the analysis, consequently, making the final dataset consisting of 74 chloroplast genomes.

Genome annotation for differentiation into LSC, SSC and IR regions

PGA (Plastid Genome Annotator) tool15 was used to perform annotation of all 74 genomes, and categorization into LSC, SSC and IR regions. Using the annotation data, the sequences of four regions were extracted using a self-written bash script followed by calculation of GC content in each of the regions. The genome with accession ID MW817090 of A. scaposum was excluded from this analysis due to discrepancy in the annotation of four regions.

Pan-plastome analysis and functional characterization of chloroplast genes

Pan-plastome analysis was performed with the tool Proteinortho 6.1.716, where translated CDS files of all 74 Aconitum genomes was provided as the input. This software detects orthologous genes within different species, comparing similarities of given gene sequences and clusters them to find significant groups. The tool was run with default parameters (percentage identity = 25%, evalue = 1e-5) along with two additional parameters –singles and –selfblast which potentially help in detecting the unique genes in the genomes if present. Calculation of beta value for Heap’s law and generation of graph was performed using R.

Synteny analysis

Synteny studies were performed to visualize the genome architecture of 74 Aconitum genomes using web-based OrganellarGenomeDRAW v1.3.1 tool17 and Genbank (.gbk) files as input. The tool converts GenBank or EMBL/ENA format to graphical maps, either circular or linear. Two different set of linear maps were generated for each genome, one representing only the accessory genes and the other representing core and accessory genes. Configuration file was edited for each set to specify required genes in the map and their corresponding colour code. Each set resulted in 74 high-quality png files of linear maps which were merged into one and ‘convert’ command was used to crop, resize and merge the individual images.

Analysis of simple sequence repeats (SSRs)

GMATA v2.01 tool was used to detect SSRs18 using the genome file of 74 genomes. Self-written script was used to extract data from the ‘.ssr.sat2’ output file and graphs plotted using the extracted data. Minimum repeated times of motif was set to 10, 5, 4, 3, 3 and 3 for mono, di, tri, tetra, penta and hexa nucleotide repeats.

Variable region identification

To identify the variable region in Aconitum chloroplast genomes with respect to the reference, BLAST alignment was performed using Circular genome viewer comparison tool (CG View CT)19. This tool utilizes GenBank-formatted files as input and operates through a two-step process (project creation and map generation) facilitated by the wrapper script build_blast_atlas.sh. During the project creation step, a structured directory system is established, consisting of directories for input files, reference files, configuration files, and output files. In the subsequent map generation step, genome FASTA files are placed in the ‘comparison_genome’ folder, while the reference genome FASTA file (in this case, Gymnaconitum gymnandrum) is placed in the reference folder. The build_blast_atlas.sh script then generates maps for both nucleotide (BLASTn) and translated coding sequence (BLASTp) comparisons. The output includes CGView XML files located in the cgview_xml folder, corresponding to the nucleotide (dna_vs_dna) and protein-coding sequence (cds_vs_cds) comparisons. The same methodology was applied using A. vilmorinianum as the reference genome for comparative analysis.

Codon usage analysis

The total coding sequences for all 74 genomes (Table 1) were filtered according to the following criteria suggested by previous report20:

  • The sequence length must exceed 300 base pairs.

  • Each sequence must initiate with a start codon (ATG) and terminate with a stop codon (TAA/TAG/TGA).

  • Intermediate stop codons must be absent within the sequence.

  • The total number of nucleotides in the sequence must be divisible by three.

Table 1 Pan-plastome analysis based functional characterization of genes within the broad functional category and subcategory. The table also shows their respective gene and their gene product.
Full size table

The GC content of the first, second, and third codon positions of the 40 protein-coding sequences that passed filtration criteria was calculated using the CUSP program from the EMBOSS package21 along with the overall codon usage frequency for all 74 chloroplast genomes. Codon usage analysis was conducted employing the CAI Calculator, which provided insights into nucleotide composition, relative synonymous codon usage (RSCU) values, codon usage frequency, and codon usage per thousand values22 followed by heatmap generation using Heatmap Illustrator (HemI 1.0)23 using the average linkage method with Euclidean distance as the clustering metric. The effective number of codons (Nc) and the expected Nc values were calculated using the software DAMBE724.

To test context dependent mutation, amino acids are categorised into eight distinct fourfold degenerate families as follows: Arg, Leu and Ser were each divided into three twofold degenerate families (Arg2, Leu2, Ser2) and three fourfold degenerate families (Arg4, Leu4, Ser4) along with five fourfold degenerate families Pro4, Thr4, Ala4, Val4, Gly4. Considering mutation as an independent single-site event, the nucleotide frequencies of the third codon position in the fourfold degenerate families will not be affected by the second and/or the first codon positions. Following the method used in analysis of codon usage in Quercus chloroplast genome25, test of independence of the third codon position in the eight fourfold degenerate families was tested. The chi-square test of independence was conducted for each of the six dataset as follows: Leu4/Pro/Arg4; Val/Ala/Gly; (Leu4 + Val)/(Ser4 + Pro + Thr + Ala)/(Arg4 + Gly); Leu4/Val; Arg4/Gly; Ser4/Pro/Thr/Ala. These analyses aimed to determine whether the nucleotide composition at the third codon position within each dataset is statistically independent of the first and second codon positions.

Synonymous and non-synonymous substitution analysis

To analyse synonymous and non-synonymous substitutions, DnaSP v6.12.03 26 was employed using nucleotide FASTA files of 40 protein-coding sequences conserved across 74 species. The resulting output included pairwise estimates of Ka (non-synonymous substitution rate) and Ks (synonymous substitution rate), along with additional metrics. For each gene, the mean Ka and Ks values were subsequently calculated. To further understand overall nucleotide diversity, DnaSP v6.12.03 was further used26 where fasta sequences of coding genes, rRNA, tRNA genes and intergenic sequences were provided as input all at once using the Batch Mode and default parameters. Pi values for each sequence were used to generate a line graph using R.

Phylogenetic analysis

The genus Aconitum has been classified into two subgenera, Aconitum and Lycoctonum, based on differences in morphology and ploidy levels. Therefore, to determine whether this morphological and ploidy level classification is supported by molecular data, two distinct phylogenetic analysis was conducted using core gene sequences and whole-genome data. For core genes-based phylogeny, 74 fasta files containing protein sequences of each core gene per genome was prepared using a self-written script with sequences extracted from translated CDS file. The extracted sequences were aligned using MUSCLE v3.8.155127 followed by concatenating all 74 aligned blocks, which was further used for phylogenetic tree building using IQ-TREE multicore version 2.2.2.328. The best model suggested by the tool is used to create the phylogenetic tree followed by visualization on the web-based tool iTOL v629. For whole genome-based phylogeny, 74 whole genome sequences were aligned using the tool MUSCLE v3.8.155127. The aligned sequences were given as input to IQ-TREE (multicore version 2.2.2.3) to generate a maximum likelihood phylogeny using the model suggested by IQ-TREE multicore version 2.2.2.328. The generated tree is then visualized on the web-based tool iTOL v629.

Calculation of intra-specific K2P distance

The Kimura 2-Parameter30 substitution model was used to find the intra-specific distance in R with the help of libraries named “ape” and “seqinr”. Input for calculating the distance was MUSCLE alignment file of all 73 genomes along with the outgroup as mentioned in the phylogenetic analysis methodology.

Result

Chloroplast genome statistics of 74 Aconitum species

Out of the 73 Aconitum genomes under study, the presence of 40 unique species depicts the diversity of the used dataset (Supplementary Table 1). Out of the total available data, 30 and 43 genome files were sourced from RefSeq and GenBank databases, respectively. The genome size of these species ranges from 151,214 to 157,688 bp with A. episcopale having the smallest and A. brachypodum having the largest genome. Variation in number of genes is not very high and ranges from 123 in A. austrokoreense, A. coreanum and A. volubile to 132 in twenty different Aconitum genomes (Supplementary Table 1). Interestingly, one out of two reported genomes A. coreanum and A. austrokoreense has 123 genes annotated and the other has 132. The protein coding genes range from 82 to 87, and tRNA genes ranges from 36 to 38, however, 8 rRNA genes are present consistently across all genomes under study. Six pseudogenes are also annotated in A. pseudolaeve, which is the highest among all the considered genomes. All 73 genomes of Aconitum had a quadripartite structure, i.e. had a large single copy region (LSC), small single copy region (SSC) and two inverted repeat regions (IRA and IRB) (Fig. 1a). As far as the total GC content is considered, there is no significant difference in GC percentage among the chloroplast genomes, which ranges from 37.99% to 38.30%. However, the GC content between different regions of the same genome varies, with IR regions having high GC percentage and LSC region having lower GC percentage (Fig. 1b, Supplementary Table 2). The GC percentage of LSC region ranges from 36.01 to 36.39%, that of SSC ranges from 32.42% to 32.85%, whereas for IRA and IRB region the range lies from 42.94% to 43.10%. This indicates that although there is difference amongst the four regions, there is no significant difference among the different species for one particular region.

Fig. 1
Fig. 1
Full size image

(a) Circular map of A. barbatum chloroplast genome. Genes represented on the outer side of the outermost track are transcribed clockwise and genes represented on the inner side are transcribed anti-clockwise. The second track represents GC content of chloroplast genome at different loci. LSC: Large Single Copy region; SSC: Small Single Copy region; IR: Inverted repeat region. (b) Jitter and box plot of GC content of four regions of chloroplast genome as well as the overall GC content of the 74 species. (c) Heap’s law model graph indicates that the pan-plastome of Aconitum chloroplast is open. The Heap’s law formula is given by V = k*n^beta; beta > 0 indicates open pan-plastome and beta < 0 indicates closed pan-plastome. Here, beta value was 0.0084.

Pan-plastome analysis revealed that Aconitum has an open pan-genome

Pan-plastome analysis revealed that the number of core genes ranged from 75 to 77 among the 74 Aconitum species, whereas the accessory genes ranged from 6 to 10 (Supplementary Table 3). To evaluate the openness of the chloroplast pan-plastome, Heap’s law formula (V = k·n^β) was applied, where V represents the number of distinct genes, n is the number of genomes, k the scaling parameter, and β the growth parameter. The β value was calculated to be 0.0084, which, being greater than zero, indicates that the Aconitum chloroplast pan-genome is open (Fig. 1c). This suggests that as additional chloroplast genomes of Aconitum are sequenced, novel or variable genes are still likely to be discovered. The gene order of both accessory and core genes was found to be well conserved across all species (Fig. 2, Supplementary Fig. 1). Functional characterization of genes (Table 1) classified all gene functions into three main categories: chloroplast envelope membrane protein genes, genes for photosynthesis and genes for transcription and translation, which were further classified into sub-categories. The Chloroplast envelope membrane protein genes include the Cytochrome b6f. group of genes, whereas the photosynthesis genes include ATP synthase (atp genes), NADH oxidoreductase (ndh genes), photosystem I (psa genes) and photosystem II genes (psb genes). The sub-categories of genes for transcription and translation include the large ribosomal subunit (rpl genes), RNA polymerase (rpo genes), small ribosomal subunit (rps genes) and translation initiation factor (infA gene). All photosynthetic genes were assigned under core genes. Five (rpl16, rpl2, rpl20, rps16 and infA) out of nine accessory genes belonged to the category of genes related to transcription and translation. A group of researchers have performed the functional and structural analysis of rpl16 using bioinformatics tools, suggesting it to be a thermo-stable, acidic and hydrophilic protein. One of the predicted counterparts of RPL16 includes RPL2 which is also an accessory protein among the genomes considered in the study. The presence of its counterpart rpl2 in A. reclinatum (MF186593) could be a possible reason behind the absence of rpl16 gene having no effect. The protein RPL20 along with other proteins initiates the 50 s ribosomal subunit assembly which binds directly to the 5’ end of the 23 s rRNA31. It has been found that rps16 is lost in many taxa from ferns to angiosperms. However, its nuclear genome counterpart is present in such taxa as it is necessary for the survival of the organism32. The investigation of nuclear genome will reveal presence or absence of such counterparts in the Aconitum genome.

Fig. 2
Fig. 2
Full size image

Accessory gene map representing synteny across 73 genomes of Aconitum genus and the outgroup Gymnaconitum gymnandrum.

Analysis of SSRs exhibits lack of their conservation in Aconitum chloroplast genomes

This analysis revealed that mononucleotide SSRs exhibited the highest frequency, followed by di-, tri-, tetra-, penta-, and hexanucleotide repeats (Fig. 3b–f, Supplementary Table 4a). The frequency of mononucleotide repeats ranged from 0.12 to 0.33 SSRs/kb, with A. finetianum displaying the highest abundance of mononucleotide repeats (Fig. 3b). The number of mononucleotide SSRs varied from 19 in A. volubile to 52 in A. finetianum. In contrast, the number of dinucleotide SSRs showed minimal variation, ranging from 10 to 17 (Fig. 3a), with a frequency distribution of 0.06 to 0.1 SSRs/kb. The frequency of trinucleotide, tetranucleotide, and pentanucleotide repeats ranged from 5 to 12, 5 to 9, and 1 to 5, respectively, with pentanucleotide SSRs absent in certain species. Hexanucleotide SSRs were identified exclusively in eight species: A. jaluense subsp. jaluense, A. austrokoreense, A. scaposum var. vaginatum, A. longecassidatum, A. tanguticum, A. ramulosum, A. stylosum, and A. delavayi. Variability in di-, tri-, and tetranucleotide SSRs was relatively low across species (Fig. 3a). The frequency distributions of tri-, tetra-, penta-, and hexanucleotide repeats were 31.78–76.26, 31.71–57.8, 0–31.8, and 0–19.09 SSRs/Mb, respectively (Supplementary Table 4a).

Fig. 3
Fig. 3
Full size image

SSR analysis for 74 Aconitum genomes showcasing the lack of conservation of SSRs in the genomes. (a) Number of SSR loci per genome: here X-axis and Y-axis depicts accession ID of the genomes and number of SSRs, respectively. Region wise frequency (SSRs/kb) of (b) mononucleotide repeats (c) dinucleotide repeats (d) trinucleotide repeats (e) tetranucleotide repeats (f) pentanucleotide repeats. It should be noted that the scale for the graphs varies according to the range of frequencies for different types of SSR.

SSR analysis across the LSC, SSC, IRA, and IRB regions revealed that the LSC, followed by the SSC and IR regions, harboured the highest number of all six types of repeat sequences, which is expected given the larger genomic span of the LSC. Two genomes, MW817090 and MT584425, were excluded from this analysis due to discrepancies in IR region annotations. The frequency of total SSRs, expressed in SSRs/kb, was highest in the LSC region. Specifically, mononucleotide repeat frequencies ranged from 1.7 to 4.6, 0 to 3.8, and 0 to 1.2 SSRs/kb in the LSC, SSC, and IR regions, respectively (Fig. 3b). Dinucleotide repeat frequencies varied between 0.3 and 0.9 SSRs/kb in the LSC, 0 to 0.8 SSRs/kb in the SSC, and 0 to 0.4 SSRs/kb in the IR regions (Fig. 3c). Trinucleotide repeats exhibited frequency distributions of 0.04 to 0.3 SSRs/kb in the LSC, 0.2 to 0.6 SSRs/kb in the SSC, and 0 to 0.15 SSRs/kb in the IR regions (Fig. 3d). Tetranucleotide repeats were identified in only one genome within the IR regions at a frequency of 0.1 SSRs/kb, whereas in the LSC and SSC regions, their frequencies ranged from 0.1 to 0.2 and 0.2 to 0.5 SSRs/kb, respectively (Fig. 3e). Pentanucleotide repeats were observed at frequencies of 0 to 0.14 SSRs/kb in the LSC region and were detected in the SSC region in only one genome at a frequency of 0.2 SSRs/kb, while they were completely absent in the IR regions (Fig. 3f). Hexanucleotide repeats were identified in six genomes within the LSC region, all exhibiting a frequency of approximately 0.03 SSRs/kb. These repeats were absent in the SSC region, whereas in the IR regions, they were found in only one genome at a frequency of 0.1 SSRs/kb (Supplementary Table 4b).

Among mononucleotide repeats, A/T-rich repeats were more prevalent than C/G-rich repeats in both the LSC and SSC regions, whereas the IR regions contained only A/T repeats. The only dinucleotide repeats present across all four regions were of the AT/AT type. The LSC region contained three distinct trinucleotide repeat motifs: AAT/ATT, ATC/ATG, and CCG/CGG, while the SSC region exhibited only the AAT/ATT motif, and the IR regions exclusively harboured AAG/CTT trinucleotide repeats. The dominant tetranucleotide repeats in the LSC region were AAAG/CTTT and AAAT/ATTT, whereas AATG/ATTC was the predominant motif in the SSC region (Supplementary Table 4c).

Intra-species distribution of simple sequence repeats (SSRs)

To assess intra-specific variation in repetitive elements, we compared the distribution of chloroplast SSRs among multiple accessions belonging to the same Aconitum species. Overall, the SSR distribution was highly conserved within several species, particularly across di-, tri-, and tetranucleotide repeat classes. For instance, A. barbatum (including A. barbatum var. hispidum and A. barbatum var. puberulum) exhibited almost identical frequencies across all SSR types, differing only slightly in mononucleotide repeats (248.8–267.9 SSRs/Mb). Similarly, A. pseudolaeve accessions showed near-identical SSR profiles, with differences of less than 0.01 SSRs/Mb across di-, tri-, and tetranucleotide repeats, and variation restricted to mononucleotides. Consistent patterns were also observed in A. coreanum, A. carmichaelii, A. kusnezoffii, and A. scaposum, suggesting a high level of plastome stability within species.

Moderate intra-specific divergence was detected in A. brachypodum, A. pendulum, A. flavum, and A. episcopale, where specific accessions deviated notably in mononucleotide or pentanucleotide frequencies, occasionally showing the gain or loss of hexanucleotide motifs. In A. delavayi, one accession possessed hexanucleotides while another lacked them, indicating subtle but meaningful plastome variability. These findings collectively demonstrate that chloroplast SSR profiles are largely uniform within species but can reveal cryptic differences among conspecific lineages.

Variable region identification and nucleotide diversity analysis reveal some tRNA genes with diverse sequence

Whole-genome BLASTN analysis was conducted for 74 chloroplast genomes, and comparative circular plots were generated using the CGView Comparison Tool against Gymnaconitum gymnandrum and A. vilmorinianum as reference genomes. The resulting plots illustrate the percentage similarity across genomes. Overall, the chloroplast genomes of the 74 Aconitum species exhibit high conservation, with more than 90% sequence identity. In the first plot, several regions display sequence identity below 96%, with a small stretch between matK and psbI showing less than 90% identity (Fig. 4a). The second plot indicates variability in the matK to psbI region in only a subset of the genomes analysed (Fig. 4b). Further examination revealed that these genomes belong to the subgenus Lycoctonum, whereas the reference genome used for the second plot belongs to the subgenus Aconitum. To gain deeper insights into genome-wide sequence variability at the nucleotide level, a nucleotide diversity analysis was performed.

Fig. 4
Fig. 4
Full size image

CG View plot of 73 Aconitum genomes based on BLASTn analysis against (a) G. gymnandrum as reference (b) A. vilmorinianum as reference, showcase sequence identity across all studied genomes.

Nucleotide diversity (Pi) analysis provides a measure of sequence variation across different genomic regions. The Pi value represents the proportion of nucleotide sites expected to differ between any two randomly selected DNA sequences, with higher values indicating greater sequence variability. Regions exhibiting high Pi values are potential candidates for marker development. This analysis illustrates the Pi values of various genes (Fig. 5a) and intergenic regions (Fig. 5b), arranged in genome order. Notably, sequences with Pi values exceeding 0.2 predominantly correspond to tRNA genes, along with a single intergenic region between psbH and petB. The tRNA genes exhibiting high nucleotide diversity include trnL-CAA, trnN-GUU, trnV-GAC, trnR-ACG, trnI-GAU, trnA-UGC, and trnI-CAU, all of which are located within the IR regions. Additionally, in non-IR regions, tRNA genes such as trnG-GCC, trnG-UCC, and trnM-CAU show elevated Pi values (greater than 0.1). These highly diverse tRNA gene sequences hold potential for marker-based species identification. Within the matK–psbI region, rps16 is the only gene exhibiting relatively higher nucleotide diversity, with a Pi value exceeding 0.05.

Fig. 5
Fig. 5
Full size image

Graph of Pi values representing nucleotide diversity of (a) genes (b) intergenic regions, depicting extent of variation in respective sequences amongst the 74 species.

Codon usage analysis

Codon usage bias, the preferential use of certain synonymous codons over others, reveals fundamental evolutionary forces shaping genomic architecture. There are several aspects of codon usage that provide insights into crucial aspects of molecular evolution such as GC content of the codons, codon usage bias measured by Nc values (Effective number of codons) and codon usage preference based on RSCU values. Previous studies on chloroplast genomes have consistently shown that GC content decreases from the first (GC1) to the third (GC3) codon positions, favouring A/T-ending codons due to mutational pressures and compositional bias25,33,34, whereas chloroplast genomes shows consistent ENc values across species, reflecting weak overall codon bias and evolutionary conservation33,34,35. Relative Synonymous Codon Usage (RSCU) analysis has identified certain amino acids absent in specific genes across chloroplast genomes, emphasizing lineage-specific translational optimization and evolutionary constraints25,36. It has also been reported in previous studies through understanding the relationship between Nc and GC3 values that GC3 exerts some influence on the codon usage pattern although other factors like selection also play important role37,38.

GC content reduces from first position of codon to the third indicating preference of A/T ending codons over G/C ending codons

The base composition at the first (GC1), second (GC2), and third (GC3) codon positions was analysed for all 40 genes across 74 chloroplast genomes (Supplementary Table 5a). While GC content varied among genes, no significant variation was observed between genomes (Supplementary Fig. 3). However, a significant difference was noted among GC1, GC2, and GC3 values (Fig. 6a, Supplementary Fig. 3). Across genomes, GC1 values ranged from 47.19% (n = 74, SD = 5.79) to 47.49% (n = 74, SD = 5.65), GC2 values from 39.27% (n = 74, SD = 5.24) to 39.73% (n = 74, SD = 5.36), and GC3 values from 28.58% (n = 74, SD = 3.39) to 28.90% (n = 74, SD = 4.02). In contrast, variation across genes was more pronounced, with GC1 values ranging from 36.14% (n = 40, SD = 0.24) to 58.55% (n = 40, SD = 0.14), GC2 values from 27.94% (n = 40, SD = 0.21) to 57.55% (n = 40, SD = 0.00), and GC3 values from 23.07% (n = 40, SD = 0.52) to 35.92% (n = 40, SD = 0.41).

Fig. 6
Fig. 6
Full size image

(a) Box plot of GC content of first, second and third position of codons (b) Gene-wise heatmap of Nc values (effective number of codons) and (c) RSCU values.

The heatmap presented in Supplementary Fig. 3 highlights the lower GC3 values, indicating a preference for A/T-ending codons over G/C-ending codons. This pattern has been commonly observed and reported across chloroplast genomes of various species, including Morus cathayana, Morus multicaulis, six Euphorbiaceae species, and three Camellia species, among others39,40,41.

Codon usage bias measured by effective number of codons reveals consistency of codon usage across genomes

The effective number of codons (Nc) metric was utilized to assess codon usage bias among different genes and across all 74 analysed genomes. Nc values range from 20 to 63, where lower values indicate a stronger codon bias, and higher values suggest a more uniform usage of synonymous codons. A lower Nc value implies that an organism preferentially utilizes a subset of synonymous codons, whereas a higher Nc value reflects a reduced bias in codon selection.

As illustrated in Fig. 6b, the rps18 and petD genes exhibit the highest codon bias among the 40 analyzed genes, whereas ycf4 and clpP display the least bias, clustering together in the heatmap. Additionally, genes such as rps14, psbA, ndhA, petB, and atpF exhibit a codon usage pattern similar to petD, forming a distinct cluster. The heatmap further highlights variations in codon bias among genes while demonstrating a largely consistent trend across the genomes under study.

Codon usage preference based on RSCU values reveals the relative absence of some amino acids in several genes across 74 chloroplast genomes

Comparative analyses of transfer RNA (tRNA) across all kingdoms have demonstrated that no single organism possesses tRNAs with anticodons complementary to all 61 sense codons42,43. This is due to the wobble hypothesis, where a single tRNA can recognize multiple synonymous codons because the third position of the codon (wobble position) exhibits flexibility in base pairing.

Relative Synonymous Codon Usage (RSCU) quantifies codon usage bias by comparing the observed frequency of synonymous codons for a given amino acid to the expected frequency under equal usage conditions. An RSCU value of 1 indicates no bias, values greater than 1 signify positive codon usage bias, while codons with RSCU values below 0.6 are considered underrepresented, and those above 1.6 are overrepresented.

RSCU analysis revealed that atpE and rpoC2 genes lack both codons for tyrosine, rps18 lacks both codons for histidine, psbA lacks both codons for lysine, and ndhC, rps18, and rps7 lack both codons for cysteine. As depicted in Fig. 6c, G/C-ending and A/T-ending codons tend to form distinct clusters, with A/T-ending codons exhibiting relatively higher RSCU values. Notably, most G/C-ending codons clustered together, along with TTT, ATA, and CTA, whereas within the major A/T-ending cluster, only TTG was grouped with the G/C-ending codons.

Codons were classified into six groups based on their RSCU values: (1) overrepresented codons (RSCU > 1.6), (2) positively biased codons but not overrepresented (1 < RSCU ≤ 1.6), (3) unbiased codons (RSCU = 1), (4) negatively biased codons but not underrepresented (0.6 ≤ RSCU < 1), (5) underrepresented codons (0 < RSCU < 0.6), and (6) unused codons (RSCU = 0). Our analysis revealed that 80.16% of thymine-ending codons and 71.25% of adenine-ending codons were positively biased, whereas only 12.19% of cytosine-ending and 14.23% of guanine-ending codons exhibited positive bias (Supplementary Table 6).

A comparison of the codon usage pattern of the psbA gene with the overall trends observed across 40 genes indicated that psbA follows a similar codon usage preference to the whole genome. However, fewer than half of the adenine-ending codons were positively biased (Supplementary Table 6), suggesting that the psbA gene exhibits a stronger preference for A/T-ending codons.

GC3 values have low level correlation with the effective number of codons suggesting some influence of GC3 on codon usage

Correlation analysis revealed that neither GC1 nor GC2 values exhibited a significant correlation with the effective number of codons (Nc) (Fig. 7a–c). Scatter plot of correlation analysis depicting the comparison between expected and observed Nc values against GC3 demonstrates a noticeable deviation of observed Nc values from expected trends (Fig. 7c). The correlation coefficient was higher for expected Nc values against GC3 than for observed Nc values, indicating that while GC3 strongly influences theoretical Nc expectations, actual codon usage patterns deviate due to additional evolutionary factors.

Fig. 7
Fig. 7
Full size image

Scatter plot of (a) Nc and ENc values against GC1; R2 = 0.06 and 0.14 (b) Nc and ENc values against GC2; R2 = 0.0 and 0.0009 (c) Nc and ENc values against GC3; R2 = 0.33 and 0.83 (d) GC1 and GC2 against GC3; R2 = 0.02; Nc is effective number of codons; ENc is expected effective number of codons.

A low positive correlation was observed between Nc and GC3 values, with observed Nc values tending to be lower than expected at higher GC3 values. This suggests that an increased proportion of GC-ending codons corresponds to reduced codon bias. Notably, the ndhA gene clustered with psbA in Fig. 6b, and its observed Nc value (50.2) closely matched the expected value (50.5) (Supplementary Table 5d). However, this trend did not hold true for psbA, implying that distinct selective pressures or mutational influences are shaping the codon usage of psbA differently from other genes.

Additionally, regression analysis of average GC1 and GC2 values against GC3 yielded an insignificant R2 value (0.02), reinforcing the notion that codon usage patterns arise from a complex interplay of mutational and selective forces (Fig. 7d).

Analysis of RSCU values indicates absence of context dependent mutation

The analysis of context dependent mutation verifies the condition that mutation is a single-site event, meaning the nucleotide frequencies of the third codon position in the fourfold degenerate families will not be affected by the second and/or the first codon positions. The result of this analysis showed that the variation of codon’s third base does not correlate with either second or first base (Supplementary Table 7). It can thus be hypothesised that for the Aconitum chloroplast genomes, the likelihood of mutation occurring at a particular site in the codon is not influenced by the neighbouring nucleotides.

Synonymous and non-synonymous substitution analysis reveal that the Aconitum chloroplast genes are under purifying selection

A nucleotide substitution that alters the encoded amino acid of a protein is termed a non-synonymous substitution (Ka), whereas a substitution that does not change the amino acid sequence is referred to as a synonymous substitution (Ks). The Ka/Ks ratio serves as an indicator of coding sequence evolution, where a value of 1 suggests neutral evolution, a value greater than 1 indicates positive or diversifying selection, and a value less than 1 implies negative or purifying selection. The Ka/Ks ratio was calculated for 40 conserved genes across 74 Aconitum species, revealing a predominant pattern of purifying selection. In this analysis, Ka values were generally lower than Ks values, with the exceptions of matK, clpP, and rpoC1 (Supplementary Fig. 4, Supplementary Table 8). This finding aligns with expectations, as all analysed genomes belong to the same genus, where strong evolutionary constraints act to preserve functional integrity and limit protein-coding sequence divergence.

Non-monophyletic clustering of samples in phylogenetic analysis

To assess whether phylogenetic analysis aligns with the morphological classification, whole-genome and core-gene phylogenies were constructed using the best-fit models suggested by IQ-TREE. The whole-genome phylogeny was inferred using the TVM + F + I + R3 model, while the core-gene phylogeny was based on the Q.mammal + F + R2 model. In both phylogenies (Fig. 8) all the accessions of subgenera Aconitum and Lycoctonum cluster together respectively, with one exception. A. flavum (GenBank accession: MT982388) was observed to cluster phylogenetically with members of the subgenus Lycoctonum, rather than grouping with its taxonomic subgenus Aconitum, as would be expected based on current classification. This anomalous placement was consistently recovered in phylogenies constructed from both whole chloroplast genome sequences and core gene datasets of 73 species (Fig. 8).

Fig. 8
Fig. 8
Full size image

Core-gene phylogeny (left) and Whole genome phylogeny (right) for 73 Aconitum genomes with G. gymnandrum as an outgroup. Both maximum likelihood phylogenies have been generated using IQTREE using its best model selection pipeline. ITol was used for visualization with subgenus mapped to the phylogenies for better understanding.

To further investigate this incongruence, Kimura 2-Parameter (K2P) genetic distances were calculated among A. flavum accessions. The pairwise comparisons revealed that MT982388 displayed markedly higher intraspecific K2P distances (ranging from ~ 0.0093 to 0.0095) when compared with other A. flavum accessions (e.g., MW839579, MW839580, MW839582, and NC_056280), whose mutual distances were significantly lower (as low as 4.5 × 10⁻5) (Supplementary Table 10).

The pattern of monophyletic clustering, although observed at the level of subgenus, was not consistently observed at the levels of series or species. Several species exhibited strong intraspecific monophyly across both the whole-genome and core-gene trees, including A. scaposum, A. barbatum, A. coreanum, and partially A. kusnezoffii. Likewise, consistent series-level monophyly was recovered for Bullatifolia, Brachypoda (partially), Scaposa, and Longecassidata. In contrast, intraspecific non-monophyly was observed in A. sinomontanum, A. carmichaelii, A. delavayi, A. vilmorinianum, and A. austrokoreense, while the series Ambigua, Stylosa, Volubilia, Inflata, and Racemulosa displayed polyphyletic or paraphyletic clustering patterns, suggesting possible taxonomic inconsistencies or historical hybridization events.

Certain accessions showed anomalous clustering patterns, such as A. kusnezoffii (NC_031422), which did not group with other conspecific samples, implying potential misidentification or chloroplast capture. Similarly, A. delavayi (NC_038097) did not cluster with other members of its designated series (Ambigua), and instead grouped distantly, indicating either sample mislabelling or complex evolutionary history within this lineage.

Discussion

Chloroplast genomes typically display a high degree of conservation in structure and gene content, a pattern also evident in the comparative analysis of 73 Aconitum chloroplast genomes. Despite this overall conservation, the pan-plastome analysis revealed an open chloroplast pangenome, as indicated by the positive β value (0.0084) obtained from Heap’s law modelling. Rather than indicating the acquisition of truly novel genes, this openness likely reflects the annotation and structural variation among plastomes in the dataset. The pattern may arise from factors such as incomplete or inconsistent gene annotations across public genomes, lineage-specific gene truncation or pseudogenization, sampling effects due to unequal taxonomic representation, or structural rearrangements near the IR boundaries. Thus, the observed “open” configuration should be interpreted as evidence of residual annotation and structural variability within an otherwise highly conserved plastome rather than ongoing gene content expansion.

Interestingly, while nine genes were identified as accessory, BLASTn comparisons showed that homologous sequences of these genes were present across all Aconitum chloroplast genomes with > 96% similarity (Supplementary Table 9). Their apparent absence in several NCBI-annotated genomes points to annotation discrepancies rather than genuine gene loss, underscoring the importance of cross-validation using alternative annotation pipelines such as Plastid Genome Annotator (PGA)15. In the PGA-annotated datasets, the accessory genes infA and ycf15 still remained undetected in all genomes; however, other accessory genes, i.e., rpl16, rpl2, rpl20, rps16, psbN and ycf1, were detected. The observed inconsistencies between NCBI and PGA outputs emphasize the necessity of manual curation and improved automated algorithms for reliable gene identification in chloroplast genome studies.

Collectively, these findings not only highlight the open nature of the chloroplast pan-genome in Aconitum but also call attention to the methodological challenges in accurate plastome annotation, both of which have significant implications for understanding chloroplast genome evolution in the Ranunculaceae family. The observed inconsistencies between the NCBI and PGA annotations in the relative presence/absence of genes in genome sequences despite their absence in annotation files, may stem from undetected mutations or sequencing errors in chloroplast genome assembly. These results strongly emphasize the need for rigorous manual validation and improved annotation algorithms to ensure thorough gene identification in chloroplast genome studies.

The relative presence of several gene groups on chloroplast genome was assessed using KEGG pathway database. This study revealed that several members of different gene groups are absent from the chloroplast genome. It was noted that out of the 14 genes of Ndh group of oxidative phosphorylation, 11 are present on the chloroplast, the pet group of genes coding for cytochrome b6f. have 6 out of 8 members on the chloroplast whereas none of the Pet genes coding for proteins of a part of the photosynthetic electron transport chain, are present on the chloroplast genome. Further, the psa and psb group of genes coding for Photosystem I and Photosystem II respectively, have 5 out of 18 and 15 out of 28 genes respectively present on the chloroplast. As organellar (chloroplast) genomes have extremely reduced their genome size, the genes absent on the chloroplast genome might be encoded in the nuclear genome. The missing genes were looked up on the gene annotation file of A. thaliana nuclear genome. It was found that ndhL of the ndh group of genes, psaD, psaF, psaG, psaK and psaO of the psa group of genes, psbP, psbQ, psbR, psbY and psb27 of the psb group of genes are present in the nuclear genome.

SSRs are well known as microsatellites and range from 1 to 6 nucleotides as repeating units. The frequency of mononucleotide repeats ranged from 0.7 to 2%, indicating a significant variability among the species. Notably, A. finetianum exhibited the highest number of mononucleotide repeats (52), while A. volubile had the lowest (19). This considerable range within the same genus suggests that certain species within Aconitum may have undergone different evolutionary pressures or replication dynamics that influenced their SSR accumulation. Both whole genome analysis of SSRs and region wise analysis points towards the fact that SSRs do not exhibit conservation in the chloroplast genomes. The distribution of repeats in the four regions of genome namely LSC, SSC and IR regions is inconsistent among the 74 genomes and in some case, inconsistencies are observed even among different strains of the same species. The inclusion of multiple accessions per species in this study enabled detection of such intra-specific variation, which would have remained unrecognized in single-sample analyses. Previous studies on chloroplast SSRs in Aconitum14,44 were based on a limited number of species and typically included only one representative genome per species, thereby restricting their ability to assess intra-specific variation. The observed uniformity of SSR patterns within most species such as A. scaposum, A. coreanum, A. kusnezoffii, A. barbatum and A. episcopale (Supplementary Fig. 5) suggests that SSR distribution may serve as a potential diagnostic feature for species identification in Aconitum. Expanding the number of sampled individuals in future studies could further clarify the extent of intra-specific plastome variation and improve the resolution of SSR-based molecular markers for taxonomic and evolutionary research.

A previous study on codon usage in Aconitum chloroplast genomes45 primarily focused on genome-wide codon usage patterns across multiple species and reported a general AT bias at the third codon position (GC3) and low codon bias, consistent with findings across other angiosperms. However, their analysis was conducted at the whole-genome level, averaging codon frequencies across all coding regions. In contrast, our study performed a gene-wise codon usage analysis, enabling the detection of gene-specific deviations that are otherwise masked in genome-wide summaries. This approach revealed distinct codon bias patterns among genes—for example, psbA exhibited stronger selection-driven codon bias than mutation-driven bias, whereas ndhA followed expected mutational trends. Additionally, while the previous study did not assess the relationship between codon bias metrics and compositional parameters at the level of individual genes, our results demonstrated low but measurable correlation between GC3 content and Nc values, highlighting the combined influence of mutational bias and selection. Furthermore, our RSCU-based analysis identified gene-specific absence of certain amino acid codons and confirmed the independence of codon position context, providing deeper insight into translation-level constraints and evolutionary pressures acting on the Aconitum plastome. Thus, this study complements and extends earlier genome-wide observations by highlighting functional gene-level nuances in codon usage evolution within Aconitum chloroplast genomes.

Nucleotide diversity (Pi) analysis across 74 Aconitum chloroplast genomes revealed that the most variable regions are predominantly tRNA genes located within the IR regions, including trnL-CAA, trnN-GUU, trnV-GAC, trnR-ACG, trnI-GAU, trnA-UGC, and trnI-CAU. A few tRNA genes outside the IR, such as trnG-GCC, trnG-UCC, and trnM-CAU, also exhibited elevated diversity, alongside the rps16 gene within the matK–psbI region. These loci, particularly the tRNA genes, represent promising candidates for species identification within this taxonomically complex genus. The unexpectedly high nucleotide diversity observed within the IR region, primarily across tRNA genes, contrasts with previous studies in Aconitum46,47, which reported this region as the most conserved. This discrepancy may stem from alignment artefacts or annotation inconsistencies across publicly available genomes, particularly at the IR boundaries, which are known to vary48. Given that even a few divergent or misassembled IR sequences can substantially elevate Pi values in a dataset of 74 genomes, the observed pattern warrants cautious interpretation. Alternatively, minor IR expansion/contraction events or pseudogenization could also contribute to genuine increases in variability. Future reannotation and validation of IR boundary regions across Aconitum plastomes would help clarify whether the observed diversity reflects biological signal or assembly noise. Although multiple variable genes and intergenic regions were identified based on nucleotide diversity, single-gene phylogenetic analyses of these loci failed to consistently resolve species boundaries across the 73 Aconitum taxa. This result highlights an inherent limitation of chloroplast fragments for species-level identification in Aconitum, likely due to recent divergence, hybridization, chloroplast capture, and incomplete lineage sorting.

The analysis of synonymous (Ks) and non-synonymous (Ka) substitutions across 40 conserved chloroplast genes from 74 Aconitum species revealed that nearly all genes are under strong purifying selection, as indicated by Ka/Ks ratios significantly below 1. Only three genes—rpoC1, clpP, and matK—exhibited Ka values marginally higher than Ks, suggesting possible localized relaxation of selection or adaptive evolution in these loci. Among them, matK is known for its high substitution rates across angiosperms49 and encodes a maturase enzyme involved in group II intron splicing50,51. The rpoC1 gene, encoding a β subunit of chloroplast RNA polymerase, and clpP, encoding a proteolytic subunit of the ATP-dependent protease complex52,53, also displayed elevated Ka values, indicating that these genes may experience relatively relaxed evolutionary constraints. Our findings are broadly consistent with earlier reports by Yanfei et al. (2023) and Zhu et al. (2025), both of which identified purifying selection as the dominant evolutionary force acting on Aconitum chloroplast genes. Yanfei et al. (2023) observed Ka/Ks ratios below 1 for most genes across seven Aconitum species, except ycf1, rpl20, cemA, and rps18, which exhibited signs of positive selection in specific taxa47. Similarly, Zhu et al. (2025) reported Ka/Ks ratios < 0.5 for most protein-coding genes in Aconitum carmichaelii cultivars, with ycf1 showing relatively higher ratios, suggesting functional diversification and potential as a phylogenetic marker54. In contrast, our large-scale multi-species comparison did not detect ycf1 under positive selection, possibly due to broader taxon sampling across subgenera that averaged out lineage-specific adaptive effects observed in narrower datasets.

Together, these studies highlight the conserved nature of Aconitum chloroplast genes, where purifying selection maintains essential photosynthetic and housekeeping functions, while a few genes such as matK, rpoC1, and clpP show signatures of relaxed or positive selection. Such genes may serve as useful candidates for investigating adaptive divergence and resolving phylogenetic relationships within Aconitum.

The phylogenetic analysis validated seed morphology-based classification of subgenus Aconitum and Lycoctonum with one exception. The anomalous placement of A. flavum (GenBank accession: MT982388) was investigated using Kimura 2-Parameter (K2P) revealing divergence of the accession with its other conspecific members. The elevated divergence, along with the unexpected phylogenetic placement, suggests that MT982388 may represent a misidentified sample. Alternatively, it could reflect deep genetic divergence within A. flavum, indicative of cryptic speciation or historic hybridization. However, given that the GenBank submission includes a voucher specimen, it is more likely that the anomaly stems from either incorrect species identification during sequencing or contamination of the submitted sample. Hence, re-examining the original voucher specimens and corresponding sequencing data is recommended to confirm their identity and ensure the accuracy of chloroplast genome records. Further, the phylogenetic tree constructed from complete chloroplast genome sequences of 73 accessions representing 40 Aconitum species revealed non-monophyly among several conspecific samples. While multiple accessions of some species clustered together as expected, others were placed in separate clades or grouped more closely with different species, suggesting complex evolutionary relationships. Notably, accessions of species such as A. pendulum, A. flavum and A. kusnezoffii exhibited such inconsistent clustering patterns. Several factors may contribute to these patterns. First, misidentification or labelling errors in public databases can lead to incorrect species assignments. Second, incomplete lineage sorting and chloroplast capture due to hybridization55,56 can obscure phylogenetic signals, particularly when relying solely on chloroplast genomes.

The clustering of A. pendulum and A. flavum accessions, despite their designation as distinct species, aligns with recent findings from population genetics and ecological niche modelling, which suggest these taxa may represent a single species complex with weak genetic differentiation, historical gene flow, and evidence of a demographic bottleneck during the Last Glacial Maximum57. Similarly, the grouping of accessions from A. kusnezoffii (NC_031422), A. jaluense subsp. jaluense (KT820668) and A. japonicum subsp. napiforme (KT820670) supports earlier morphological and ecological studies from Mt. Sobaek in Korea, which indicate extensive hybridization and repeated introgression among these taxa58. Interestingly, the accessions mentioned above that form anomalous cluster in the core genes phylogeny, mainly belong to the Republic of Korea (Supplementary Fig. 6, Supplementary Table 11). These observations underscore the need for integrative taxonomic approaches in Aconitum, combining molecular data from both nuclear and organellar genomes with detailed morphological and ecological analyses to resolve species boundaries and evolutionary histories more accurately.

Our phylogenetic results also provide molecular support for the earlier taxonomic inference proposed by Gao et al. (2012), who, based on karyotypic and morphological evidence, suggested that A. angustius should be transferred from series Lycoctonia to series Volubilia within subgenus Lycoctonum59. They observed that A. angustius differs from A. sinomontanum var. sinomontanum in its often decumbent to twining stem habit, narrower and recurved upper sepal, and generally paler (whitish to pale purple) flowers—traits characteristic of series Volubilia. The consistent clustering of A. angustius with Volubilia members such as A. finetianum, A. longecassidatum, and A. quelpaertense in our chloroplast phylogeny provides independent molecular validation of Gao et al.’s conclusion, strengthening the case for its formal reassignment to Volubilia. Furthermore, A. angustius appears morphologically intermediate between A. sinomontanum and A. puchonroenicum. While A. angustius and A. sinomontanum share similar leaf and floral architecture, the twining habit and narrow, recurved upper sepal of A. angustius align it more closely with Volubilia. In contrast, A. puchonroenicum—though non-twining—also exhibits a mostly glabrous stem and leaves, racemose inflorescence, elongated and recurved upper sepal, glabrous petals with curved spurs, and overall floral structure consistent with the Volubilia group rather than the typically erect, purple-flowered Lycoctonia series members60,61. Collectively, both morphological and phylogenetic evidence suggest that A. angustius and A. puchonroenicum share derived characters characteristic of series Volubilia, supporting their reassignment from series Lycoctonia.

Our phylogenetic observations are largely consistent with those reported by Hong and Yang (2017), who, based on ITS sequence analyses, demonstrated non-monophyly within several series of subgenus Aconitum, including Stylosa, Volubilia, Inflata, and Ambigua62. The present study strengthens these findings using whole chloroplast genome data, providing robust molecular support for the complex and unresolved relationships among these series. Additionally, our results reveal non-monophyly of series Brachypoda, as A. coreanum does not cluster with other members of the series but instead forms an independent clade. This pattern aligns with earlier chemotaxonomic evidence from Xiao et al. (2006), which proposed that A. coreanum should be segregated from series Brachypoda owing to its distinct secondary metabolite profile—comprising predominantly C20-diterpenoid alkaloids, in contrast to the C19-diterpenoid alkaloids characteristic of other Brachypoda species63.

Overall, these observations highlight that while chloroplast phylogeny robustly supports subgeneric boundaries, finer-scale relationships at the series and species levels are confounded by hybridization, incomplete lineage sorting, or potential sample misidentification. The congruence between molecular and morphological data for A. angustius and A. puchonroenicum underscores the need for continued taxonomic re-evaluation of Aconitum series boundaries using integrative approaches combining plastome, nuclear, and morphological datasets. In conclusion, this study emphasizes the highly conserved yet subtly variable nature of chloroplast genomes within the Aconitum genus. The identification of accessory genes and annotation discrepancies underscores the limitations of current automated tools and the need for manual curation and better gene-calling and annotation tools. Codon bias patterns and Ka/Ks analyses reveal the interplay between mutational pressure and selection, particularly in genes like psbA. Phylogenetic incongruities suggest complex evolutionary histories shaped by hybridization and misidentification, calling for integrative taxonomic approaches.