Species of Tulipa (Liliaceae) are of great horticultural importance and are distributed across Europe, North Africa, and Asia. The Tien Shan Mountain is one of the primary diversity centres of Tulipa, but the molecular studies of Tulipa species from this location are lacking. In our study, we assembled four Tulipa plastid genomes from the Tien Shan Mountains, T. altaica, T. iliensis, T. patens, and T. thianschanica, combined with the plastid genome of T. sylvestris to compare against other Liliaceae plastid genomes. We focussed on the species diversity and evolution of their plastid genomes. The five Tulipa plastid genomes proved highly similar in overall size (151,691–152,088 bp), structure, gene order, and content. With comparative analysis, we chose 7 mononucleotide SSRs from the Tulipa species that could be used in further population studies. Phylogenetic analyses based on 24 plastid genomes robustly supported the monophyly of Tulipa and the sister relationship between Tulipa and Amana, Erythronium. T. iliensis, T. thianschanica, and T. altaica were clustered together, and T. patens was clustered with T. sylvestris, with our results clearly demonstrating the relationships between these five Tulipa species. Our results provide a more comprehensive understanding of the phylogenomics and comparative genomics of Tulipa.

1. Introduction

Plastid, a unique semiautonomous organelle of green plants, serves as a protagonist in photosynthesis and carbon fixation and provides essential energy for plants [1, 2]. Plastid DNA has been widely used in evolutionary biology analyses because of its uniparental inheritance, relatively stable genome structure, and gene content [37]. A plastid genome has a quadripartite circular structure consisting of two copies of inverted repeat (IR) regions, a large single copy (LSC) region, and a small copy region (SCR) in most seed plants [810]. Comparative genomics of whole plastid genomes has been used to generate genetic markers for molecular identification [11, 12]. Plastid genomic data provide new and robust insights into phylogenetic analysis and genetic variation detection of plants [13].

Tulips (Tulipa L.) are famous ornamental and cut flowers due to their beautiful and colorful corolla [1416]. Tulipa is a member of Liliaceae sensu APG IV [17], subfamily Lilioideae, tribe Tulipeae Kostel [18, 19]; the tribe Tulipeae including Gagea Salisbury, Amana Honda, Erythronium L., and Tulipa [20] is sister to the tribe Lilieae [21]. The close relationship between Gagea and Tulipa, Amana, and Erythronium is generally accepted; however, the phylogenetic relationships among Tulipa, Amana, and Erythronium have remained controversial due to inconsistencies across studies and lack of strong support [22, 23]. Amana used to be treated as a group in Tulipa [2426] but is now generally accepted as a separate genus [2729]. Several studies have supported a sister relationship between Tulipa and Amana [30, 31], whereas others clustered Tulipa and Erythronium together [22, 29]. A close relationship between Amana and Erythronium has been suggested; however, the Tulipa data used in these studies was limited [3133].

Tulipa includes more than 100 species [34, 35] and are distributed across the temperate regions of Europe, North Africa, and Asia. The Tien Shan, Pamir-Alay, and Caucasus Mountains are considered the primary diversity centres of Tulipa [36, 37]. Thirteen species of Tulipa are recorded and described in China, and 11 species are distributed in the Xinjiang Uygur Autonomous Region [37], which is home to the Tien Shan Mountain. Several studies have utilised a few plastid regions and internal transcribed spacer makers to analyse the phylogeny and evolution of Tulipa in the Middle East [38, 39] and Europe [29]. However, the phylogeny and evolution of Tulipa at Tien Shan Mountain are limited. Two Tulipa plastid genomes have been published [40, 41], but genomics analysis for Tulipa is lacking. Therefore, comprehensive studies including the specimens from Tien Shan Mountain and accurate analysis of the plastid genome are required to enable further Tulipa molecular studies.

We collected four Tulipa species from Tien Shan Mountain and reported the complete plastid genome sequences of these four Tulipa species to address gaps in Tulipa phylogenetic research. Combining previously reported plastid genomes of T. sylvestris and other Liliaceae species, we performed comparative genomics and phylogenetic analyses. We aimed to resolve the phylogenetic relationships between five Tulipa species and other genera and characterize and compare the plastid genomes of Tulipa species to detect the genetic variation. Our conclusion will contribute to an understanding of Tulipa plastid phylogenomics and provide genetic resources for tulip research.

