The Complete Chloroplast Genomes of Primula obconica Provide Insight That Neither Species nor Natural Section Represent Monophyletic Taxa in Primula (Primulaceae)

The genus Primula (Primulaceae) comprises more than 500 species, with 300 species distributed in China. The contradictory results between systematic analyses and morphology-based taxonomy make taxonomy studies difficult. Furthermore, frequent introgression between closely related species of Primula can result in non-monophyletic species. In this study, the complete chloroplast genome of sixteen Primula obconica subsp. obconica individuals were assembled and compared with 84 accessions of 74 species from 21 sections of the 24 sections of the genus in China. The plastome sizes of P. obconica subsp. obconica range from 153,584 bp to 154,028 bp. Genome-wide variations were detected, and 1915 high-quality SNPs and 346 InDels were found. Most SNPs were detected in downstream and upstream gene regions (45.549% and 41.91%). Two cultivated accessions, ZP1 and ZP2, were abundant with SSRs. Moreover, 12 SSRs shared by 9 accessions showed variations that may be used as molecular markers for population genetic studies. The phylogenetic tree showed that P. obconica subsp. obconica cluster into two independent clades. Two subspecies have highly recognizable morphological characteristics, isolated geographical distribution areas, and distinct phylogenetic relationships compared with P. obconica subsp. obconica. We elevate the two subspecies of P. obconica to separate species. Our phylogenetic tree is largely inconsistent with morphology-based taxonomy. Twenty-one sections of Primula were mainly divided into three clades. The monophyly of Sect. Auganthus, Sect. Minutissimae, Sect. Sikkimensis, Sect. Petiolares, and Sect. Ranunculoides are well supported in the phylogenetic tree. The Sect. Obconicolisteri, Sect. Monocarpicae, Sect. Carolinella, Sect. Cortusoides, Sect. Aleuritia, Sect. Denticulata, Sect. Proliferae Pax, and Sect. Crystallophlomis are not a monophyletic group. The possible explanations for non-monophyly may be hybridization, polyploidization, recent introgression, incorrect taxonomy, or chloroplast capture. Multiple genomic data and population genetic studies are therefore needed to reveal the evolutionary history of Primula. Our results provided valuable information for intraspecific variation and phylogenetic relationships within Primula.


Introduction
Species are fundamental units of biodiversity, but the boundaries of species are still confusing [1], especially in species-rich lineages that share significant similarities, morphological characteristics, and ecological regions. Generally, prior species delimitation is based on morphological traits, which may fail to discriminate polymorphic taxa [2]. Due to the rapid development of next-generation sequencing technologies, chloroplast genomes have become easier to obtain. As "super-barcodes", chloroplast genomes have been successful in solving taxonomic issues of closely related plant taxonomy [3], such as in Ottelia acuminata, which consists of six phenotypic varieties, of which the chloroplast genome was successfully divided it into two conspecific varieties [4], thus providing new insight into the use of the chloroplast genome as a super-barcode for species delimitation in taxonomically difficult plant taxa. the intraspecific variation in P. obconica subsp. obconica; (2) decipher the taxonomic delimitation of the P. obconica complex; (3) reveal the phylogenetic relationships of 21 sections of Primula in China and test morphology-based taxonomy with molecular data.

Sample Collection, DNA Extraction, and Sequencing
Thirteen wild accessions were collected from the Sichuan, Yunnan, Hubei, and Hunan provinces. Samples were collected across China, including the main distribution areas of the species. Sichuan, Hubei, and Hunan included two populations, and each population contained two wild accessions. Two Chinese-cultivated accessions were collected from the Wuhu (Anhui) flower market (Table S1). All specimens were deposited in the Qinghai-Tibetan Plateau Museum of Biology, Northwest Institute of Plateau Biology, Chinese Academy of Sciences (HNWP). The silica gel-dried leaves were sent to Genepioneer Biotechnologies (Nanjing, China). DNA was extracted using a plant genomic DNA Kit (DP305), and paired-end reads of 2 × 150 bp for all samples were sequenced using NovaSeq 6000 (Illumina Inc., San Diego, CA, USA). The clean data were from 1.38 to 2.93 GB.

