Comparative analysis of complete plastid genomes from Lilium lankongense Franchet and its closely related species and screening of Lilium-specific primers

Plant materials and DNA extraction

Leaf samples were collected from each field site (Table S1) and then immediately dried with silica gel to preserve them for DNA extraction. The total genomic DNA was extracted via the NuClean Plant Genomic DNA Kit (ComWin Biosciences, Jiangsu, China) following the manufacturer’s instructions.

Plastid genome sequencing, assembling, and annotation

All complete plastid genomes were sequenced using the Illumina Novaseq 6000 Platform (Novagene, Beijing, China) with an average paired-end read length of 150 bp. In order to assess the quality of sequenced raw reads, we used the FastQC (Andrews, 2014) v0.11.7. The plastid genome related reads were filtered via mapping after they were assessed for quality and all the raw reads were associated with the previously reported plastid genome sequences in Liliaceae using Bowtie 2 (Langmead & Salzberg, 2012) v2.3.4.3. The published sequence with the highest overall similarity to the reads was used as the seed sequence for further assembly using NOVOPlasty v2.7.2 (Dierckxsens, Mardulyn & Smits, 2017) to obtain a draft sequence. The draft sequence was imported into Geneious v11.0.4 software (Kearse et al., 2012) and was compared with the optimal reference sequence and mapped reads to check and correct for mismatches. The annotation of the plastid genome was also performed in Geneious. Every annotated gene was manually edited for start/stop codons and intron/exon boundaries to correct errors and ensure the accurate annotation of the genome. Finally, circular plastid genome maps were drawn using OGDRAW (Lohse et al., 2013).

Simple sequence repeats (SSRs) and IR borders analysis

Perl script MISA (Thiel et al., 2003) was used to detect SSRs loci in the six plastid genomes of Lilium with thresholds of 10, 5, 4, 3, 3, and 3 for mono-, di-, tri-, tetra-, penta-, and hexanucleotides, respectively. The IR regions were discovered using the plugin Repeat Finder v1.0 in Geneious and were manually checked. The IR plot of the six Lilium species plastid genomes was obtained using an IRscope online program to compare the IR borders and neighboring genes of the six Lilium species plastid genomes (Amiryousefi, Hyvönen & Poczai, 2018).

Lilium-specific primers

We generated multiple sequence alignments of the six newly assembled sequences and all available plastid genome sequences in Lilium (downloaded from NCBI in October 2019) using MAFFT v1.4.0 under the automatic model selection option (Katoh & Standley, 2013). We developed primers specific for the highly mutational regions in Lilium using Geneious after alignment (Kearse et al., 2012). The primers were synthesized by Sangon Biotech, China. 12 primers pairs were designed (Table S2) through a series of experiments and two segments of gene spacers were selected with the goal of identifying those primer pairs that would produce the best results. Primer sequences and polymerase chain reaction (PCR) amplification results were summarized in Table S3. PCR of trnC-petN was conducted as follows: initial denaturation at 94 °C (3 min) followed by 33 cycles with denaturation at 94 °C (40 s), annealing at 60 °C (50 s), extension at 72 °C (1 min). PCR for rpl32-trnL was: initial denaturation at 94 °C (3 min), followed by 33 cycles of denaturation at 94 °C (40 s), annealing at 54 °C (45 s) and elongation at 72 °C (1 min). All reactions ended with a final elongation at 72 °C (10 min) followed by a holding step at 4 °C. All PCR products were sequenced in both directions by Sanger sequencing (Sanger & Coulson, 1975) (performed at Sangon Biotech, China), and then assembled and edited using Seqman v7.1.0 software (DNAstar package; DNAStar Inc., Madison, WI, USA) (Burland, 2000) for obtaining consensus sequences. PCR products of these intergenic spacer regions were obtained for 11 additional accessions belonging to three Lilium species, including those from geographically remote populations. L. lankongense was represented by five populations (from Sichuan, Yunnan, and Tibet), L. duchartrei was represented by four populations (from Sichuan), and L. lophophorum was represented by two populations (from Sichuan and Yunnan); each of these species had one individual per population.