2. Materials and Methods

2.1. Plant Materials, DNA Extraction, and Sequencing

Fresh leaves of four Tulipa species, T. altaica, T. iliensis, T. patens, and T. thianschanica, were collected from Yuming county (Xinjiang Uygur Autonomous Region, China) and dried with silica gel, then stored at -80°C. Total genomic DNA was extracted from leaf material with a modified CTAB method [42] and then sequenced on an Illumina Novaseq2500 sequencer (Illumina, San Diego, CA, USA) by Biomarker Technologies, Inc. (Beijing, China).

2.2. Plastid Genome Assembly, Annotation, and Analysis

The plastid genomes were assembled using raw data by NOVOPlasty 2.7.2 [43], and the plastid genome of T. sylvestris (MT261172) was selected for seed input and the reference sequence. Genome annotation and IR region search were processed by PGA [44]. Geneious R11 [45] was used on manual modifications to accurately confirm the start and stop codons and the exon-intron boundaries of genes based on comparison with other Liliaceae plastid genomes. The circular plastid genome map of Tulipa was drawn by the OGDRAW1 program [46]. The total GC content and GC content of each region (IR, LSC, SSC) were analysed by the program Geneious R11.

2.3. Contraction and Expansion of IRs and SSRs

Four plastid genomes of the tribe Tulipeae, Amana edulis (NC034707), Erythronium japonicum (MT261155), Erythronium sibiricum (NC035681), and Gagea triflora (MT261157), were downloaded from the GenBank for comparative analysis with five Tulipa species.

