In total, 23 species of Swertia L. were selected, four of which were sequenced using Illumina sequencing; the remaining 19 sequences were obtained from GenBank. Fresh young leaves of S. tetraptera, S. franchetian, S. przewalskii, and S. bifolia were sampled from Mengyuan county (101.32′ E, 37.62′N, 3,208 m), Huangzhong county (101.63′ E, 36.57′N, 2,510 m), Qilian county (99.61′E, 38.83′N, 3,234 m), and Qilian county (102.22′E, 37.45′N, 3,135 m), respectively, all in the Qinghai province of China. Voucher specimens were deposited in the QTPMB (Qinghai–Tibetan Plateau Museum of Biology) with the voucher numbers QHGC-2011, QHGC20190821, QHGC-2013, and QHGC-2014, respectively. The leaves were dried and preserved in a silica gel.

The improved cetyltrimethylammonium bromide method was used to extract the total DNA of Swertia L. plants (Doyle, 1991). Agarose gel electrophoresis and a NanoDrop 2000 microspectrophotometer were used to measure the purity and concentration of the DNA. After Illumina PE library was constructed, high-throughput sequencing was completed by Beijing Biomarker Technologies Co., Ltd. Moreover, 150bp paired-end sequencing was performed using Illumina HiSeq (TM) 2000. Raw sequencing data were transformed into sequenced reads (raw data) by performing a base calling analysis of the raw image files. Raw reads data obtained by sequencing were filtered using ngSQCToolkit_v2.3.3 software (Patel and Jain, 2012) to remove low-quality regions and obtain clean reads. The results were then stored in the FASTQ format.

Chloroplast genome assembly was performed using the iterative organelle genome assembly pipeline (Bakker et al., 2016). The chloroplast genome of S. mussotii (NC_031155) was used as the reference sequence. SPAdes v3.6.1 software was used for ab novo splicing under default parameters and to generate a series of contigs (Prjibelski et al., 2020). Contigs larger than 1,000 bp were used for chloroplast genome assembly. Complete chloroplast genome sequences were constructed by matching and linking contigs (Kearse et al., 2012) and filling the gaps after assembly using second-generation sequencing technology.

The plastid genome annotator tool was used for the functional annotation of Swertia L. chloroplast genomes; the start codon, stop codon, and other problematic sites in the annotation result were adjusted manually (Qu et al., 2019; Tian et al., 2021). The annotated chloroplast genome data were exported in Gb format, and the chloroplast genome maps of the four Swertia L. species were drawn using OGDRAW (Marc et al., 2013) software. The sequence data and gene annotation information were then uploaded to the NCBI database. The GenBank accession numbers were NC_056357 (S.franchetiana), ON164641 (S.tetraptera), ON017794 (S.przewalskii), and ON018645 (S.bifolia).

We used CodonW1.4.2 software to confirm the relative synonymous codon usage (RSCU) and amino acid usage frequency.

The chloroplast DNA rearrangement analyses of the 23 Swertia L. species were carried out using Mauve alignment (Darling et al., 2004). To show interspecific variation, after annotating the files using Python 3.10.1, the chloroplast genomes of another 22 species of Swertia L. were compared using the online software mVISTA (Frazer et al., 2004) and S. wolfgangiana as a reference genome. Variations were detected using the Shuffle-LAGAN model. The percentages of variable characters in the coding and non-coding regions were calculated using the method developed by Zhang et al. (2011). IRscope software (Amiryousefi et al., 2018) was used to visually analyze the contraction and expansion of the four boundaries of the 23 species of Swertia L.