Phylogenetic analysis

We combined the 27 published plastid genome sequences (downloaded from NCBI in October 2019) and ITS sequences from the NCBI and the six plastid genome sequences and corresponding ITS sequences measured by this study (GenBank accession number: MT260888 –MT260893, leaf samples shown in Table S1). Two sequences in Fritillaria were treated as the outgroups. We adopted the maximum likelihood (ML) and Bayesian inference (BI) methods to analyze the phylogenetic relationship of Lilium. The universal primers ITS4 and ITS5 (White et al., 1990) were used to amplify the ITS according to the standard PCR protocols of Gao et al. (2013). A total of three data sets, namely 33 complete plastid genome sequences, all shared protein-coding genes (CDS) (only containing one IR). The 33 plastid genome sequences and ITS sequences of the species corresponding to these 33 plastid genome sequences were used to reconstruct the phylogenetic tree of Lilium. We extracted the shared CDS from 33 plastid genome sequences as follows: all CDS were extracted from all 33 sequences in the software Geneious v11.0.4, and shared 71 CDS (excluding additional copies in the IR) of the 33 species were selected and manually sorted. The CDS belonging to one species were concatenated to generate 33 sequences for phylogenetic analysis. The three data sets were aligned by MAFFT v1.4.0 under the automatic model selection option, trimmed via trimAl v1.2 (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009) with parameters Trimal -in *.fas -noallgaps -fasta -out *.fas, and were compiled into three alignment matrices. The three alignment matrices were used to search for the best-fit substitution model using ModelFinder plugin in PhyloSuite v1.2.1 software before ML and BI analyses were conducted (Zhang et al., 2020). ML analysis based on the SYM+G (ITS) and GTR + I + G (plastid genome sequences and shared concatenated 71 CDS sequences) model was conducted using RAxML v8.2.8 (Stamatakis, 2014) with 1000 bootstrap replicates. BI analysis based on the GTR + I + G model was conducted using the Mrbayes v3.2.6 plugin in PhyloSuite v1.2.1 software. The Monte Carlo Markov chains (MCMCs) were run 1 × 10⁸ generations and the first 30% of trees were discarded as burn-in. ML tree used Bootstrap support (BS) and BI tree used posterior probability (PP) to evaluate the feasibility of each branch.

Plastid features of Lilium species

The complete plastid genome of the six species in Lilium was deposited in GenBank. These plastid genomes ranged from 152,307 bp (L. stewartianum) to 152,611 bp (L. lankongense) in length, with the minimum and maximum differences being 15 and 304 bp, respectively (Table 1 and Fig. 2). All six plastid genomes displayed a typical quadripartite structure, consisting of the LSC (81, 870–82, 150 bp) and SSC (17, 382–17, 542 bp) regions separated by a pair of IRs (26, 425–26, 577 bp). In the six Lilium plastid genomes, the overall GC content was 37%. Both contained 132 genes, with 86 protein-coding genes, 38 tRNA genes, eight ribosomal RNA genes (Table 1 and Table S3).

Simple sequence repeats (SSRs) analysis

The SSRs loci in the six plastid genome sequences of Lilium were detected using MISA perl script. We detected a total of 340 SSRs in the six Lilium plastid genomes, with different numbers of SSRs per species (Table 2). 53, 54, 57, 58, 63, and 55 SSRs were observed in Lilium lankongense, L. duchartrei, L. stewartianum, L. matangense, L. lophophorum, and L. nanum (Table 2, Table S4). L. lankongense (53 SSRs) and L. lophophorum (63 SSRs) had the lowest and highest number of SSRs, respectively. All SSRs were categorized as mononucleotide, dinucleotide, trinucleotide, tetranucleotide, pentanucleotide, and hexanucleotide (Table 2). The total number of mononucleotide repeats consisting of A or T bases was the largest compared to other types of microsatellites (Fig. 3A) and was greater than the sum of the other types. There were no hexanucleotide repeats found in the six Lilium species. In the total SSRs loci, most of the repeats were located in the LSC region, followed by the SSC region and IR regions (Fig. 3B). Additionally, 20 common SSRs were found among these six plastid genome sequences (Table S4), and the number and location of these SSRs were different.