The IR/SC borders with full annotations were compared between the five Tulipa species and with the other four tribe Tulipeae species using the program IRscope (https://irscope.shinyapps.io/irapp/) [47]. Simple sequence repeats (SSRs) were detected using Perl script MISA [48] with the following minimum number (threshold) settings: 10, 5, 4, 3, 3, and 3 repeat units for mono-, di-, tri-, tetra-, penta-, and hexanucleotide SSRs, respectively.

2.4. Phylogenetic Analyses

To reconstruct phylogenetic relationships between the five Tulipa species and other Liliaceae species, a total of 20 plastid genome sequences were downloaded from GenBank, and Smilax china (Smilacaceae, HM536959) was selected as an outgroup. The alignment of 24 plastid genome sequences was generated by MAFFT v7.402 [49] with the default parameters set. The best-fit model selected by ModelFinder was GTR+G. A maximum likelihood (ML) tree was constructed using RAxML 8.0 [50] with 1000 bootstrap replicates.

2.5. Codon Usage Analysis and SNP Analyses

All 84 protein-coding sequences extracted from nine plastid genomes were used to analyse codon usage, which was undertaken with the CodonW v1.4.2 program (J. Peden, http://codonw.sourceforge.net). The plastid genome sequences of nine tribe Tulipeae species were used for SNP analyses. The alignment of all plastid genome sequences was generated by MAFFT v7.402. SNP analysis was conducted to determine the nucleotide diversity of the plastid genomes using DnaSP v5, with the following parameters: 200 bp of step size and 600 bp window length [51]. Results of the SNP analysis were illustrated using TBtools v1.087 software [52].

3. Results and Discussion

3.1. The Plastid Genomes of Five Tulipa Species

The total plastid genome sizes of the five Tulipa species ranged from 151691 bp (T. altaica) to 152088 bp (T. patens). All five plastid genomes showed the typical quadripartite structure (Figure 1) and, like other angiosperms, consisted of a pair of IR regions. The G+C content of the five species in whole genomes (36.6-36.7) and LSC (34.5-34.6), SSC (30.0-30.2) was nearly identical but in the IR regions was higher (42.0%). Details of genome features are given in Table 1. The annotated genome sequences of T. altaica, T. iliensis, T. patens, and T. thianschanica were deposited in the GenBank under the accession numbers MW077741, MW077740, MW077739, and MW077738, respectively (Table 1).

The plastid genomes of five Tulipa species contained 134 genes. Of these 134 genes, 112 genes were nonredundant including 78 protein-coding genes, 4 ribosomal RNA (rRNA) genes, and 30 transfer RNA (tRNA) genes, and four genes were pseudogenes (Table 2). The 134 genes had 18 duplicated genes located in the IR region, including six coding genes (ndhB, rpl2, rpl23, rps7, rps12, and ycf2), four rRNA genes (rrn4.5, rrn5, rrn16, and rrn23), and eight tRNA genes (trnA-UGC, trnH-GUG, trnI-CAU, trnI-GAU, trnL-CAA, trnN-GUU, trnR-ACG, and trnV-GAC).

Four pseudogenes (ycf1, rps19, and two ycf68) were found in the five plastid genomes (Table 2). The rps19 and ycf1 genes were located in the boundary area of the IR regions, and their protein-coding ability was lost due to partial gene duplication [6, 10, 40, 53]. Whether ycf68 and ycf15 genes lost abilities or occur as pseudogenes has already been discussed in several studies [40, 53, 54]. In this study, the ycf15 gene was not annotated due to its short length. The infA gene, which codes for translation initiation factor 1, was lost in all five Tulipa plastid genomes because of a missing base. The deletion of the infA gene also occurred in Amana and Erythronium [40], which were sister relationships with Tulipa, and many other seed plants, such as Smilax (Smilacaceae) [40] and Alstroemeria (Alstroemeriaceae) [33].

3.2. Inverted Repeats Contraction, Expansion, and SSR Analysis

The IR/SC boundary regions of the Tulipa plastid genomes were compared to the closely related plastid genomes, Amana, Erythronium, and Gagea. Typically, the lengths of IR regions are different among various plant species [53], while the lengths of Tulipa plastid genome IR regions were similar (26307 bp-26341 bp) but larger than Amana (25633 bp), Erythronium (25765 bp and 26001 bp), and Gagea (25521 bp) plastid genome IR regions (Figure 2). Furthermore, the expansion and contraction at the IR regions were the primary cause of size variation in plastid genomes and played an important role in the evolution of the genome [5557]. After comparing the location and adjacent genes of IR regions between nine plastid genomes, we found that the gene number and order were conserved, but some distinct differences existed at the boundaries (Figure 2). The boundary of the LSC and IRb regions was located at rps19 genes in eight plastid genomes and positioned at the noncoding region between rps19 and rpl2 genes in the Amana edulis plastid genome. The ndhF and ycf1 genes traversed the regions of IRb/SSC and SSC/IRa, with 38/40 bp of the ndhF gene and 1587 bp/1589 bp of the ycf1 gene located at the IR region in the Tulipa plastid genome. In general, the length and structure of IR regions were similar in Tulipa genomes but showed obvious differences with other tribe Tulipeae genomes.

Simple sequence repeats (SSRs) in the plastid genomes are suitable molecular makers and have been widely used in evolutionary and ecological studies due to their high variation [5864]. Given that SSRs have high polymorphism at the species level and commonly show intraspecific variation, SSRs were used as important molecular markers to reconstruct phylogenetic relationships [6567]. The results of SSR analysis of Tulipa and its close relatives are shown in Figure 3 and Table S1-S2. There were six categories of SSRs (mono-, di-, tri-, tetra-, penta-, and hexanucleotide repeats) found in the plastid genome of nine species, and the mononucleotide repeats were the most frequent (Figure 3(a), Table S1). The highest percentages of SSRs located in the LSC region were the maximum (69.44%-79.17%), and 4%-10% SSRs were distributed in IR regions (Figure 3(b), Table S2). A total of 41 types were detected in nine plastid genomes (Figure 3(c)), where bases T and A were the dominant elements. In this study, we manually chose 7 mononucleotide SSRs located in trnK-rps16, psbK-I, accD-psaI, and psaJ-rpl33 regions and atpF, rpoC1, and petB genes as effective polymorphic SSRs between Tulipa species based on three critical criteria outlined in previous research [40]. These 7 mononucleotide SSRs could be used in the further population studies of Tulipa (Table S3).

3.3. Codon Usage Bias and SNP Analyses

The results of codon usage frequency and relative synonymous codon usage (RSCU) from 84 protein-coding sequences from nine tribe Tulipeae species are presented in Figure 4 and Table S4. The 84 protein-coding sequences were similar across the nine species. The total codon number from the nine species plastid genomes ranged between 21170 (Erythronium sibiricum) and 21284 (T. thianschanica) (Table S4). Among all amino acids, leucine and cysteine were the most and the least frequent, on average, 2165 (10.19%) and 333 (1.57%), respectively (Figure 4 and Table S4). The third codon position occupied by the A or T base was the most common in all nine tribe Tulipeae species, which is also found in other plastid genomes in seed plants [6871]. Codon usage bias was related to gene expression and differed between species [72, 73]. Our results of codon usage bias will be important for understanding the molecular evolution mechanisms of Tulipa and its relatives.

The single-nucleotide polymorphism (SNP) analyses of the alignment for the nine tribe Tulipeae plastid genomes showed that the IR regions were more conserved than the LSC and SSC regions, where the SNP number was low (Figure 5), and was in agreement with previous reports for the angiosperm plastid genome [74, 75]. The noncoding regions were more variable than the coding-protein regions. In the sequence alignment of the nine tribe Tulipeae plastid genomes, six noncoding regions, rps16-trnQ, trnE-trnT, accD-psaI, rpl32-trnL, rps15-ycf1, and rps4-trnT, and ycf1 and ndhA genes were highly variable.

3.4. Phylogenetic Analysis

The phylogenetic relationship of Tulipa and other Liliaceae species was reconstructed based on 24 plastid genomes, representing 14 genera (Figure 6). The maximum likelihood (ML) tree strongly supported Tulipa as a monophyletic genus that was sister to Amana and Erythronium (1/100%), which was the same as previous research [29, 40]. Tulipa, Amana, Erythronium, and Gagea formed a monophyletic clade (1/100%). The ML analytical result of Tulipa and other Liliaceae species based on plastid genomes was in accordance with APG IV [17]. Previous studies could not resolve the phylogenetic relationships of Tulipa, Amana, and Erythronium because of insufficient data. In our ML tree, inferred from 23 Liliaceae species’ plastid genomes, Amana and Erythronium clustered together with strong support (1/100%), and Tulipa was a sister to the Amana/Erythronium clade. The phylogenomics reconstructed distinct relationships among Tulipa and other genera.

In Tulipa, T. iliensis and T. thianschanica were very close and were sisters to T. altaica. We also found that T. patens was sister to T. sylvestris. In previous research [26, 38], Tulipa was divided into four subgenera: Clusianae, Tulipa, Eriostemones, and Orithyia. T. iliensis and T. altaica belong to subgenus Tulipa [37], yet T. iliensis, T. thianschanica, and T. altaica have similar morphological characteristics, and our results confirmed this grouping. T. patens and T. sylvestris were within the subgenus Eriostemones, and T. patens was once treated as a varietas of T. sylvestris [37]. Our phylogenetic tree based on plastid genomes demonstrated clear relationships between the five Tulipa, which was in accordance with a previous classification [37], and found that T. thianschanica should belong to subgenus Tulipa. Our results provide a better understanding of the evolution and molecular biology of Tulipa.

4. Conclusions

In this study, the plastid genome sequences of fourTulipa species were reported. We described the comparative characteristics of nine tribe Tulipeae plastid genomes and the phylogenetic relationships of 23 Liliaceae plastid genomes. We found that Tulipa plastid genomes were highly similar in overall size, structure, IR/SC boundary, SSRs, and codon usage bias. The phylogenetic tree identified a clear sister relationship between Tulipa and Amana, Erythronium and clear relationships between the five Tulipa species. Our study supplements the molecular data of Tulipa and provides a better understanding of Tulipa plastid genome evolution.

Data Availability

The assembled plastid genome sequences of the four species were submitted to NCBI with the accession numbers MW077741 (T. altaica), MW077740 (T. iliensis), MW077739 (T. patens), and MW077738 (T. thianschanica). Users can download the data as a reference for research purposes only.

Conflicts of Interest

There was no conflict of interest.


We acknowledge Hai-Ying Liu, Yi-Qi Deng, Sheng-Bin Jia, and Qiu-Pin Jiang for their help in material collection. This work was supported by the National Natural Science Foundation of China (Grant Nos. 31872647 and 32070221); the National Specimen Information Infrastructure, Educational Specimen Sub-Platform (Grant No. 2005DKA21403-JK); and the Fourth National Survey of Traditional Chinese Medicine Resources (Grant No. 2019PC002).

Supplementary Materials

Table S1: number of different SSR categories detected in nine species. Table S2: the frequency of identified SSRs in LSC, IR, and SSC of nine species. Table S3: seven polymorphic SSRs between Tulipa species. Table S4: the codon numbers of amino acids in nine plastid genomes. (Supplementary Materials)