Genome Assembly and Annotation
The de novo assembly of the clean data was performed on a GetOrganelle pipeline [29], with -t 1, -R 10, using the P. obconica subsp. obconica (MK344754) plastome as a reference. Bandage was used to confirm that the assembly graph was circular [30]. Assembled genomes were annotated using CPGAVAS2 [31]. OrganellarGenomeDRAW version 1.3.1 was used to visualize the plastomes of P. obconica [32].

Comparative Analysis of the Chloroplast Genomes
A comparative analysis of the chloroplast genomes of sixteen wild and cultivated accessions was carried out using mVISTA online software (https://genome.lbl.gov/vista/ mvista/submit.shtml, accessed on 30 August 2021) with a shuffle-LAGAN model, and HBD1 as a reference. To compare the IR/SC boundaries of different subspecies within P. obconica, the IRscope web service (https://irscope.shinyapps.io/irapp, accessed on 7 September 2021) was used to visualize the result.

Population Structure Analysis of P. obconica subsp. obconica
From the SNPs and InDels used in the aforementioned step, 1361 SNPs and 321 InDels were analyzed after filtering with the parameters mac 3 and max-missing 0.5. Sixteen sequenced accessions were converted into a concatenated sequence in Fasta format using PGDSpider version 2.1.1.5 [38]. The ML tree based on SNPs was constructed by IQ-TREE version 2 [39], with the K3P model. PCA was conducted by Plink version 1.9 [40]. Gene flow was inferred using TreeMix version 1.12 [41].

Phylogenetic Analysis
To study the relationships between Primula, we used 100 accessions of 74 species from 21 sections of the 24 sections of the genus Primula in China. An ML phylogenetic tree was performed, using Pentaphylax euryoides (MF179498) as an outgroup. Seventy-seven protein-coding genes were extracted, aligned, and merged by PhyloSuit version 1.1.2 [42]. Maximum likelihood (ML) estimation with 1000 bootstrap replications was conducted using IQ-TREE version 2 [39], with a TVM + F + R3 model from the result of ModelFinder.

General Features of the P. obconica subsp. obconica Chloroplast Genome
The 16 accessions of P. obconica subsp. obconica had a conventional quadripartite structure characteristic of angiosperm plastomes, with an LSC, SSC, and two IR regions ( Figure 1). The chloroplast genomes ranged from 153,584 bp (SCB1) to 154,028 bp (HNL2 and HND1). The LSC and SSC regions have a relatively higher length variation compared with the IR regions, for example, the length of the LSC region was 84,712 bp and 84,403 bp in SCD1 and YN1 chloroplast genomes, respectively. The SSC region ranged from 17,743 bp (HNL2, HND1) to 17,347 bp (SCD1) ( Table S2). A total of 128 genes were detected, including 83 protein-coding genes, eight ribosomal RNAs, and 37 transfer RNAs (Table S3). There were six protein-coding genes with the duplication genes and 22 genes that contained one intron. Two of these genes (ycf3 and clpP) contained two introns. The ycf3 and clpP genes were also trans-spliced into three exons (Table S4).

Chloroplast Genomic SNP Variations
The clean data of 16 accessions were mapped onto the Primula obconica subsp. obconica (MK344754) chloroplast genome reference. A total of 1915 high-quality SNPs and 346 InDels were detected (Table S5). ZP2 harbored the largest number of genomic variants (1002), while the YN population (YN1 and YN2) possessed the fewest number of SNPs (335). The LSC region had a relatively higher number of SNPs compared with other regions. In the protein-coding gene region, there were 15 protein-coding genes which had an abundance of SNPs (matK, atpF, rpoC2, rpoC1, rpoB, psaA, ycf3, clpP, rpl16, ndhF, ccsA, ndhD, ndhA, ndhH, and ycf1). Moreover, ycf1 contained the most abundant number of SNPs (133) (Figure 2A). In the intergenetic region, matK-rps16, rps16-psbK, psbI-atpA, psbM-psbD, rps4-ndhJ, ndhC-atpE, rbcL-psaI, psbE-petL, and ndhF-rpl32 had more sequence variations, and the psbM-psbD intergenic spacers had the most abundant number of SNPs ( Figure 2B). The results of SnpEff showed the total count of effects at 21371, and the plastome had an SNP density of 79 per bp. The ratio of non-synonymous to synonymous mutations was 1.1667. Most SNPs were detected in downstream and upstream gene regions (45.549% and 41.910%) ( Table S6). One hundred and thirty SNPs were shared among all accessions. The Yunnan and Sichuan populations (YN, SCD, SCB) shared more variations in their genomic sequences than other accessions. The HN population possessed no species-specific SNPs, and the YN population possessed the most species-specific SNPs ( Figure 2C). Most of the shared SNPs were detected in the ycf1 gene (9), but when compared with the overall SNPs detected in the ycf1 gene (133), it seems that the ycf1 gene had more variations at the intraspecific level ( Figure 2D).

Simple Sequence Repeat Analysis
A total of 401 SSRs were identified, including 317 mononucleotide, 31 dinucleotide, 39 tetranucleotide, and 14 pentanucleotide repeats ( Figure S1A). Two cultivated species (ZP1 and ZP2) possessed the highest number of SSRs (51 and 49). Moreover, ZP1 and ZP2 possessed the greatest number of mononucleotide repeats when compared with other accessions. Mononucleotide repeat units (A/T) were the most abundant. Only one mononucleotide repeat unit (C), and one pentanucleotide repeat unit (TATAA) were highlighted in the cultivated accession of ZP1. Mononucleotide repeats were detected in ZP1, ZP2, HND1, and HNL1, respectively. Moreover, tetranucleotide repeat units (AAAT) only appeared in HBE1. However, no SSR motifs for tri-and hexanucleotide repeats were identified ( Figure S1B). We found 15 SSRs that were shared in the nine accessions, twelve of which showed length variations between accessions (Table S7).

IR Expansion and Contraction
Studies on the expansion and contraction variation in junction regions between separate species are common but rarely focused on the intraspecific level. Although differently geographically distributed accessions of Primula subsp. obconica exhibited the same pattern in content and the number of genes in the IR/LSC and IR/SSC boundary regions, some differences still existed (Figure 3), such as the ndhF gene spanning the IRb/SSC border and the length ranging from 2177 to 2192 bp. The ycf1 gene spanned the SSC/IRa border, from 4533 bp to 4573 bp in the SSC region, and from 1004 bp to 1017 bp in the IRa region. The rps19 gene was located in the LSC region and extended 5 bp in the IRb region, except for P. obconica subsp. parva and P. obconica subsp. begoniiformis with 3 bp extended.

Population Structure Analysis
In total, 1361 SNPs and 321 InDels were used to investigate genomic evolution among the 16 Primula obconica subsp. obconica after filtration. We constructed an ML phylogenetic tree based on SNP data ( Figure 5A). Two major branches are shown in the phylogenetic trees. The Sichuan (SCB and SCD) and Yunnan populations (YN) were gathered in one clade, implying a close relationship between the two populations. The Hubei populations (HBE and HBD), the cultivated accessions (ZP1 and ZP2), and the Hunan populations (HND and HNL) showed closely-related relationships. A PCA analysis based on the SNPs and InDels was also performed to investigate the population genetic structure in 16 P. obconica subsp. obconica accessions ( Figure 5B,C). The 16 P. obconica subsp. obconica accessions could be divided into two groups, which were consistent with topology A of the phylogenetic tree. To understand the history of the divergence, TreeMix was applied to investigate the migration events within the P. obconica subsp. obconica accessions, with the SCD population as the outgroup taxa ( Figure 5D). The TreeMix result showed a relatively high gene flow that appeared in the HBE and Hunan populations (HNL and HND). A weak gene flow from SCB to the Hunan populations and YN to HBE was observed. Furthermore, we did not observe gene flow between the cultivated accessions, and the SCD and HBD populations.

Phylogenetic Analyses
An ML tree was constructed using 77 protein-coding genes that were shared by 100 accessions (Figure 6). The vast majority of the nodes in the phylogenetic trees received more than 90% bootstrap support values. Three major branches are shown in the phylogenetic trees, both with identical tree topologies.
Minutissimae forming a monophyletic group in cluster II.
The rest of the sections (Sect. Proliferae Pax, Sect. Petiolares, Sect. Amethyatina, Sect. Crystallophlomis, Sect. Ranunculoides) formed cluster III, in which Sect. Proliferae Pax and Sect. Crystallophlomis did not form a monophyletic section. Four species (P. chrysochlora, P. helodoxa, P. miyabeana, P. wilsonii) in Sect. Proliferae Pax formed a clade sister to Sect. Amethyatina. One species is of special interest in cluster III: P. handeliana is ascribed to Sect. Crystallophlomis, but it formed a distinct clade with Sect. Crystallophlomis and was more closely related to Sect. Ranunculoides. The monophyly of Sect. Ranunculoides is strongly supported in this study.

Discussion
Plastome variations within species have rarely been reported. The length of plastomes, the gene number, and the content within species were conserved. The plastomes of P. obconica all contain 128 genes, including 83 protein-coding genes, eight ribosomal RNAs, and 37 transfer RNAs. Prior research has indicated that the length variation in the whole plastome is mainly affected by the difference of LSC region [7], but our results showed that both the LSC and SSC regions play an important role in plastome size variation within species. The LSC region of HBE1 was 309 bp less in size than HBL2 and HBD1. The SSC region among HBL2, HBD1, and SCD1 varied by 396 bp in length, and the IR regions were more conserved, with only 48-bp length variations at the intraspecific level. The length varied between different accessions of P. obconica subsp. obconica, which may be affected by geographical isolation.
Although plastomes have the properties of uniparental inheritance and a low frequency of recombination, these features have also prompted plastomes to maintain highly conserved features; however, mutations at the intraspecific level have still occurred, such as in Eragrostis tef (12 InDels and 19 SNPs) [43], Selaginella tamariscina (1641 InDels and 1213 SNPs) [44], and five Goodyera schlechtendaliana (414-2133 InDels and 200-844 SNPs) [45]. A total of 1915 high-quality SNPs and 346 InDels were detected in 16 accessions of P. obconica subsp. obconica. The cultivated accession ZP2 harbored the largest number of SNPs (1002), and wild accessions of YN (YN1 and YN2) had the fewest number of SNPs (335). In contrast, the genomic variants detected in wild rubber trees (193) were much larger than those detected in cultivated (91) accessions [46]. Many of the morphological characteristics shown in the cultivated accessions were different from the wild populations; for example, the leaves and flowers were larger than those of the wild accessions. Moreover, the inferred ML tree displayed more distant relationships between ZP and P. obconica subsp. obconica (MK344754) than other accessions. Among the protein-coding genes, 15 genes had a relatively higher number of SNPs (matK, atpF, rpoC2, rpoC1, rpoB, psaA, ycf3, clpP, rpl16, ndhF, ccsA, ndhD, ndhA, ndhH, and ycf1). Moreover, SNPs were especially abundant in the ycf1 gene (133), which was also reported in Ricinus communis [7]. The matK, psbA, ndhF, and ycf1 genes were used as cpDNA barcodes to identify closely related species [47][48][49][50]. Moreover, the ycf1 gene was considered as the core barcode of land plants [50]. These cpDNA barcodes not only greatly varied between the interspecies of the genus Primula but were also varied among the different varieties of P. obconica [51]. In total, 1915 SNPs were detected in 16 plastome genomes of P. obconica subsp. obconica, but only 131 SNPs were common, which showed higher variation at the intraspecific level of P. obconica subsp. obconica. These samples were collected from different geographic locations, and climatic variation may be the main driver triggering genomic variations [52]. The SnpEff annotation results showed that most of the SNPs were distributed in the downstream and upstream regions (45.549% and 41.910%), which was consistent with the Triticum-Aegilops complex plastome genomes [53].
Due to their high abundance and variability, especially of the repeat and motif structures of SSRs in the cp genome, the SSRs have been used for population genetics and evolution studies [54,55]. Seven individuals from different geographic locations and two cultivated accessions were analyzed. A total of 401 SSRs were detected, and two cultivated accessions had the highest abundance of SSRs, which was consistent with cultivated rubber [46]. Our results are consistent with previous research, which described that mononucleotide repeats representing the most common repeat type of SSRs and chloroplasts generally consisted of short polyA/T repeats [56,57]. Furthermore, 15 SSRs shared in the nine accessions were detected, and 12 of these had length variations between those accessions. These variations may be helpful for population genetic analysis in P. obconica subsp. obconica.
The expansion and contraction between IR and the single-copy (SC) boundary regions are considered the main causes of the cp genome variation [60]. There has been considerable research indicating that genes in the IR and single-copy boundary regions vary in different taxa [61,62], such as in cultivated and wild Hevea brasiliensis, in which there are variations in gene number and the length of ycf1 spanning the boundary of the SSC region [46]. In our study, the junctions between the IR and single-copy boundary regions were highly conserved between the wild and cultivated P. obconica subsp. obconica. There were apparent variations in some genes spanning the boundary of the SC/IR region. These variations are likely associated with phylogenetic signals, for example, the ndhF spanning the boundary of the IRb/SSC region was consistent with the clade of the phylogenetic tree in P. obconica subsp. obconica.
Southwest China is considered the biodiversity center of Primula and the place of origin of other land plants [63,64]. Migration events between Southwest China and East China may occur along several mountain ranges, such as the Qin-Daba Mountains and Dalou Mountains [65]. Considerable research efforts have demonstrated that P. obconica originated from the Yunnan and Sichuan provinces [66,67]. Two major branches of P. obconica were shown in phylogenetic trees and PCA analysis, which were based on data regarding SNPs, InDels, and 77 shared protein-coding genes. Our results were consistent with the finding that P. obconica was split into the Sichuan-Central group and Yunnan-Eastern group [67]. There are a series of high mountains that may be a geographical barrier between Yunnan and Sichuan, which may have triggered intraspecies genetic divergence between the Yunnan lineages and Sichuan lineages [68]. Our phylogenetic tree showed that P. obconica subsp. obconica was separated into two groups, consistent with the result of our SSR analysis [69]. Non-monophyletic species make it difficult to infer phylogenetic relationships [70]. Whether monophyly can be used to define species is also controversial and non-monophyletic species have also been considered as "wrong" taxonomy [71]. However, non-monophyletic species occur in many taxa of the family Primulaceae, such as Dodecatheon dentatum, Dodecatheon jeffreyi, Primula veris, and Primula vulgaris [72,73]. The explanations for non-monophyletic species are hybridization, ancestral polymorphism, or both processes [72]. Moreover, geographical isolation and environmental heterogeneity each play a crucial role in species diversification [74]. Based on the insufficiency of our data, it is not possible to reliably infer the reasons behind the non-monophyletic occurrence of P. obconica subsp. obconica.
Many factors can affect genetic diversity within species, such as geographic distribution, breeding system, and isolation by gene flow [75]. Gene flow plays an important role in the speciation and adaptation of land plants [76]. The extent of gene flow between species could be confused with the interspecific boundaries and can bring difficulties to species delimitation [77]. In contrast, frequent gene flow has occurred within species, which could reduce intraspecific divergence, trends which are helpful for species delimitation [78]. High levels of genetic variation are common in the genus Primula, such as P. sikkimensis (Shannon's index H SP : 0.5576, expected heterozygosity Hj: 0.4032) [75], and P. interjacens (H SP : 0.4618) [79]. According to previous research, there is a lower genetic diversity within populations than between populations of P. obconica. Prior research has implied that the genetic differentiation in P. obconica was mainly caused by geographical barriers and seed dispersal mechanisms (dispersed by gravity, ants, or rodents) [67,80]. In this study, a relatively high gene flow had appeared from the Hunan populations (HNL and HND) to HBE. The two populations had closely related relationships and geographical distances. Low-frequency gene flow was detected in clades A and B, which confirmed the results of the geographical barriers and seed dispersal mechanisms, exacerbating the process of divergence in P. obconica subsp. obconica. P. obconica subsp. werringtonensis is distributed in the west of Sichuan, with a margin of usually sinuate-lobulate that is distinguished from P. obconica subsp. obconica. P. obconica subsp. begoniiformis possesses slender petioles, a leaf blade of ovate-rotund to suborbicular shape, and a margin of crenate-lobulate that differs from P. obconica subsp. obconica. P. obconica subsp. parva and P. obconica subsp. fujianensis all have the characteristic of scapes that are shorter than leaves. The four subspecies all have very small populations and are restricted to the Yunnan, Sichuan, and Fujian provinces of the Chinese mainland. A previous study inferred that P. obconica subsp. fujianensis originated from P. obconica subsp. obconica, and the east of Yunnan populations originated from the Sichuan population. Moreover, the variation in genetics between the four subspecies and P. obconica subsp. obconica was not obvious [25]. In our phylogenetic tree result, only P. obconica subsp. fujianensis and P. obconica subsp. werringtonensis were nested within P. obconica subsp. obconica, suggesting the two subspecies had genetic variation overlapping the P. obconica subsp. obconica populations. Additionally, P. obconica subsp. parva and P. sinolisteri subsp. sinolisteri were clustered in a sister clade, while P. obconica subsp. begoniiformis and P. ambita were closely related. The two subspecies have highly recognizable morphological characteristics, isolated geographical distribution areas, and a distinct phylogenetic relationship compared with P. obconica subsp. obconica. We propose the elevation of the two subspecies of P. obconica to separate species, P. begoniiformis and P. parva.
Our phylogenetic tree showed multiple non-monophyletic sections of Primula, such as Sect. Obconicolisteri, Sect. Monocarpicae, Sect. Cortusoides, Sect. Carolinella, Sect. Aleuritia, Sect. Denticulata, and Sect. Crystallophlomis. Three possible explanations may aid in the interpretation of these phenomena. Firstly, genetic drift, ancestral polymorphism, polyploidization, and recent introgression lead to sister species becoming reciprocally monophyletic. Gene exchange among interspecies may lead to the combination of genetic components of different ancestors, in which case, phylogenetic analysis between hybrids will be more challenging. Previous research has demonstrated that natural hybridization is common and has been detected between interspecies in Primula [21,81,82]. The main pollinators of Primula species were bees and butterflies, there were 12 families and 22 species of flower-visiting insects of Primula lithophila [83]. Multiple pollinators have increased the opportunity of Primula outcrossing. Polyploidization has been detected in Primula, suggesting there was different ploidy in or between sections [84]. Furthermore, polyploidization has also caused non-monophyletic species in the genus Aconitum [85][86][87]. Phylogenetic genetic analysis revealed tetraploid P. egaliksensis (taxonomically ascribed to Sect. Armerina) were clustered to Sect. Aleuritia [88]. Polyploidization and hybridization between other natural sections will need extensive analyses to reveal whether these events are responsible for the non-monophyly of the phylogeny. Secondly, it may be a wrong taxonomy. Morphological classification possesses strong subjectivity based on taxonomists. Besides, phenotypes are susceptible to environmental effects, such as convergent evolution, leading to two distant species evolving similar morphological characteristics. The calyptrate capsule that occurred in the genera Pomatosace, Anagallis, Soldanella, and Bryocarpum in Primulaceae has been attributed to convergent evolution [89]. Furthermore, P. secundiflora has long been treated as a member of Sect. Sikkimensis based on the campanulate corolla, but phylogenetic analysis suggested it should be placed in Sect. Proliferae Pax [90]. Recently, P. subansirica (taxonomically ascribed to Sect. Sikkimensis) was found to belong to Gesneriaceae [91]. Thirdly, chloroplast capture through hybridization may occur in these sections. Chloroplast capture is common in plants and also has been detected in the genus Primula [92][93][94].
Richard (2003) recognized seven subgenera and thirty-eight sections in Primula worldwide [23]. Our phylogenetic tree was largely inconsistent with morphology-based taxonomy. Twenty-one sections were divided into three main clades in our results. Clade I contained major sections of subgen. Auganthus; only one section belonging to subgen. Carolinella was nested within subgen. Auganthus. Our phylogenetic analysis implied neither subgen. Auganthus nor subgen. Carolinella are monophyletic taxa, consistent with previous analyses [10,95,96]. The subgen. Carolinella with its calyptrate capsule can be distinguished from other subgenera. Previous research has implied that the calyptrate capsule might have evolved multiple times in Primula [96]. We agree with this view because P. calyptrata and P. kwangtungensis were separated from each other and distantly related. Previous studies have also been performed on Sect. Obconicolisteri. These results were inconsistent with our study, due to insufficient samples and not being combined with other sections in previous phylogenetic analyses [10,16,69]. Our results demonstrated that Sect. Obconicolisteri is not a monophyletic group. Three species of this section were clustered with Sect. Monocarpicae and Sect. Cortusoides. To make Sect. Obconicolisteri monophyletic, these three species should be excluded. P. dumicola and P. oreodoxa were treated as members of Sect. Obconicolisteri, but appear closer to Sect. Monocarpicae. Regarding the two species with calyx broadly campanulate and P. dumicola with calyx margin entire, these characteristics are more similar to Sect. Monocarpicae. Besides, P. calyptrata (taxonomically ascribed to Sect. Carolinella) and P. handeliana (taxonomically ascribed to Sect. Crystallophlomis) possess a distant phylogenetic relationship with their natural morphology sections. Phylogenetic analyses based on protein-coding genes were also performed by Ren (2018), implying that P. calliantha and P. woodwardii did not have the closest relationship [97]. Molecular data and morphological cluster analysis based on population samples will be needed to infer whether it is a wrong taxonomy, or can be explained by chloroplast genome capture, introgression, or other evolution events.
Clade II included subgen. Primula (Sect. Primula) and subgen. Aleuritia. The subgenus Aleuritia is a large group in Primula with 19 sections. Our results were largely congruent with previous research that subgen. Aleuritia is paraphyletic [98,99]. P. gemmifera (taxonomically ascribed to Sect. Aleuritia) is endemic to China, mainly distributed in southern Gansu, western Sichuan, and northeastern Tibet. In our results, P. gemmifera is grouped with Sect. Capitatae Pax and Sect. Denticulata and form a well-supported clade. The non-monophyly of Sect. Aleuritia is consistent with the conclusion made by Ren (2015) [100]. However, the monophyly of sections in Aleuritia is also supported by other research [88,101], though these results were due to insufficient samples and a lack of comparative analysis with relevant groups. Furthermore, the monophyly of Sect. Sikkimensis is well supported in our result.

Conclusions
Although chloroplast genomes have the properties of uniparental inheritance and low frequency of recombination, these features have prompted plastomes to possess highly conserved features. However, variations at the intraspecies level have still occurred. Sixteen plastome genomes of P. obconica subsp. obconica were identified to have different lengths, which were mainly caused by variations in the LSC and SSC regions. Highly genome-wide variations were detected. Most SNPs were detected in downstream and upstream gene regions (45.549% and 41.91%). Multiple sequence divergent regions and SSRs found in this study would provide useful information for population genetic analyses in P. obconica. The phylogenetic tree showed that P. subsp. obconica can be split into two independent clades. The taxonomy of the P. obconica complex was redefined in this text. We elevated two subspecies of P. obconica to separate species. Our phylogenetic tree was largely inconsistent with the previous morphology-based taxonomy. Twenty-one sections of Primula were mainly divided into three clades. The monophyly of Sect. Auganthus, Sect. Minutissimae, Sect. Sikkimensis, Sect. Petiolares, and Sect. Ranunculoides were well supported in the phylogenetic tree. The Sect. Obconicolisteri, Sect. Monocarpicae, Sect. Carolinella, Sect. Cortusoides, Sect. Aleuritia, Sect. Denticulata, Sect. Proliferae Pax, and Sect. Crystallophlomis were not a monophyletic group. Several sections were represented by only one species, the phylogenetic relationships of these sections will need a larger number of samples in future studies.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13040567/s1, Figure S1: The presence of simple sequence repeats (SSRs) in the P. obconica subsp. obconica chloroplast genomes. (A) The number of repeats types; (B) Frequency of variations SSRs types, Figure S2: Comparison of 5 subspecies of P. obconica plastome genomes using mVISTA, Table S1: Details of accessions used in the study, Table S2: List of genes in the plastome genomes of P. obconica subsp. obconica, Table S3: The lengths of introns and exons for the splitting genes in P. obconica subsp. obconica, Table S4: SNPs identified in P. obconica subsp. obconica plastomes, Table S5: SNP annotation based on 16 plastome genomes of P. obconica subsp. obconica using SnpEff,

Conflicts of Interest:
The author declare no conflict of interest.