IR expansion and contraction

The IR boundaries of the six Lilium species were shown in Fig. 4. The rps19, ndhF, ycf1, and psbA genes were distributed in the J_LB, J_SB, J_SA, and J_LA regions, respectively (Fig. 4). The order and types of genes were very conservative. Despite the similar lengths of these six species IR regions (ranging from 26,425 to 26, 577 bp) there is a slight difference in the IR expansion and contraction. The rps19 gene crosses the LSC and IR_B regions (J_LB), with 140–145 bp extension into the IR_B region. The J_SA line intersects the ycf1 gene and the SSC and IR_A regions are the same in Lilium stewartianum and L. matangense (4,327 bp in SSC and 1,226 bp in IR_B) but are different in other species (4,323 to 4,333 bp in SSC and 1,226 to 1,244 bp in IR_B) (Fig. 4). The psbA gene in the J_LA regions remained completely within the LSC region, 78 to 83 bp away from J_LA line. The changes in the boundary transformations between the LSC/IR and J_SA regions were relatively stable, whereas the J_SB regions varied among the six species to some extent. The J_SB line is located between ycf1 and ndhF. Compared to the six Lilium species, the J_SB region showed that ndh F expanded in three species (L. stewartianum, L. matangense, L. lophophorum), but contracted in two species (L. duchartrei, L. nanum) (Fig. 4).

Development of Lilium-specific primers

Two of the 12 primer pairs designed to sequence the plastid genome of the genus Lilium and tested on three species of Lilium, amplified consistently and easily produced good quality sequences. The two primer pairs, listed in Table 3 and Fig. 5, were tested on 11 individuals belonging to Lilium lankongense, L. duchartrei, L. lophophorum, including those from geographically remote populations (Table S5). Comparisons of intergenic spacer regions among L. lankongense, L. duchartrei, L. lophophorum plastomes provided evidence of high variability within trnC-petN, rpl32-trnL regions. The 620 bp trnC-petN and 793 bp rpl32-trnL regions were aligned for 11 individuals (Table S6, Supplementary Nucleotide alignment trnC-petN and Nucleotide alignment rpl32-trnL). We found nine, four, and one mutation sites in the trnC-petN region in L. lankongense, L. duchartrei, and L. lophophorum, respectively. Sequences in the rpl32-trnL region aligned data set were observed in the different number of variable sites in L. lankongense (10), L. duchartrei (4), L. lophophorum (1) (Table 4 and Table S6).

Phylogenetic analysis

