Thirty-five individuals of A. annua were cultivated in Huairou District, Beijing, from seeds collected from four countries to test and verify the polymorphic loci (Supplementary Table 1). Fresh leaves were snap-frozen with liquid nitrogen and stored at −80°C (Eppendorf, Hamburg, Germany). The total DNA of A. annua individuals were extracted by the modified cetyltrimethylammonium bromide (3 × CTAB) method (Allen et al., 2006), and the DNA quality and concentration were measured by electrophoresis in 1.0% agarose gel and the NanoDrop2000 ultra-micro ultraviolet spectrophotometer (Thermo Scientific, MIT, United States). Primers were designed by Primer Premier 5 (Supplementary Table 2). PCR amplification was performed on 25 μL reaction mixtures containing 2X Pro Taq Master Mix (dye plus) 12.5 μL, nuclease-free water 8.5 μL, 1.0 μL each of 10 μM forward and reverse primers, and genomic DNA 1.0 μL. The PCR reaction conditions were listed in Supplementary Table 1. The amplified products were detected by electrophoresis in 1.0% agarose gel. The synthesis of primers, and sequencing of amplification products were conducted by Sangon Biotech Guangzhou branch office.

WGS (whole genome sequencing) datasets of LQ-9 and HAN1 strains, and whole genome resequencing datasets of 36 individuals were obtained from previous study (Liao et al., 2022a) (Supplementary Table 3). The quality of raw sequencing data was evaluated by FastQC 0.11.5 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc) and low quality bases and reads were trimmed and removed by Skewer (Jiang et al., 2014). The LQ-9 cp genome was assembled and as the reference genome for later analysis. WGS reads of LQ-9 were mapped to A. annua cp genome (GenBank Accession Number: MF623173) and mapped reads were extracted as cp-like reads. The extracted reads were then assembled into contigs by ABySS 2.0.0 (https://github.com/bcgsc/abyss) (Jackman et al., 2017). Finally, the cp sequence contigs were ordered and concatenated based on the collinearity with reference cp genome sequences. The initial gene annotation was conducted with plann 1.1.2 (Huang and Cronk, 2015) and then validated by BLAST and manually correction. The transport RNA (tRNA) genes were identified with tRNAscan-SE software (Lowe and Chan, 2016). Circular gene maps of the A. annua cp genomes were generated using Chloroplot software (Zheng et al., 2020). The cp genomes assembled in this study have been deposited in the Global Pharmacopoeia Genome Database (Liao et al., 2022b) at http://www.gpgenome.com/species/92 under the “SuperBarcodes” section.

A total of 117 complete cp genome sequences from 38 Artemisia species, and one complete cp genome from Chrysanthemum were downloaded from the NCBI GenBank (Supplementary Table 4) for phylogenetic analysis. Among which, 115 cp genomes were reassembled since the direction of SSC regions were opposite to A. annua cp genome in this study. Multiple sequence alignments were performed with 156 cp genomes (Supplementary Tables 3, 4) using MAFFT 7.313 (Katoh and Standley, 2013). Maximum likelihood (ML) phylogenetic tree was performed by RAxML 8.2.11 (Stamatakis, 2014) with 1,000 bootstrap replicates under the GTR + G model. Neighbor-joining (NJ) tree was constructed using MEGA 7 (Kumar et al., 2016) under the Kimura two-parameter model with 1,000 bootstrap replicates.

Whole cp genomes of 38 A. annua individuals were aligned using MAFFT 7.313 (Katoh and Standley, 2013). For rearrangement analysis, all aligned sequences were constructed with Geneious Prime program (Kearse et al., 2012). The IRscope program was used to evaluate the expansion and contraction of inverted repeat region (IR) (Amiryousefi et al., 2018). Comparative analysis of cp genomes was performed among Artemisia species and visualized by R. Furthermore, nucleotide diversity (Pi) was estimated by sliding window analysis using DnaSP 6.0 with a 600 bp window length and 200 bp step size (Rozas et al., 2017). High throughput sequencing (HTS) data of 38 A. annua individuals were aligned to the LQ-9 reference cp genome with Bowtie2 2.4.4 (Langmead and Salzberg, 2012). The variants were called and filtered using BCFtools (http://samtools.github.io/bcftools/) with a minor allele count higher than 1. False positives caused by sequencing errors were excluded by visualizing the sequencing data of mutation sites by Integrative Genomics Viewer (IGV) (Robinson et al., 2011).

A total of 38 complete A. annua cp genomes were assembled with WGS data of 38 distinct individuals released by previous study (Supplementary Table 3) (Liao et al., 2022a). To explore the phylogenetic and evolutionary relationships of A. annua populations, a phylogenetic tree was constructed using ML and NJ methods based on 156 complete cp genomes of 38 Artemisia species, and C. morifolium was included as the outgroup (Supplementary Table 4). The topological congruence between the ML and NJ phylogenetic trees was observed, with the majority of nodes exhibiting robust support values (>99%) (Figure 1; Supplementary Figure 1). Thirty-eight species were divided into two main branches, the first branch only comprising A. keiskeana (Figure 1, clade 1), while the second branch was further divided into seven well-supported clades (Figure 1, clades 2–8). Seriphidium was considered to be an independent genus (Lin et al., 2011; Haghighi et al., 2014), while this was not supported by the phylogenetic relationship of cp genomes as Seriphidium species were close related with A. schrenkiana and A. scopiformis. The phylogenetic tree also showed that all cp genomes from the same species clustered together, except A. argyi and A. selengensis. All A. annua cp genomes clustered into a same clade and shared the most recent common ancestor with A. fukudo and A. nakaii, which was consistent with the previous results of Shen et al. (2017) and Jin et al. (2023). Intriguingly, all A. annua strains were divided into two subgroups, forming two distinct lineages. Geographically, the first lineage (subgroup 1) containing a mixture of southern and northern strains was defined as a broad lineage, and the second (subgroup 2) was defined as a northern lineage containing only northern strains.

Thirty-eight assembled A. annua cp genomes were further annotated and analyzed. All cp genomes possessed the typical quadripartite structure which consisted of a large single-copy (LSC) region, a small single-copy (SSC) region, and two copies of IR regions (Figure 2). All cp genomes had lengths ranging from 150, 953 bp (LQ-9) to 150, 974 bp (RS16) with GC contents ranging from 37.47% to 37.48%. The LSC regions had lengths ranging from 82, 701 bp (RS18) to 82, 785 bp (RS24), and the SSC regions ranged from 18, 267 bp (RS24 et al.) to 18, 349 bp (RS18), while the lengths of IR regions were 24, 956 bp in all A. annua strains (Table 1), similar to the previously reported Artemisia cp genomes (Shen et al., 2017; Kim et al., 2020; Jin et al., 2023). Gene annotation showed that all cp genomes contained 131 genes, including 87 protein-coding genes, 36 tRNA genes and eight ribosomal RNA (rRNA) genes (Table 1). Of these, 17 genes (seven protein-coding genes, six tRNA genes and four rRNA genes) were duplicated in the IR regions in all A. annua strains (Supplementary Table 5). In addition, 17 intron-containing genes were identified, of which 15 contained one intron and two contained two introns (Supplementary Table 5).

The IR border structure was conserved in A. annua that no IR contraction was observed. In each A. annua individual, the junctions between IRs and LSC and SSC were flanked by rps19 and ycf1, respectively (Supplementary Figure 2), in concordance with other Artemisia species (Kim et al., 2020; Jin et al., 2023). Comparison of genome structure showed that A. annua cp genomes were highly conserved, and no significant gene rearrangements were observed (Supplementary Figure 3). On the whole, the A. annua cp genomes were structural conserved.

In this study, the intraspecific polymorphisms of A. annua cp genomes were assessed in inter-individual levels. Complete cp genomes of 38 A. annua individuals were aligned using the LQ-9 as the reference genome. In total, 60 polymorphic sites were identified after excluding false positives caused by sequencing errors, including 19 singleton variable sites and 41 parsimony informative sites, among which, 48 SNPs were identified in LSC, 11 in SSC and one in IR regions, respectively (Figure 3A; Table 2). The SSC regions were found to have the highest polymorphisms with a SNP site density of 8.8/10 kb, followed by LSC regions (5.3/10 kb), and the SNP density of IRs regions was only 0.4/10 kb (Table 2), revealing that the IR regions were more conserved than the single-copy regions, which is consistent with the comparative analysis of the whole cp genome of Artemisia species (Figure 4). Furthermore, the Pi values of 38 A. annua cp genomes were calculated using DnaSP v6.0, and the Pi values ranged from 0 to 0.0013 with an average of 0.00007. Based on DNA polymorphisms, seven highly diverged regions (Pi > 0.0007) were identified, including rps16-trnQ-UUG, petN-psbM, psbM-trnD-GUC, rpl20-clpP, clpP, trnL-UAG-rpl32 and rpl32 (Figure 3B), similar results have been observed in other Artemisia species (Shen et al., 2017; Kim et al., 2020; Jin et al., 2023). Fifteen protein-coding genes were found to containing variations, of which the ycf4 gene had the highest Pi value (0.0007) and the ndhF gene had the lowest Pi value (0.000002) (Figure 3C; Supplementary Table 6).

In addition, the WGS data of 38 individuals were aligned to LQ-9 reference genome with Bowtie2 for polymorphic sites validation. All 60 polymorphic sites were supported by mapped reads (Figures 3A, D). Notably, most of the polymorphic sites were single nucleotide polymorphisms (SNPs) between A. annua individuals, but partial single nucleotide polymorphism (pSNP) which contained more than one alternative bases in one site defined by James et al. (2009) were also found within individuals (Figure 3D). The polymorphic loci were confirmed by Sanger sequencing using the total DNA of 35 A. annua individuals (Figure 3F; Supplementary Table 1). In the Sanger sequencing results, the pSNP site generally contained multiple nested-peak which were identified as degenerate bases, and the SNP site showed clean single-peak representing complete base mutations (Figure 3E).

A total of 14 haplotypes were found in 38 A. annua individuals, which can be further divided into two distinct lineages as shown above (Figure 1; Figure 3G). Germplasm resources collected from the same location clustered together excluding Gansu, Hubei, and Xizang populations, but strains from different populations were also observed to have the same haplotype. Besides, polymorphism analysis showed that the subgroup 2 (the northern lineage) exhibited more unique polymorphic sites than the subgroup 1 (the broad lineage) and they shared only two polymorphic sites (Table 3). Other Artemisia species shared fifteen and four polymorphic sites with the northern lineage and the broad lineage, respectively (Table 3; Supplementary Figure 4). We proposed that the broad lineage may be the stable and specific evolutionary lineage during the speciation of A. annua, and the northern lineage may be the post-forming lineage. The above results showed that A. annua cp genomes were highly polymorphic, and the cp variants dataset can be used as the reference dataset for the identification and evolutionary analysis of A. annua strains.