P. kingianum with different morphological characteristics and other verticillate leaf-type Polygonatum species were obtained from Yunnan Province, China ( Supplementary Table 1 ). Fresh and healthy leaves from these plants were collected and stored in silica gel for DNA extraction. The voucher specimens were identified by Professor Yunheng Ji and deposited in the herbarium of Kunming Institute of Botany (Chinese Academy of Sciences). In addition, 12 Polygonatum plastomes were downloaded from the NCBI database for subsequent analysis ( Supplementary Table 2 ).

Total DNA was extracted from approximately 20 mg silica-gel dried leaf tissues using the modified CTAB method (Porebski et al., 1997). Sequence libraries were constructed using approximately 5 µg of purified genomic DNA with a TruSeq DNA Sample Prep Kit (Illumina, Inc., San Diego, CA, United States) according to the manufacturer’s recommendation. The NGS QC tool kit tool was used to trim the raw reads to obtain clean data (Patel and Jain, 2012).

We used NOVOPlasty v4.1.0 to perform cp genome assembly with a k-mer of 39 (Nicolas et al., 2017). Gene matK sequence from P. kingianum (GenBank: JX185477) as the seed for iterative extension of contigs to recover the whole plastome of each accession. The positions of start and stop codons were adjusted manually with a BLAST search against the NCBI protein database. All tRNA genes were identified by the tRNAscan-SE 2.0 with default parameters (Schattner et al., 2005). Thereafter, the boundaries of IR (inverted-repeat) and SC (single-copy) regions were defined by Geneious v. 9.0.2 (Kearse et al., 2012). Organellar Genome DRAW (Lohse et al., 2007) was used to draw the cp genome maps (http://ogdraw.mpimp-golm.mpg.de/). Additionally, all annotated cp genome sequences were submitted to GenBank on the NCBI website (https://www.ncbi.nlm.nih.gov/).

To verify the expansion or contraction of cp gene regions within verticillate leaf type of the Polygonatum species, the alignment was visualized using online mVISTA (https://genome.lbl.gov/vista/index.shtml) (Mayor et al., 2000) in Shuffle LAGAN mode, with the annotated cp genome of P. kingianum (Genbank number: MW566455) as a reference.

The protein-coding genes and intergenic spacers of the 20 plastomes (16 Polygonatum species) were extracted in Geneious v9.0.2 (Kearse et al., 2012) and aligned using the MAFFT v7.221 software (Katoh and Standley, 2013). Then, nucleotide diversity values (Pi) were calculated to estimate the divergence level of shared genes and intergenic spacers using DnaSP v5 (Librado and Rozas, 2009). The IRscope online program (https://irscope.shinyapps.io/irapp/) was used to evaluate IR expansion and contraction of plastome genetic architecture (Amiryousefi et al., 2018).

Online program (https://webblast.ipk-gatersleben.de/misa/): MISA software was used to detect simple sequence repeats (SSRs) in the cp genomes (Beier et al., 2017). Parameter setting thresholds of 10, 5, 4, 3, 3, and 3 were set for mono-, di-, tri-, tetra-, penta-, and hexanucleotide SSRs, respectively.

The amount of codon and relative synonymous codon usage (RSCU) were calculated by CodonW v.1.4.2 with default parameters. Moreover, the heatmap of the RSCU values was shown by using TBtools software (Chen et al., 2020). Following the previous research methods, protein-coding genes (PCGs) of less than 300 nucleotides in length and the repetitive gene sequences were removed to reduce the deviation of the results (Fang et al., 2023).

All 21 available plastome sequences, including 12 published sequences downloaded from NCBI and 9 sequences obtained in this study, were included in the phylogenetic analysis. We constructed the phylogenetic trees by Maximum likelihood (ML) and Bayesian inference (BI) using the entire cp genomes. P. cyrtonema, the relative of P. kingianum, was used as the outgroup to root the phylogenetic trees. The ML tree was built using the CIPRES Science Gateway website (https://www.phylo.org/) under a General Time Reversible + Gamma (GTR + G) model and estimated with 1,000 bootstrap replicates (Stamatakis, 2006). The BI analysis was performed with Mrbayes 3.2.7 under the best substitution models and parameters (Posada and Crandall, 1998; Ronquist and Huelsenbeck, 2003; Posada and Buckley, 2004). The resulting phylogenetic trees were visualized by FigTree v1.4.4 software (http://tree.bio.ed.ac.uk/software/figtree/).

The size of newly sequenced plastid genomes within this study varied from 155,437 to 155,977 bp ( Supplementary Table 1 ). All of them exhibited typical quadripartite structures that were similar to the pattern observed in other typical cp genomes of angiosperms, consisting of a pair of inverted repeat regions from 26,278 bp to 26,394 bp, a long single copy region 84,406 bp to 84,742 bp, and a short single copy region 18,415 bp to 18,554 bp ( Supplementary Table 1 ; Figure 1 ). The GC contents of the nine plastomes were very similar and ranged from 37.6% to 37.7% ( Supplementary Table 1 ). In this study, 133 genes were annotated in these cp genomes, including 86–87 protein-coding genes, 38 genes encoding transfer RNA, and 8 genes encoding ribosomal RNA. All encoded genes were divided into four categories according to their functions, including self-replication genes, genes related to photosynthesis, other genes, and unknown function genes ( Table 1 ). Among these genes, 18 genes were intron-containing genes, 15 genes containing one intron (rpl2, rpl16, rpoC1, rps16, trnA-UGC, trnG-UCC, trnI-GAU, trnK-UUU, trnL-UAA, trnV-UAC, atpF, ndhA, ndhB, petB, and petD), and others containing two introns (rps12, clpP, and ycf3) ( Table 1 ). All infA genes were found as pseudogenes in verticillate leaf types of the Polygonatum species except the P. sibiricum. This is because a G base was inserted into the infA sequence, which caused the infA gene to contain in-frame stop codons. All sample cp genome sequences were deposited in the GenBank database of the NCBI ( Supplementary Table 1 ).

The expansions and contractions of IR borders were considered to be evolutionary events and the main cause of cp genome length changes (Wang et al., 2008). We analyzed the border regions of 16 verticillate leaf types of Polygonatum species plastomes to find similarities and differences ( Figure 2 ). In our study, the rps19 and rpl22 genes were located in the IRb/LSC borders. In the SSC/IRa border, rps19 and psbA genes were found. The ndhF crossed the IRa/SSC region. The ycf1 gene was a cp-genome-encoded giant open reading frame that crossed the IRb/LSC region. Both the pseudogene fragment ψycf1 and the ndhF gene were located at the IRb/SSC border and partially overlapped in the verticillate leaf type of the Polygonatum species. In P. sibiricum, the two rps19 genes extended into the LSC region by 27 and 60 bp, respectively. However, the rps19 genes of other species were entirely located in the IR regions.

In order to further analyze the results, the 16 verticillate leaf types of the Polygonatum species (20 plastid sequences) were aligned and visualized to determine the overall variations by the program mVISTA. In agreement with other studies (Huang et al., 2021; Wu et al., 2021; Yang et al., 2021), the result of the alignment revealed that the coding regions are more highly conserved than non-coding regions. Most of the variations detected were found in non-coding sequence areas ( Figure 3 ).

The nucleotide polymorphism (Pi) values of the 20 cp genomes in this study were computed to detect hyper-variable regions using the DnaSP v5 software (Librado and Rozas, 2009). Figure 4 shows that the IR regions were more conserved than the LSC and SSC regions. The average Pi value of shared PCGs ranged from 0 to 0.01113, and the average Pi value of intergenic spacers ranged from 0 to 0.03551 ( Supplementary Table 3 ; Figure 4 ). These hotspots (PCGs, Pi > 0.0036; Intergenic spacers, Pi > 0.02) were selected including ccsA, psbF, ycf1, psbM, rps19, rpl16, rrn4.5–rrn5, trnG-GCC–trnfM-CAU, ccsA–ndhD, atpA–atpF, and rps19–psbA, which were potentially useful for species identification in P. kingianum with related verticillate leaf type of the Polygonatum species.

SSRs are useful for detecting genetic diversity and polymorphisms at the population, intraspecific, and cultivar levels, as well as for species-distinguishing (Zaloglu et al., 2015; Guo et al., 2018; Guang et al., 2019). The total number of simple sequence repeats were identified in cp genomes by Perlscript MicroSAtellite (MISA) (Beier et al., 2017) and ranged from 64 to 81, with the most abundant type of SSRs being the mono-nucleotide repeats ( Figure 5B ; Supplementary Table 4 ). In each case, the 37–53 mononucleotide SSRs, 13–17 dinucleotide SSRs, 2–5 trinucleotide SSRs, 7–10 tetranucleotide SSRs, and 2 pentanucleotides SSRs were identified. There were no hexanucleotide repeats detected in the verticillate leaf type of the Polygonatum species. Most repeats of all six microsatellite types were of the A/T motif rather than the G/C motif ( Figure 5A ; Supplementary Table 4 ). SSRs composed of G/C were not found in the plastomes of P. griffithii, P. tessellatum, P. hookeri, P. zanlanscianense, P. cathcartii, P. sibiricum, P. curvistylum, P. kingianum (OQ532971), P. kingianum (MW566455), and P. kingianum (OQ532967), while trinucleotide SSR AGC/CTG was only present in P. curvistylum, P. stewartianum, P. hookeri, and P. cirrhifolium. The tetranucleotide SSR AAAG/CTTT was only detected in the cp genome of the P. stenophyllum.

A total of 56 protein-coding genes were selected for codon usage analysis. A total of 64 codons were observed, of which 61 codons encoded 20 amino acids and 3 stop codons ( Figure 6 ; Supplementary Table 5 ). These sequences were 72,513–72,885 bp in length and encoded 24,171–24,295 codons. Of these amino acids, leucine had the highest number of codons, whereas cysteine had the least ( Supplementary Table 5 ). Figure 6 presents the relative synonymous codon usage (RSCU) of each plastome. RSCU values differ slightly in these plastid sequences. Additionally, the results showed that the RSCU values of 31 codons were higher than 1 and that 29 codons presented A/U- endings. Tryptophan and methionine were the amino acids without codon bias (RSCU = 1).

We employed the whole cp genome sequences to build the phylogenetic tree to determine the phylogenetic relationship. Maximum likelihood (ML) and Bayesian inference (BI) were used to perform phylogenetic analyses with P. cyrtonema as an outgroup. The topologies of the ML and BI trees were identical ( Figure 7 ). P. kingianum with different morphological characteristics clustered together with high support. P. kingianum was closely related to P. stenophyllum and P. sibiricum (BS = 100, PP = 1.00). P. griffithii and P. cathcartii clustered into a subclade and formed a sister relationship with the subclade of P. oppositifolium and P. tessellatum. P. punctatum and P. annamense clustered into a subclade; P. verticillatum and P. zanlanscianense clustered into a subclade; and P. hookeri, P. prattii, P. stewartianum, P. curvistylum, and P. cirrhifolium clustered into a subclade.

In this study, the cp genomes of nine samples were sequenced by next-generation sequencing (NGS) methods. Similar to previous studies (Huang X. et al., 2022; Zhang et al., 2022; Yang et al., 2023), the newly sequenced cp genomes of the samples were characterized, which displayed the typical quadripartite structure of flowering plants, and cp genomes contained more AT content than GC content. The nine cp genomes were similar both in size and composition. In the current study, we found that most of the infA genes in these cp genomes were pseudogenes except P. sibiricum. Pseudogenization is common in cp genomes ( Abdullah et al., 2021;Li X. et al., 2021; Scobeyeva et al., 2021). In the present study, base insertion was responsible for the pseudogenization of the infA genes. All in all, the cp genomes were highly conservative similar to the majority of plants.

These cp genomes underwent a comparative analysis with sequences from Genebank. We compared the boundary between the SC and two IR regions and their neighboring genes. There was no significant difference at the boundary of most species except in the LSC and IR junctional regions of P. sibiricum. Previous research demonstrated that many rps19 genes of the genus Polygonum were located entirely in the IR region (Wang et al., 2022; Xia et al., 2022). In contrast, the two rps19 genes of P. sibiricum were located partly in the LSC Region, possibly due to IR contraction. Sequence divergence among plastid genomes comparisons based on mVISTA and nucleotide diversity showed that the low-identity region was mainly distributed in the non-coding region of the SC region compared with the IR region through multiple sequence alignment. Mutational hotspots of the cp genome are often developed to distinguish closely related species or genera, which can serve as specific DNA barcodes (Dong et al., 2012). Our study identified eleven regions with a higher evolution rate, namely, ccsA, psbF, ycf1, psbM, rps19, rpl16, rrn4.5–rrn5, trnG-GCC–trnfM-CAU, ccsA–ndhD, atpA–atpF, and rps19–psbA, all of which might be potentially used as candidate markers for identifying P. kingianum and its related verticillate leaf type of the Polygonatum species.

Cp microsatellites (cpSSRs) exist widely in cp genomes and have the advantages of codominant inheritance, high polymorphism, and abundance (Huang et al., 2018). Therefore, cpSSR markers were widely used for species authentication, evolutionary studies, population genetics, and so on (Dzialuk et al., 2009; Jayaswall et al., 2015; Singha et al., 2017; Mahajan et al., 2022).

In our study, we found that cpSSRs abundances in different species were varied. Even within the same species (P. kingianum), the number and type of cpSSRs were different. A or T mononucleotide repetition was the most common primary repetitive type. C or G mononucleotide repetition existed in some species and was rare in number. Previous studies have suggested that this might be due to the convenience of A-T conversion over G-C in the cp genome (Huang Y. et al., 2022). The number of cpSSRs of five kinds of SSRs were discovered: mononucleotide > dinucleotides > tetranucleotide > trinucleotide > pentanucleotide. The cpSSRs detected in this study would be helpful for the identification of the verticillate leaf type of the Polygonatum species, for germplasm resource evaluation, and for the conservation of P. kingianum.

Codon usage biases are found in all eukaryotic and prokaryotic genomes and refer to the unequal use phenomenon of synonymous codons in an organism (Zhou et al., 2016). Codon usage biases are associated with the gene expression level and accuracy of the translation, gene length, gene translation initiation signals, protein amino acid composition, protein structure, tRNA abundance, mutation frequency and patterns, and GC composition (Romero et al., 2000; Xu et al., 2011). These results suggested that all individuals share similar codon usage patterns, with only slight differences in the number of codons and the RSCU value. The AT-rich bias was strongest in the third codon position, and the results of our study are similar to those of other species with cp genome codon usage biases (Huang et al., 2020; Dong et al., 2021; Tang et al., 2021). The research on codon preferences could help us provide a theoretical basis for codon modification of exogenous genes, accuracy of prediction about new members of the cp gene family, and identification of unknown genes.

Phylogenetic analysis is of great significance for clarifying genetic relationships and species identification (Lan et al., 2022). Despite their morphological diversity, P. kingianum could be clustered into one clade by the phylogenetic tree based on cp genomes. We inferred that the limited geographical distribution, short-distance pollination by insects, and long-distance seed dispersal by birds might contribute to the frequent intraspecies gene flow of P. kingianum. These factors could decrease genetic differences among different morphological P. kingianum, improving the diversity and adaptability of the species. The cp genome of P. griffithii was reported for the first time in this study. The plastome phylogeny found that P. griffithii and P. cathcartii were closely related, which confirmed the speculation of FOC (Flora of China). All verticillate leaf types of the Polygonatum species in this research paper, including P. kingianum, have different degrees of overlap in geographical distribution, and many of them also have morphological intermediates types (e.g., P. sibiricum, P. verticillatum, P. cirrhifolium, and P. zanlanscianense). These species were morphologically indistinguishable. Our high-resolution phylogenetic trees not only distinguished these species but also provided some insight into the verticillate leaf type of the Polygonatum species relationship.