We constructed the phylogenetic tree (containing 31 Lilium species, two Fritillaria as the outgroups) from the three data sets of the complete plastid genome sequences, the shared 71 CDS sequences, and ITS sequences using the maximum likelihood and Bayesian inference methods. The results are shown in Figs. 6 and 7, and Fig. S1. The phylogenetic trees obtained by ML and the BI methods had identical topologies. The topologies of the phylogenetic tree constructed by ML and BI methods were largely the same based on the complete plastid genome sequences and the shared 71 CDS sequences (Fig. 6 and Fig. S1). The trees strongly supported that Lilium was a monophyletic group. The support values of phylogenetic trees constructed by the whole plastid genome sequences were slightly improved when compared with the phylogenetic trees constructed by shared CDS sequences. The phylogenetic relationships of most species were relatively clear and had high support values. Two groups were identified among the 31 Lilium species. We considered traditional sections (Comber, 1949; Wang & Tang, 1980) and used Nomocharis or Lilium duchartreis (including L. duchartrei and L. lankongense) as a clade, respectively, for clarification. Group I was composed of 14 species and was divided into the Sinomartagon II + Martagon I, Leucolirion I, and Liriotypus clades. Group II was made up of the remaining 17 species, belonging to the Lilium duchartreis + Lophophorum, Sinomartagon I, Nomocharis, Martagon II, Leucolirion II, and Pseudolirium clades. The ML tree and the BI tree also indicated that Sect. Sinomartagon, Sect. Martagon, Sect. Leucolirion and Sect. Pseudolirium did not form monophyletic groups. Sect. Sinomartagon, Sect. Martagon and Sect. Leucolirion each formed two clades. L. lankongense, L. duchartrei, L. stewartianum, L. matangense, L. lophophorum, and L. nanum were clustered into the Lilium duchartreis + Lophophorum branch with strong support (BS = 100%, PP = 1.00). The phylogenetic trees (ML tree and BI tree) reconstructed by ITS sequencing also strongly suggested that Lilium was a monophyletic group (Fig. 7). The main structure of Lilium was divided into four groups, Group A included the Leucolirion I, Sinomartagon II, Sinomartagon I, and Martagon clades, Group B contained the the Lophophorum, Pseudolirium, and Leucolirion II clades, Nomocharis and Liriotypus clades made up Group C, and the Lilium duchartreis clade alone formed Group D. The phylogenetic tree, established by the 33 plastid genomes or the shared 71 CDS, was not completely consistent with the phylogenetic tree constructed by ITS. In the ITS phylogenetic tree, four species in Sect. Martagon and four species in Sect. Pseudolirium were respectively monophyletic clades. L. lankongense and L. duchartrei belonged to Lilium duchartreis clade, while L. stewartianum, L. matangense, L. lophophorum, and L. nanum belonged to another clade (Lophophorum clade), and the Lilium duchartreis clade and the Lophophorum clade were in two different groups. Compared with the phylogenetic trees formed by complete plastid genome sequences and its shared 71 CDS sequences, the ITS phylogenetic tree was more in line with the previous classification based on morphology.

Comparative analysis of Lilium plastid genomes

The plastid genomes of the six Lilium species were compared and analyzed. Similar to other Lilium species (Du et al., 2017; Kim & Kim, 2013; Liu et al., 2018a; Liu et al., 2018b), the plastid genomes of these six species all have LSC and SSC regions, and a pair of IRS (Table 1). Moreover, these six plastid genomes have the same protein-coding genes, tRNAs, and rRNAs (Table 1), and have the typical structure of land plant plastid genomes (Wicke et al., 2011). Their lengths varied from 152,307 to 152,611 bp (Table 1), indicating a slight variation. The LSC and IR boundaries and the J_SA region of the six plastid genomes were relatively stable with small differences, while the J_SB region was relatively obvious. The ndh F gene of the J_SB region showed signs of contraction and expansion. Typically, the IR (IR_A and IR_B) length differs among plant species. In angiosperms, the size variations in the plastid genome are attributed to the expansion and contraction of the IR region and single-copy (SC) border regions (Chumley et al., 2006) and play a vital role in its evolution (Wang et al., 2008). According to our results, the IR expansion and contraction may act as a primary mechanism for the variation in the length of the six plastid genomes. Similar results have been shown in most land plants (Gao et al., 2019; Hansen et al., 2007; Huang et al., 2014).

Microsatellites (SSRs) also known as short tandem repeats (STRs), are usually composed of one to six nucleotide repeat units (Chen et al., 2006; Powell, Machray & Provan, 1996), and are found in the majority of eukaryotic genomes (Tautz & Renz, 1984). SSRs are suitable molecular markers because of their high rate of polymorphism, and they are mostly co-dominant, and highly transferable (Scott et al., 2000). The SSRs in the plastid genome, which were viewed as complements to the nuclear microsatellite markers, were useful in research of phylogeography, population genetics, systematics, and breeding (Gaudeul et al., 2011; Hashimoto et al., 2004; Sánchez-Pérez et al., 2005; Xue, Wang & Zhou, 2012). The repeat units in plastid SSRs of these Lilium genomes were mostly comprised of adenine (A) or thymine (T) repeats, but were rarely guanine (G) or cytosine (C) repeats (Fig. 3A), which resulted from the AT-richness of the plastid genome (Zhang et al., 2016). These results are identical to other Lilium species (Du et al., 2017; Liu et al., 2018b). We identified a total of 340 SSR information sites. These SSR information sites can provide a reference for the subsequent development and application of SSRs molecular markers within Lilium.