The online software REPuter (Kurtz et al., 2001) was used to detect repeats in the chloroplast genome, such as forward (F), reverse (R), complementary (C), and palindromic (P). The minimum repetition was set to 30 bp and minimum repetition sequence length distance to 3. In addition, the online program Tandem Repeats Finder was used to detect tandem repeats (Benson, 1999). MISA software (https://pgrc.ipk-gatersleben.de/misa/) was used to predict simple sequence repeat (SSR) in chloroplast genome, and the parameters were set as follows: mononucleotide unit repetition number ≥10; dinucleotide unit repetition number ≥5; trinucleotide unit repetition number ≥4; and tetraconucleotide, pentanucleotide, and hexanucleotide unit repetition number ≥3 (Beier et al., 2017).

In this study, 23 species were used to construct a phylogenetic tree based on Bayesian inference (BI) (Ronquist and Huelsenbeck, 2003), using Gentianopsis paludosa (NC_050656) as the outgroup. Mafft v7.205 software was used to compare the sequences and remove irregular sequences at both ends (Kazutaka and Standley, 2013). Before building the BI tree, PAUP and MrModeltest were jointly run through MrMTgui. The Akaike information criterion results showed that the best model for BI analysis was GTR + I + G, with a random tree as the starting tree. Starting with four Markov chains, that is, three hot chains and one cold chain, we saved one tree every 100 generations, calculated 9,000,000 generations, discarded the first 25% preheated (Burn-in) trees, and used the remaining trees to calculate the Bayesian posterior probability (PP) of the consistent tree and each branch.

Based on the obtained chloroplast genome sequences, the divergence times of Swertia L. species was estimated using the Markov Monte Carlo algorithm (MCMC) molecular sequence Bayesian analysis in BEAST V1.7 (Drummond et al., 2012). First, BEAUti in the software package of BEAST was used to set the parameters of the sequence file in the Nexus format, and the optimal nucleotide substitution model was GTR + I + G, which was selected by MrModeltest. The uncorrelated relaxed clock method was used for the branch lengths with a Gama distribution. Due to the lack of fossil evidence for Swertia L. plants, the time was set at 15 Ma (million years), which was from the published literature (Chassot et al., 2001; Von Hagen and Kadereit, 2002; Cao et al., 2021), and the standard variance was 1.0. After a burn-in of 10,000,000 steps, all of the parameters were collected once every 1,000 steps up to 1,00,000,000 Markov chain Monte Carlo (MCMC) algorithm steps. Then an XML format file was generated. The XML format file was imported to BEAST software. The convergence of the MCMC results was detected by using the Tracer v 1.5 program to check that the chain was balanced; we then used the Tree Annotator v 1.7.5 program to obtain the best tree merging and Figtree v 1.4.4 (Rambaut, 2018) was used to view the resulting tree.

The chloroplast genome lengths of the Swertia L. species ranged from 149,036 bp to 153,691 bp, with an average length of 152,385 bp (Table 1). S. bimaculata had the longest chloroplast genome, differing from other species in Swertia L by 0.06–4.715 kb. As can be seen from the comparison of chloroplast sections, such differences mainly occurred in the large single-copy (LSC) and IR regions. The chloroplast genome length of angiosperms is generally 115–165 kb and that of Gentianaceae is 137–154 kb, which is consistent with the length characteristics of angiosperms and Gentianaceae (Li et al., 2018; Dong et al., 2020). Compared with other genera of Gentianaceae, the average chloroplast genome length of Swertia L. was similar to that of Halenia (153 kb), but shorter than that of Paedera (154 kb) (Dong et al., 2020). The chloroplast genomes of the Swertia L. species contained two reverse repeats, IRa and IRb, which divided the whole genome into four parts; the remainder comprised LSC and small single-copy (SSC) regions (Figure 1). The chloroplast genomes of the Swertia L. species had the ring tetrad structure typical of angiosperm chloroplast genomes (Palmer, 1985), which made the chloroplast genome highly conserved. The lengths of the LSC regions varied from 80,432 bp to 84,156 bp, with a total GC content of 32.18%–36.35%. The GC content of the SSC region was 31.25%–33.66%, and the total length ranged from 17,887 bp to 18,395 bp. The pair of IRs had a length range of 25,069–25,890 bp and GC content of 42.16%–44.38% (Table 1).

Similar chloroplast DNA GC compositions were found in all of the Swertia L. species (Table 1), demonstrating high species similarity. The IR regions had a higher GC content than the LSC and SSC regions; this has also been reported in other plants (Choi and Park, 2015; Guo et al., 2020). The IR region contained rRNA and tRNA genes, which accounted for the high DNA GC content of this region (Doorduin et al., 2011; Asaf et al., 2017; Shen et al., 2017).

Most of the chloroplast genomes of angiosperms encode 74 proteins, but some genes have been captured, rearranged, and lost across different families, genera, and species (Millen et al., 2001; Kim et al., 2009). The results of our study showed that S. bimaculata, S. cordata, S. diluta, S. erythrosticta, S. franchetian, S. kouitchensis, S. leducii, S. macrosperma, S. mussotii, S. punicea, S. souliei, S. vertickllifolia, and S. wolfgangiana had 133 genes comprising 87 protein-coding genes, 38 tRNA genes, and eight rRNA genes. S. cincta, S. dichotoma, S. nervosa, and S. pubescens lacked the rps16 gene found in the chloroplast genomes of other species of Swertia L. Thus, these four chloroplast genomes consisted of 132 genes. The ycf15 gene in the two reverse repeats was lost in S. przewalskii and S. bifolia, implying that their chloroplast genomes contained 131 genes. Our result was different from the previous result obtained for Gentianaceae (Dong et al., 2020), which showed that the chloroplast genome of Gentianaceae had 67–80 protein-coding genes, 30 tRNA genes, and four rRNA genes. This difference mainly arose due to gene deletion between genera. For example, the loss of ndh genes, including ndhA, ndhC, ndhG, ndhH, ndhI, ndhJ, and ndhK, was common to all Gentianaceae species. In addition, four pseudogenes (ψrps16, ψrps19, ψinfA, and ψycf1) were present in the chloroplast genomes of the Swertia L. species. Previous studies have shown that Gentianaceae plants generally have the same four pseudogenes; our results confirm these previous observations. The ψinfA pseudogene likely appeared due to transfer or loss during species evolution (Millen et al., 2001; Zhou et al., 2016). The appearance of the ψrps19 and ψycf1 pseudogenes is likely due to their location at the boundary of the chloroplast gene region, which experiences a boundary effect (Li et al., 2018). The second missing exon in the ψrps16 pseudogene was first detected in Gentiana macrophyllum (Ni et al., 2016) and in non-parasitic species of the Chrysanthemum branch (APG IV). Since then, the ψrps16 pseudogene has been detected in the chloroplast genomes of several Gentianaceae members, the structures of which are similar across all species.

The functions of the major genes in the chloroplast genomes of Swertia L. were roughly classified into three categories (Table 2): chloroplast self-replication–related genes, photosynthesis-related genes, and other genes (Saski et al., 2005). Genes related to photosynthesis and self-replication accounted for the majority of the chloroplast genome.

Further analysis of the chloroplast genes of Swertia L. showed that they were similar to those of other plants and that most did not contain introns (Du et al., 2018; Guo et al., 2020). In this study, only 16 genes (rps12, trnK-UUU, atpF, rpoC1, ycf3, trnL-UAA, trnV-UAC, clpP, petB, petD, rpl16, rpl2, ndhB, trnI-GAU, trnA-UGC, and ndhA) in the chloroplast genomes of Swertia L. contained introns, and all of them contained one intron except for the clpP and ycf3 genes, which had two introns (Table 2). The rps12 gene in the chloroplast genomes of Swertia L. experienced trans-splicing, in which the 3′ end was in the IR region and 5′ end was in the LSC region. This phenomenon has been observed in the majority of other land plants (Du et al., 2018).

The preference of 59 synonymous codons was evaluated using RSCU (Wu et al., 2007). Based on the statistical analysis, the number of codons in the Swertia L. species varied from 49,696 to 512,30. Leucine (Leu; 4,988–5,394 codons), isoleucine (Ile; 3,730–4,277 codons), and phenylalanine (Phe; 3,498–3,641 codons) were the three amino acids with the highest coding rates in the Swertia L. species chloroplast genomes. Only 663–719 codons encoded tryptophan (Trp), which had the lowest coding rate among all of the amino acids (Supplementary Table S1).

Repetitive sequences are the main sources of duplication, deletion, and rearrangement in the chloroplast genome (Li and Zheng, 2018). In this study, four kinds of repetitions were counted: forward, palindromic, tandem, and reverse. The results showed that the distributions and numbers of repeats in the 23 chloroplast genomes were similar and conserved (Figure 2; Supplementary Table S2). Tandem units were the most repeated type (605), followed by forward (260), palindromic (209), and reverse repeats (4) (Figure 2C). There were interspecific differences in the tandem repeats, but the ratio of forward to palindromic repeats was about 1:1. Reverse repeats only existed in S. cincta, S. leducii, and S. macrosperma. The lengths of the repeat units were mainly 8–39 bp (Figure 2A). The majority of repetitive sequences were scattered across intergenic or intronic regions, with only a few distributed across gene regions such as ycf3, ycf2, ndhE, psaB, accD, petB, ndhA, psbA, accD, rps18, rps16, psbK, clpP, ycf1, atpH, and rps2 (Supplementary Table S2). S. bimaculata had the most repeat sequences (76) of all the analyzed Swertia L. species, followed by S. leducii (67); S. bifolia had the fewest repeat sequences (34) (Figure 2B).

As a classical molecular marker, simple repeat sequence (SSR) has been widely used in the analysis of population genetic evolution. We analyzed the simple repeat sequence (SSR) in the chloroplast genomes of 23 species of Swertia L. and the result showed that the numbers of SSR ranged from 35 to 61. S. tetraptera had the most SSRs (61) and S macrosperma had the fewest SSRs (38). Moreover, the numbers and types of mononucleotide, dinucleotide, trinucleotide, tetranucleotide, pentanucleotide, and hexanucleotide repeats were also different in the 23 species of Swertia L. (Figure 3; Supplementary Table S3). Mononucleotides were the most common repeat type. The proportion of mononucleotides in all SSRs ranged from 50.00% to 82.22% in 23 species of Swertia L. This finding is in accordance with the previous observation (Kuang et al., 2011). In total, 70 dinucleotides were detected in 23 species, which were AT/TA, accounting for 3.23%–10.53% of the SSRs. In total, 76 trinucleotides and 133 tetranucleotides were found in the 23 complete cp genomes. A total of 15 pentanucleotides were discovered in chloroplast genes of 23 species in Swertia L. Only S. cordata (2), S. dichotoma (1), S. franchetiana (1), S. mussotii (1), S. nervosa (3), and S. tetraptera (3) had hexanucleotides. In addition, compound SSRs accounted for 2.17%–10.87% of the 23 genomes. The richness of SSRs and the count of SSRs were different within Swertia L. thus these may be helpful molecular marker for species identification. However, adopting SSRs to clarify ecological and evolutionary processes has yet to be fully implemented (Ebert and Peakal, 2009). The results of this study will provide a basis for the study of chloroplast SSR markers in the future and lay a foundation for the study of the genetic relationship and diversity of this genus.

Oligonucleotide repeats are widely found in the plastome (Ahmed et al., 2012, 2013; Abdullah et al., 2021). These repeats have an effect on generating mutations and have been suggested as a proxy for mutational hotspots (Ahmed et al., 2012; Ahmed et al., 2012; Abdullah et al., 2020; Abdullah et al., 2021). Abdullah et al. (2020) proposed that the co-occurrence of repeats with substitutions was up to 90%, whereas 36%–91% co-occurrence was found at the genus level. In the present study, 10 highly polymorphic loci were found. Among these, five loci belong to the regions where repeats are present, including psaA-ycf3 and rps15, which showed the highest incidence of polymorphisms. Here, our findings support the use of repeats as a proxy, and this approach may also be helpful for the identification of suitable polymorphic loci for phylogenetic inference of other taxonomically complex genera. This approach is promising since the plastome of a single species can be used to identify polymorphic regions. Repeated coding regions and IR regions need to be avoided, however, due to the purifying selection pressure of protein-coding genes (Henriquez et al., 2020) and the fact that copy-dependent repair mechanisms (Zhu et al., 2016) lead to low rates of mutation.

The chloroplast genomes of the 23 Swertia L. species were relatively conserved, with four parts of the genomes being arranged in consistent sequences (Figure 4) and no rearrangement found in gene organization after verification (Figure 5). Moreover, there was a higher degree of variation in non-coding regions than in the coding regions of the chloroplast genome of Swertia L. In the non-coding regions, the percentage of variations ranged from 13.14% to 81.84% (Figure 6B), averaging 49.02%, whereas in coding regions, the percentage of variations ranged from 0.35% to 31.27%, averaging 9.10% (Figure 6A). The SSC region variability of the 23 species in Swertia L. was higher than that of the LSC and IR regions in both coding (7.96%, 2.19%, and 17.16% for LSC, IR, and SSC regions, respectively) and non-coding regions (49.00%, 42.44%, and 54.23% for LSC, IR, and SSC regions, respectively). The degree of variation was lowest in the IR region, indicating a high degree of conservatism. These results were consistent with those obtained for other angiosperms (Dong et al., 2013; Guo et al., 2020). In addition, some genes (rpoC1, ccsA, ndhI, ndhA, and rps15) exhibited higher variability than other genes in the 23 species of Swertia L. Some of the non-coding regions with high sequence divergence were trnH-GUG-psbA, psaA-ycf3, cemA-petA, ycf15-trnL-CAA, and ccsA-ndhD. These genes and hotspot regions can either be used in phylogenetic analyses or serve as potential DNA molecular barcodes (Zhang et al., 2011; Maier et al., 1995; Diekmann et al., 2009).

The chloroplast genome has two IR regions, which form four boundaries with LSC and SSC regions: IRb-LSC, IRb-SSC, IRa-LSC, and IRa-SSC. When the ancient genome evolved, the IR boundary expanded and contracted, causing some genes to enter IR regions and some to enter the single-copy regions, with different levels of sequence replication at each species boundary. As can be seen from Figure 7, the four boundaries of the chloroplast genomes of the Swertia L. species were relatively well-conserved. The rps19 gene spanning the LSC and IRb regions was present at the IRb-LSC boundary in all 23 Swertia L. chloroplast genomes. This gene was mainly located in the LSC region at the same bases, except in S. cordata (85), S. cinata (118), and S. pubescens (118). The IRa-LSC boundaries in most of the Swertia L. chloroplast genomes occurred between the rps19 pseudogene in the IRa region and the trnH gene in the LSC region; however, the rps19 pseudogene was absent in S. bifolia, S. przewalskii, S. nervosa, and S. multicaulis. The IRb-SSC boundaries in the Swertia L. chloroplast genomes varied greatly. This boundary was located in the overlapping region of the ycf1 pseudogene and ndhF gene in 11 Swertia L. chloroplast genomes, with the IRb-SSC boundary in six Swertia L. chloroplast genomes crossing the overlap region and extending 5–100 bp to the ndhF gene. The ycf1 pseudogene in eight Swertia L. chloroplast genomes was present in the IRb region, along with a terminal from the IRa-SSC border. In addition, the ycf1 pseudogene was lost in the IRb-SSC boundaries of the chloroplast genomes in S. tetraptera, S. nervosa, and S. multicaulis (Figure 6). The IRa-SSC boundary was located in the ycf1 gene in all of the species, but the length of the ycf1 gene fragment in the IRa region differed to some extent and ranged from 988 bp to 1,004 bp. The length of this fragment was about 5,400 bp in most Swertia L. species, except for S. nervosa and S. souliei. The ycf1 gene in the S. nervosa chloroplast genome was present in the SSC region, with a terminal 126 bp from the IRa-SSC border. The total length of the ycf1 gene in the S. souliei chloroplast genome was 1,013 bp, with only 10 bp located in the SSC region. The sliding of the IRa-SSC and IRb-LSC boundaries in the chloroplast genomes of vascular plants generally occurs in different genera or even within the same genus, resulting in large variations in chloroplast genome length across different plants. The IRb-LSC boundaries of the Swertia L. species were largely located within the rps19 gene and, as mentioned earlier, the IRa-LSC boundary was located between the rps19 gene of the IRa region and the trnH gene of the LSC region. However, in monocotyledon plants such as those in the Orchidaceae and Poaceae families, the boundaries are extended and the rps19 and trnH genes are located in the IR regions (Tang et al., 2011; Hu, 2020). Both genes changed from one to two copies, whereas in barley and sorghum, boundary shrinkage occurred, resulting in two copies of the rps19 and trnH gene in the LSC region (Tang et al., 2011). The IRb-SSC boundary was located in the ycf1-ndhF overlap region in 11 Swertia L. species, which is consistent with the observations from many species of cruciferous plants (Li et al., 2017), for example, Aethionema grandiflorum, Arabidopsis thaliana, Barbarea verna, Brassica napus, Cakile arabica and so on. The extension of the IRb-SSC boundary to the ndhF gene due to boundary expansion has also been detected in Arabidopsis thaliana (Tang et al., 2011), in which the IRa-SSC boundary is located in the ycf1 gene. In A. thaliana, the fragment lengths of the ycf1 gene in the SSC and IRa regions are different due to either contraction or expansion of the boundary. In rice, wheat, maize, and other plants, this boundary is located on the ndhH gene, further indicating that the boundary between dicotyledons and monocotyledons is quite different (Melodelima et al., 2013).

People have different views on the mechanism of contraction and expansion of the IR region (Guo et al., 2020). DNA double-strand breaks (DSBs) are currently considered the main molecular mechanism underlying IR expansion. As large contractions in the IR region are rare, the DSB theory may also underlie IR region contraction.

The maximum likelihood and Bayesian methods were used to construct phylogenetic trees for the chloroplast genomes of the 23 Swertia L. species. The topological structures of the phylogenetic trees obtained using the two methods were similar (Figure 8). Phylogenetic analysis showed that all 23 species of Swertia L. in conjunction with those of G. paludosa, formed a well-supported clade, indicating that the genus Swertia L. was not monophyletic. This result is supported by previous studies (Chassot et al., 2001; Struwe and Albert, 2002; Von Hagen and Kadereit, 2002; Favre et al., 2010; Xi et al., 2014; Cao et al., 2021). In addition, the well-supported clade was divided into two major clades (A and B), corresponding to the subgen. Swertia (A) and subgen. Ophelia (B).

Seven species in Swertia L. were clustered into a subgen. Swertia phylogenetic tree at the base, which showed a close genetic relationship. Ho et al. (1994) believed that this subgenus was a relatively primitive group of Swertia L. consisting of a perennial herb with ancestral traits such as a single stem and unbranched large flowers. Within the same clade, the four species of sect. Swertia (S. souliei, S. bifolia, S. wolfgangiana, and S. erythrosticta) formed a single clade (A1) and two species (S. cordata and S. nervosa) of sect. Ophelia formed another clade (A2). These two were sister branches, further supporting the division of these groups by Ho et al. (1994). Clade B had two branches: B1 and B2 subclades. The B1 subclade contained S. bimaculata, which belonged to sect. Ophelia. This clade also included an S. dichotomy–S. tetraptera branch. S. bimaculata and S. dichotomy–S. tetraptera were sisters. The plants in these two branches were closely related (100 bootstrap), indicating a common ancestor. The B2 subclade contained sect. Ophelia, sect. Platynema, sect. Poephila, and sect. Macranthos. In this subclade, S. leducii was differentiated first and located at the base. Furthermore, two parallel branches were then isolated: sect. Ophelia and sect. Platynema; sect. Poephila and sect. Macranthos. S. multicaulis, from subgen. Poephila, and S. verticillifolia, from sect. Macranthos, were first clustered into a small clade and then into a large clade with the three species of sect. Platynema and one species of sect. Ophelia. This differed from the morphological classification. Sect. Platynema was at the top of the B2 subclade, indicating that it was located in a comparable evolutionary position of the phylogenetic tree of Swertia L. The clustering results partially validated the results obtained by Ho et al. (1994), who showed that sect. Platynema and sect. Kingdon-Wardia (Marq.) were the most evolved groups of the genus and characterized by extremely enlarged filaments at the base, a single glandula in each corolla lobe, and diminished tassels. In the present study, sect. Kingdon-Wardia (Marq.) was not included in the phylogenetic tree, making it impossible to show its systematic position. However, sect. Platynema and sect. Kingdon-Wardia (Marq.) were clustered together and located in the same relative evolutionary branch of Swertia L. in a study by Xi et al. (2014). From what has been discussed before, the division of the two subgenera (subgen. Swertia and subgen. Ophelia) and five sections (sect. Ophelia, sect. Platynema, sect. Poephila, sect. Swertia, and sect. Macranthos) is partially supported by molecular data. However, the systematic positions of other sections and species in Swertia L. derived from molecular data differed from the morphological classification. Inconsistencies between different data types, specifically between morphological and molecular data, remain a major problem of systematics (Lee, 2001). Such inconsistencies have been reported and discussed for many plant and animal groups, such as Rubiaceae, Loganiaceae, Isothecium, and Dendrocolaptinae (Bremer and Struwe, 1992; Irestedt et al., 2004; Draper et al., 2007). Pisani et al. (2007) argued that despite the widespread inconsistencies between morphological and molecular data, both data types were equally important in estimating phylogenetic relationships and that molecular data could not be considered more reliable. The results of this study were roughly equivalent to those of previous studies that used different gene fragments and species to examine the phylogeny of Swertia L. indicating a conflict between the morphological classification system and molecular data, which can be explained from the perspective of evolution. The formation of new species is a slow process, usually occurring over thousands of years. Variations due to natural selection and genetic drift become fixed in a group, driving the formation of new species that eventually differ from two recent common ancestors, that is, species derived from two recent common ancestors, both morphologically discontinuous and reproductively isolated, are monophyletic (Liu, 2016). Driven by the drastic changes in the geology and climate of the Qinghai–Tibetan Plateau, the ancestors of Swertia L. evolved rapidly and showed abundant morphological diversity, such as in the shape and length of the corolla and number and location of nectaries, nectary appendages, and corolla throat appendages. However, this taxon has not accumulated enough sequence variation for a molecular phylogenetic analysis over a relatively short period of time. Moreover, mutations in gene sequences have not been fixed in the population by genetic drift. In addition, the uniparental inheritance of the plastome may also confound phylogenetic inference. Previous studies have shown that the phylogeny based on plastome and mitochondria sequences contradicted with nuclear due to uniparental inheritance of these genomes (Vargas et al., 2017; Abdullah et al., 2021). Therefore, more genetic markers (nuclear) and more taxa of Swertia L will be needed to further explore the phylogenetic relationships in this genus.

Tracer v 1.5 was used to check the analysis values of each parameter, and it was shown that the number of MCMC iterations calculated by BEAST had met the effective sample size (ESS), which was greater than 200. The BEAST analysis was based on the phylogenetic trees of chloroplast genomes of 23 species of Swertia L.(Figure 9), and the numbers at each branch node of the phylogenetic tree were the divergence times (Ma) of the corresponding groups. The result showed that the estimated divergence between Swertia L. and Gentianopsis occurred at 29.60 Ma. We therefore inferred that Swertia L formed at 29.60 Ma, corresponding to the early Miocene of the Tertiary. Meanwhile, the divergence between subgen. Swertia and subgen. Ophelia appeared at 14.69 Ma. In addition, the estimated divergence time in 23 species of Swertia L. was between 12.40 and −0.05 Ma. The formation of S. franchetiana, S. mussotii, S. punicea, S. kouitchensis, S. diluta, S. pubescens, S. cincta, S. dilatata, S. hispidicalyx, S. souliei, S. bifolia, S. wolfgangiana, and S. przewalskii were at 0.05–1.33 Ma (the Quaternary), and S. macrosperma, S. erythrosticta, S. nervosa, S. cordata, S. tetraptera, S. dichotoma, S. bimaculata, S. verticillifolia, S. multicaulis, and S. leducii were formed at 2.72–12.40 Ma (end of Tertiary).

In the present study, the formation of Swertia L. was dated back to 29.60 Ma, which was slightly earlier than other studies (Chassot et al., 2001; Von Hagen and Kadereit, 2002; Cao et al., 2021). Geologic evidence demonstrated that the turn of the Oligocene and Miocene was a crucial period of the tectonic evolution of the Qinghai–Tibetan Plateau(QTP), the central part of the QTP rose to a height of nearly 3,000 m in the Early Miocene, the cooling effect made by QTP uplift resulted in the transition of QTP from tropical and subtropical environment to a warm and cool environment consistent with the temperate climate, and the further development of herbaceous plants began in the Early Miocene (Deng et al., 2019). During this period, a primitive group of Swertia L. plants appeared, represented by subgen. Swertia L., which was characterized by perennial herbs, single stems, unbranched, and large but few flowers (Cao et al., 2021).

During the following 20 Ma to 10 Ma period, the QTP was further uplifted, and the Himalayan mountains and Tianshan Mountains were significantly elevated, which strongly changed the atmospheric circulation. Meanwhile, the global temperature decreased from the optimum temperature in the middle Miocene of the third century, resulting in a cool and dry climate (Miao et al., 2012). During this period, Swertia L. plants appeared as annual herbs with strongly branched stems and many small flowers, represented by subgen. Ophelia. The new taxa produced a large number of seeds during their life cycle and thus were better adapted to changing environments (Cao et al., 2021). When the climate was suitable, the new species gave rise to a large number of offspring, which has the potential for a great deal of variation.

Since 10 Ma, the QTP has been further uplifted in the late Miocene and Pliocene periods, and the Himalayas have blocked almost all the warm and wet air masses from the Indian Ocean, and the QTP has become cold and arid. Since 4 Ma, the QTP has been affected by Quaternary glaciation (Li et al., 1999; Mulch et al., 2006). The complex landform and rapidly changing climate resulted in many isolated small populations of Swertia L. which underwent radiation differentiation due to differentiated selection and random factors, forming new species adapted to local environment in a relatively short period of time. This process of radiation differentiation eventually led to the diversity of Swertia L. plants today. In this study, 13 species of the 23 Swertia L. species were formed at Quaternary. This group is the most richly differentiated and most widely adapted in Swertia L. with distribution in both plateau and plain.