Lilium-specific primers

Noncoding regions of the plastid genome play an important molecular role in systematic, phylogeographic, and population genetic studies of angiosperms (Perret et al., 2003; Xie et al., 2019a; Yang et al., 2016). Previous researchers, including Lai et al. (2016) tried to amplify the trnV-ndhC and petL-petE regions of the plastid genome but failed to amplify the corresponding fragments of all individuals. Jiang (2017) was unable to find suitable ptDNA fragments by screening 33 published universal primers for their phylogeographic study. No suitable ptDNA fragments from the genus Lilium have been found for use in population studies to date. In this study, 12 primers were designed according to the conserved regions of published Lilium plastid genomes to amplify a wide range of Lilium species. Two of the 12 pairs of primers (trnC-petN and rpl32-trnL) could amplify the products in all experimental individuals. The trnC-petN and rpl32-trnL regions proved to be effective in phylogeny and phylogeography in other genera (Liu, 2017; Lo Presti & Oberprieler, 2011; Yu et al., 2011) and the two non-coding regions of trnC-petN and rpl32-trnL may be used for intraspecies research of Lilium. The other 10 primer pairs (Table S2) successfully amplified some individuals, such as trnQ-rps16, trnF-ndhJ, and trnE-trnT, which had similar proportions of variant sites in amplified sequences and were applied to population-level studies in other genera (Chen et al., 2013; Dane & Liu, 2007; Qiu et al., 2009). We suggest that trnC-petN and rpl32-trnL with the two new primer pairs can be used in subsequent population studies (Table 3). Future studies should focus on discovering better primers for other potential regions to amplify their products. The two pairs of primers (trnC-petN and rpl32-trnL) were selected based on their plastid genomes and provide basic data for the genetic diversity and molecular breeding of this important wild flower of Lilium.

Phylogenetic analysis

We constructed a phylogenetic tree based on 3 data sets (the plastid genome sequences, the shared 71 CDS sequences, and ITS sequences). The tree supported Lilium as a monophyletic group with subgroups that include Sect. Sinomartagon and Sect. Leucolirion, which were both polyphyletic (Figs. 6, 7, and Fig. S1). Sect. Martagon and Sect. Pseudolirium were both monophyletic in the ITS tree, but none of them formed a monophyletic group in the plastid genome and the shared 71 CDS trees (Figs. 6, 7, and Fig. S1). Givnish et al. (2020) also sequenced the plastomes and 440 nuclear loci of 67 species of Lilium, and the results are similar to ours. Previous studies have also shown a conflict between the nuclear gene tree and the plastid gene tree in Lilium. This discordance may be due to hybridization in Lilium (Gao et al., 2013; Gao, Harris & He, 2015; Givnish et al., 2020; Gong et al., 2017; Huang et al., 2018). Due to the limited dispersion of seeds and the pollination by wind or pollinators, hybridization in Lilium is likely to occur among neighboring species. Hybridization may play an important role in Lilium flower morphology (Gao, Harris & He, 2015; Givnish et al., 2020). Furthermore, there are paralogous homologs in the ITS sequences of Liliaceae that may not reflect the true interspecific relationships in the phylogenetic tree (Day et al., 2014; Gao et al., 2012b; Huang et al., 2018). The 31 species of Lilium studied could not exclude the paralogous homology of the ITS sequences. Finally, L. lankongense, L. duchartrei, L. stewartianum, L. matangense, L. lophophorum, and L. nanum in the ITS tree were split into two different clades. The phylogenetic trees constructed by Huang et al. (2018) based on the ITS of 140 accessions in Liliaceae showed that the six species occur together in a clade. Aside from those results, the phylogenetic trees constructed with 440 nuclear genes and their spacers by Givnish et al. (2020) supported five species (except L. matangense), including L. lankongense, L. duchartrei, L. stewartianum, L. lophophorum, and L. nanum as a monophyletic group. The nuclear gene phylogenetic tree in this study only involved 33 ITS sequences (the number of samples was small), the information reflected by the ITS sequences was limited, and some clades (such as the Pseudolirium and Sinomartagon I clades) had low support values. The aforementioned three points are also one of the reasons for discordance between the ITS tree and the plastid gene tree.

In sections delimited by Comber (1949), which were based on morphology, L. lankongense, L. duchartrei, L. stewartianum, L. lophophorum, and L. nanum belong to Sect. Sinomartagon. L. lankongense, L. duchartrei belong to Sect. Sinomartagon 5a, while L. stewartianum, L. lophophorum, and L. nanum belong to Sect. Sinomartagon 5c (Comber, 1949). Sect. Sinomartagon Comber, which is the largest group of the genus Lilium, includes rather disparate species, and various researchers had different treatments for the division of its subsection (Comber, 1949; Du et al., 2014; Gao, Harris & He, 2015; Givnish et al., 2020; Nishikawa et al., 2001). Sect. Sinomartagon Comber was divided into the Sinomartagon, Lophophorum, and Lilium duchartreis clades based on the ITS tree (Fig. 7). The species in the Lophophorum clade contained L. stewartianum, L. matangense, L. lophophorum, and L. nanum; the Lilium duchartreis clade was composed of L. lankongense and L. duchartrei. The plastid gene trees (the plastid genome sequences, the shared 71CDS sequences) and the ITS tree both supported that Sect. Sinomartagon Comber was not monophyletic (Figs. 6, 7, and Fig. S1). In the plastid and ITS trees, L. lankongense was sister to L. duchartrei, L. stewartianum was evolutionarily close to L. matangense. Different phylogenetic relationships were shown among L. lankongense, L. duchartrei, L. stewartianum, L. matangense, L. lophophorum, and L. nanum. These six species were a monophyletic group in the plastid gene trees, and two clades were formed in the ITS tree. Previous studies have shown that L. lankongense and L. duchartrei was a monophyletic group, while L. stewartianum, L. matangense, L. lophophorum, and L. nanum were clustered into another clade in the ITS tree (Gao et al., 2013; Gao, Harris & He, 2015). Huang et al. (2018) constructed a phylogenetic tree based on the ITS sequence, three ptDNA fragments, and the combined (ITS + ptDNA) data sets. All trees supported that the six species were a monophyletic group (Sect. Sinomartagon I). Phylogenetic trees constructed with 440 nuclear genes and their spacers by Givnish et al. (2020) showed that the five species (L. lankongense, L. duchartrei, L. stewartianum, L. lophophorum, and L. nanum) aggregated into a single clade. Hence, L. lankongense, L. duchartrei, L. stewartianum, L. matangense, L. lophophorum, and L. nanum in the phylogenetic trees of ITS did not form a clade, which may be due to the insufficient number of ITS samples and the limited information contained in the ITS sequence. A close relationship has been shown among the six species (L. lankongense, L. duchartrei, L. stewartianum, L. matangense, L. lophophorum, and L. nanum) reflected in the plastid genome datum.

The revolute tepals of L. lankongense, L. duchartrei, L. stewartianum, L. matangense are different from those of the campanulate flowers of L. lophophorum and L. nanum (Fig. 1). However, the morphological characteristics of Lilium do not reflect the actual interspecies relationships. Convergent evolution, divergent evolution, retention and loss of ancestral characteristics, hybridization, and non-one-to-one correspondence between genes and morphological traits may affect plant’s morphology. The results of these five factors may be different morphologies of closely related species and similar morphologies of distantly related species, which confuses the judgment of interspecies relationships (Givnish et al., 2020). The habitats of L. lophophorum and L. nanum are mainly alpine grasslands and the adjacent niche at higher altitudes. These habitats often experience torrential downpours and have strong ultraviolet rays. The nodding and campanulate flowers can protect the reproductive structures against heavy rainfall and ultraviolet rays (Gao, Harris & He, 2015). The other four species found in the regions of lower altitudes have flowers with a larger opening, due to the revolute tepals, which may be more beneficial to pollination. The common floral characteristics of these six Lilium species are the inner tepals with projections on both surfaces of the nectaries, however the shape of the projections is not exactly the same. The projections of L. lankongense and L. duchartrei are papillose, L. stewartianum has cristate projections, and L. matangense, L. lophophorum, and L. nanum have fimbriate projections (Fig. S2). L. lankongense and L. duchartrei have habitats in lower elevations than the other four species and their inner and outer tepals have papillose projections. The other four species have habitats at higher elevations and only have projections in the inner tepals. L. stewartianum grows in open and rocky places on limestone mountains or along valleys and forest margins. Its habitat is distinct from those of the other three species, and the cristate projections are also different. L. matangense grows on grassy slopes, L. lophophorum and L. nanum grow mainly on bushy slopes and alpine grasslands. These three species have similar habitats and their projections are fimbriate. We speculate that the adaptation to different habitats may cause the dissimilarities of floral morphology among L. lankongense, L. duchartrei, L. stewartianum, L. matangense, L. lophophorum, and L. nanum.

Although the flower morphology of L. lophophorum and L. nanum was different from the other four species, they shared similar karyotypes that were measured using MAML (Altinordu et al., 2016). Gao, Zhou & He (2011) found that L. lophophorum, L. nanum, and L. matangense held a 3A type and relative asymmetry index and there were multiple, approximative parameters of the karyotype between L. stewartianum and L. lophophorum (Wan et al., 2011). The phylogenetic trees in this study demonstrated that the four species were more closely related (Figs. 6, 7, and Fig. S1). L. matangense, L. lophophorum, and L. nanum had karyotypes that resembled a clade formed by L. duchartrei and L. lankongense (Gao et al., 2012a). Therefore, it is reasonable to conclude that these six species constituted a clade. L. duchartrei and L. lankongense, which are closely related (Jefferson-Brown & Howland, 1995; MacRae, 1998), were thought to be one species (Wang & Tang, 1980), but were later separated due to the different colors of the tepal (L. lankongense: pink tepal, L. duchartrei: white tepal), the obvious geography variations (within the Daxue Mountains), and the significant genetic differentiation of the 22 populations (Gao et al., 2013; Liang & Tamura, 2000; Shen, Zhou & He, 2014). There are obvious micromolecular differences in the leaf dermis such as the inner margin and the waxy ornamentation of the outer stomatal rim, and the characteristics of the cuticular membrane under the electron microscope (L. lankongense: the pattern of anticlinal wall of leaf epidermal cells are all straight and the guard cells have no T pieces; L. duchartrei: the pattern of anticlinal wall of leaf epidermal cells are wavy and the guard cells have T pieces at the two poles) (Zhou, 2008). Our research supports the classification of L. lankongense as an independent species and agrees with the classification of Flora of China (Liang & Tamura, 2000). The two species have been clustered in the same clade based on the internal transcribed spacer region or plastid genes (Du et al., 2014; Gao et al., 2012a; Givnish et al., 2020; Lee et al., 2011; Nishikawa et al., 2001) (Figs. 6, 7, and Fig. S1). L. duchartrei also showed a sister relationship with L. lankongense in the phylogenetic trees of ITS and plastid genomes, which is in agreement with results from previous studies (Figs. 6, 7, and Fig. S1). We reconstructed the phylogenetic relationship of Lilium with 33 plastid whole genomes. However, we only conducted a preliminary analysis of the relationship of some species in Lilium due to the limited plastid genome datum. More plastid genomes should be obtained to better clarify the classification at the section level of Lilium and the relationship among species.