Skip to main content

Morphological and ecological divergence of Lilium and Nomocharis within the Hengduan Mountains and Qinghai-Tibetan Plateau may result from habitat specialization and hybridization



Several previous studies have shown that some morphologically distinctive, small genera of vascular plants that are endemic to the Qinghai-Tibetan Plateau and adjacent Hengduan Mountains appear to have unexpected and complex phylogenetic relationships with their putative sisters, which are typically more widespread and more species rich. In particular, the endemic genera may form one or more poorly resolved paraphyletic clades within the sister group despite distinctive morphology. Plausible explanations for this evolutionary and biogeographic pattern include extreme habitat specialization and hybridization. One genus consistent with this pattern is Nomocharis Franchet. Nomocharis comprises 7–15 species bearing showy-flowers that are endemic to the H-D Mountains. Nomocharis has long been treated as sister to Lilium L., which is comprised of more than 120 species distributed throughout the temperate Northern Hemisphere. Although Nomocharis appears morphologically distinctive, recent molecular studies have shown that it is nested within Lilium, from which is exhibits very little sequence divergence. In this study, we have used a dated molecular phylogenetic framework to gain insight into the timing of morphological and ecological divergence in Lilium-Nomocharis and to preliminarily explore possible hybridization events. We accomplished our objectives using dated phylogenies reconstructed from nuclear internal transcribed spacers (ITS) and six chloroplast markers.


Our phylogenetic reconstruction revealed several Lilium species nested within a clade of Nomocharis, which evolved ca. 12 million years ago and is itself nested within the rest of Lilium. Flat/open and horizon oriented flowers are ancestral in Nomocharis. Species of Lilium nested within Nomocharis diverged from Nomocharis ca. 6.5 million years ago. These Lilium evolved recurved and campanifolium flowers as well as the nodding habit by at least 3.5 million years ago. Nomocharis and the nested Lilium species had relatively low elevation ancestors (<1000 m) and underwent diversification into new, higher elevational habitats 3.5 and 5.5 million years ago, respectively. Our phylogeny reveals signatures of hybridization including incongruence between the plastid and nuclear gene trees, geographic clustering of the maternal (i.e., plastid) lineages, and divergence ages of the nuclear gene trees consistent with speciation and secondary contact, respectively.


The timing of speciation and ecological and morphological evolutionary events in Nomocharis are temporally consistent with uplift in the Qinghai-Tibetan Plateau and of the Hengduan Mountains 7 and 3–4 million years ago, respectively. Thus, we speculate that the mountain building may have provided new habitats that led to specialization of morphological and ecological features in Nomocharis and the nested Lilium along ecological gradients. Additionally, we suspect that the mountain building may have led to secondary contact events that enabled hybridization in Lilium-Nomocharis. Both the habitat specialization and hybridization have probably played a role in generating the striking morphological differences between Lilium and Nomocharis.


The Hengduan Mountains (H-D Mountains) are located in southwestern China east of the Qinghai-Tibetan Plateau (QTP) and represent one of the world’s most biodiverse regions [1]. Many endemic vascular plant species of the H-D Mountains exhibit high levels of morphological and ecological divergence from their closest, more widespread allies. Thus, the endemics are often treated within their own genera. However, molecular phylogenetic studies have revealed that the some of these endemic genera are nested within the widespread ones. Examples include representatives of Asteraceae (Sinacalia), Brassicaceae (Solms-laubachia), Liliaceae (Lloydia), Primulaceae (Pomatosace), Genetianaceae (Lomatogoniopsis), and Amaryllidaceae (Milula) (see more detail information in Table 1, [28]). The contrasting morphological diversity and nested phylogenetic status of genera in the H-D Mountains may result from extreme habitat specialization and/or hybridization events. The H-D mountains provide many unique habitats due to their topographic complexity [9], while repeated phases of uplift of the mountain range may have enabled opportunities for hybridization [10, 11] via secondary contact. Continued research is needed to better understand the mechanisms driving morphological diversity of vascular plants within the H-D Mountains.

Table 1 Morphologically distinctive plant species that are endemic to the QTP but phylogenetically indistinct (i.e., nested within) from allies

The Lilium-Nomocharis complex represents an exceptional study system for morphological diversification and hybridization in the H-D Mountains. Nomocharis Franchet. is endemic to the H-D Mountains and adjacent QTP. Nomocharis appeared somewhat similar to Lilium when the former was first described in 1889 [12, 13] but was erected as a new genus because of its highly distinctive open-plate flowers and dark-colored tepal bases with special structures (Figs. 1 and 2) [1215]. Currently, there are eight recognized species of Nomocharis, of which seven are circumscribed in two traditional sections [14, 15], and one is a recently described hybrid species, N. gongshanensis Y. D. Gao & X. J. He [16]. Recent molecular phylogenetic studies show strong support for Nomocharis nested within Lilium [16, 17]. In contrast to Nomocharis, Lilium comprises approximately 120 species and is widespread throughout the Northern Hemisphere, including areas within the QTP and H-D mountains [1820].

Fig. 1
figure 1

Pictures of Nomocharis aperta in western Yunnan: (a-c), population from Zhongdian, Yunnan showed spot variation; (c-e), population of Fugong, Yunnan showed variations in tepal color; (f-h), habits of N. aperta under different habitats; (i-j), anatomical pictures showed two types of N. aperta from Zhongdian and Fugong, as well as a comparison of outer and inner tepals

Fig. 2
figure 2

Pictures from western China showing Nomocharis: (a-c), N. basilissa; (d-f), N. farreri; (g-i), N. gongshanensis; (j-l), N. meleagrina

The goals of our present study are to use a molecular phylogeny as a framework to 1) determine whether the timing of morphological and ecological evolutionary events in Nomocharis are consistent with phases of uplift in the H-D Mountains and QTP, and 2) detect additional hybridization events with the Lilium-Nomocharis species of the H-D Mountains and QTP.


Phylogenetic analyses

A large ITS dataset confirmed the phylogentic position of Nomocharis within Lilium and showed no major differences compared with previous studies (e.g., [16, 17, 21]). Our extensive sampling of Nomocharis enabled us to resolve three sublclades within the genus: Eunomocharis, Ecristata, and the Non-Nomocharis lilies (Lilium species, N-N, hereafter). The Eunomocharis and Ecristata subclades are congruent with traditional classifications based on morphology [13]. The N-N lilies are morphologically divergent from Nomocharis and have characteristics more like other Lilium (Fig. 3). Nomocharis and the N-N lilies are sister to a clade comprised of Lilium sect. Liriotypus (i.e., European lilies) and that these two clades are sister to the rest of Lilium (Additional file 1: Figure S1).

Fig. 3
figure 3

Maximum credibility tree showing monophyletic clade of Nomocharis and its relatives reconstructed using Bayesian analysis of ITS data and Lilium species from around the world. The position of this clade is indicated on the tree (for details see Additional file 1: Figure S1). Support values shown on braches; Bayesian posterior probabilities (PP) on left and parsimony bootstrap (BS) on right. Clade names based on Balfour [12]

Major clades of the plastid consensus trees were the same in the Bayesian and MP reconstructions, so we present only the Bayesian consensus (Fig. 4). The plastid data resolved two large clusters consisting of seven major clades (Fig. 4). Cluster I (PP = 1.00, BS = 99 %) comprised two major clades of species of Lilium that are primarily distributed throughout the Sino-Japanese Forest subkingdom [22]. Cluster II (PP = 1.0, BS = 90 %) contained Nomocharis and species of Lilium that occur within the H-D Mountains and adjacent Himalayas.

Fig. 4
figure 4

Maximum credibility tree resulting from a Bayesian analysis of combined plastid DNA. Clade names based on Comber [23] and Liang [19]. Distributional areas of clades indicated by color. Support values shown on braches; Bayesian posterior probabilities (PP) on left and parsimony bootstrap (BS) on right. Lineages identified in network (Fig. 5) were also marked for references. The Sinomartagon I clade is highlighted for its conflicting position compared to the ITS result in Additional file 1: Figure S1

Within the plastid phylogeny, Nomocharis formed a poorly resolved grade with species of the Sinomartagon and Leucolirion clades. Most of the species of Sinomartagon that associated with Nomocharis and the N-N lilies occur in the Sinomartagon I clade in the ITS topology and represent all Sinomartagon species that inhabit the H-D Mountains and QTP [23, 24]. Despite poor resolution of Nomocharis within the plastid phylogeny, the genus roughly comprised its traditionally recognized sections, sects. Ecristata and Eunomocharis. A clade of Ecristata included N. aperta accessions and N. saluenensis, which have been have been historically treated in the section. The Ecristata clade also contained clones N. gongshanensis, which is the hybrid species, L. nepalense, and N. meleagrina, which is morphologically similar to species of Eunomocharis by having whorled leaves and has traditionally been circumscribed in that section. A grade of sect. Eunomocharis also included one accession of N. aperta (Franchet) E.H. and Lilium yapingense, an N-N lily species.

Overall, Nomocharis and the N-N lilies exhibited poorly resolved relationships within cluster II of the plastid phylogeny and did not form a monophyletic group.

Statistical parsimony network

Our parsimony network was complex but relatively well resolved (Fig. 5). Interior haplotypes and their descendants appear to represent eight lineages, most of which are present in the dichotomously branching plastid phylogeny (Fig. 4). The network supported the plastid tree topology in showing that geographically proximal species have more closely related haplotypes irrespective of morphological similarities or classification in traditional subgenera. Notably, the plastid tree and network also agreed in the placement of Nomocharis. In the network, Nomocharis was divided into two lineages, II and IV, and separated by Lineage III in which N-N lilies were included (Fig. 5). Haplotypes of the Nomocharis and the N-N lilies of lineages III and IV exhibit a shared history with Sinomartagon and Leucolirion species of lineage VI and VII as well as with species of a Lilium clade (lineage VIII, compare to Fig. 4).

Fig. 5
figure 5

Parsimony network conducted by TCS [58] using combined plastid DNA matrix. Sixty-six haplotypes were identified and clustered in eight lineages with different colors. Circle sizes correspond to the number of taxa possessing the haplotype. Species names are abbreviated by the generic first letter and two or three letters of the species epithet (Table 2). Inferred haplotypes (not present in the data set) are depicted as black lines, and unnamed dots indicated the missing interior haplotypes. The Sinomartagon I clade was highlighted for its conflict position compared to the ITS result in Additional file 1: Figure S1

Divergence time estimate and biogeography inferences

We performed divergence time dating using two secondary calibration points applied to our ITS plastid dataset. According to dating using the plastid dataset, and we inferred that the last shared ancestor of the Lilium-Nomocharis occurred around 13.19 Mya and Nomocharis evolved 6.5 Mya (Fig. 6). The ITS dataset recovered a slightly older age of approximately 14 Mya for the last shared ancestor of Lilium-Nomocharis and ca. 12 Mya for the evolution of Nomocharis (Fig. 7). Overall, the ITS dates for major diversification events are older than the plastid dates (Figs. 6 and 7).

Fig. 6
figure 6

Ultrametric chronograms showing divergence time dating and biogeographic results based on the combined plastid DNA phylogeny. Scale bar at bottom indicating branch length of 2 Mya. Mean divergence age given on nodes. Bars on nodes indicate the 95 % HPD for divergence ages. Pie charts show probabilities of ancestral area reconstructions, colors of pie slices defined in legend. The bottom chart summarized the biogeographic event through time. The Sinomartagon I clade was highlighted for its conflict position compared to the ITS result in Additional file 1: Figure S1

Fig. 7
figure 7

The ancestral state reconstructions of leaf, flower, and ecological characters. Pie charts show probabilities of ancestral area reconstructions, colors of pie slices defined in legend. Reconstructions of a, leaf arrangement, b, stigma:stamen ratio, c, corolla shape, d, corolla orientation with respect to the ground, and e, elevational range

The results from Bayesian Binary Method (BBM) of biogeographic analysis show that the last shared ancestor of Lilium-Nomocharis arose in the H-D Mountain region (B: 78.4 %; Fig. 6), while the results from the DEC method in Lagrange support a broader ancestral area within the H-D Mountains and the adjacent Sino-Japanese Floristic Subkingdom (SJFS; BC: 21.4 %; Fig. 6). The results obtained from BBM and DEC may not be incongruent because no significant geographic boundary separated the H-D Mountains and the SJFS areas until at least late Miocene (~7 Mya), which is the earliest date postulated for the H-D Mountain uplift [25, 26]. Lilium-Nomocharis began intensive diversification in the late Miocene (ca. 11–5 Mya, Fig. 6 or ca. 13–6 Mya, Fig. 7). The three Nomocharis lineages, Eunomocharis, Ecristata, and the N-N lilies, originated approximately between ca. 8 Mya (ITS, Fig. 7) and 6 Mya (plastid, Fig. 6) and underwent diversification during the late Pliocene beginning ca. 7–4 Mya (Figs. 6 and 7 respectively).

Ancestral state reconstruction (ASR)

We performed our ancestral state reconstructions using a reduced ITS dataset and they showed that floral characters were more phylogenetically dependent than vegetative ones (Fig. 7). Leaf arrangement patterns showed the greatest lability within clades (Fig. 7a). Overall whorled leaves arose at least four times in Lilium, including two shifts to whorled leaved within Nomocharis and the N-N lilies occurring approximately 4 Mya and 2.5 Mya, respectively (Fig. 7). Our results show that nodding flowers with recurved tepals and roughly equal stigma and stamen lengths are most likely the ancestral condition for Lilium (Fig. 7b, c, d). Ancestors of Nomocharis had longer stigmas than stamen, and this feature also was a synapomorphy within the sympatric Sinomartagon I clade (Fig. 7b). However, one species of Nomocharis, N. saluenensis, experienced a reversion to the roughly equal condition about 1 Mya (Fig. 7b). There appeared to be a correlation between floral orientation and corolla shape; namely that species with campanifolium and recurved petals have nodding flowers, and species with flat open and funnel/trumpet shaped flowers are horizon in orientation (Fig. 7c, d). This seems to be true among modern species and reconstructed ancestors. Recurved and campanifolium petals and the nodding habit evolved in the last shared ancestor of the N-N lilies around 7.5 Mya, and distinguish them from Nomocharis, which retained flat/open flowers and horizon orientation (Fig. 7c, d). The elevation reconstruction indicate that the ancestors of Nomocharis and the N-N lilies occurred at low (<1000 m) elevations and that radiations into different elevations habitats occurred around 5.5 Mya in the N-N lilies and around 3.5 Mya in the Ecristata clade of Nomocharis (i.e., including N. aperta accessions and N. saluenensis; Fig. 7e).


Morphological divergence and habitat specialization

Traditionally, classification of Lilium has focused primarily on floral morphology, especially orientation of the flowers with respect to the ground and corolla shape. Thus, nodding flowers and campaniform corollas have been used to support a close relationship between the N-N lilies, which include L. nepalense, L. souliei, L. paradoxum, L. saccatum and L. yapingense (Additional file 3: Figure S3, Additional file 4: Figure S4), and sect. Lophophorum (e.g., Lilium nanum, Additional file 4: Figure S4h, k, and L. lophophorum, Additional file 3: Figure S3d, e, f, of sect. Lophophorum), which shares the same floral features [23]. However, our ITS phylogeny is in contrast to traditional classification of the N-N lilies with sect. Lophophorum and shows that the N-N species are nested within Nomocharis, which is otherwise monophyletic (Figs. 3, Additional file 1: Figure S1). The N-N lilies share few apparent morphological traits in common with Nomocharis and, in particular, lack the unique floral characters that have classically been used to delimit Nomocharis from Lilium.

N-N lilies and traditional Nomocharis may exhibit morphological dissimilarities despite their close evolutionary relationships due to habitat specialization. The N-N lilies may have expanded their habitats into diverse elevations around 5.5 Mya that became available after the last QTP orogeny, which occurred ca. 7 Mya [27, 28] (Fig. 5e). Similarly, uplift of the H-D Mountains probably provided new habitat for an ancestor of the Ecristata clade of Nomocharis. Within the QTP, the N-N lilies tend to occupy higher elevations than the Nomochrais species of the H-D Mountains. Differential adaptations to elevation may explain the strikingly different floral morphology of Nomocharis and the N-N species [29]. In particular, the N-N lilies live almost exclusively in alpine meadows. Thus, N-N lilies are exposed to torrential downpours in alpine meadows compared to traditional Nomocharis species, which grow in the herbaceous layer beneath bamboo canopies (Additional file 5: Figure S5b, h) [19, 20]. The N-N lilies may have evolved nodding flowers ca. 7.5 Mya during QTP uplift and campaniform corollas as advantageous protections for their delicate reproductive structures against harsh precipitation conditions [30, 31]. Although the nodding, campaniform flowers probably provide protection from rainfall for the N-N lilies, they may also have reduced pollen transfer efficiency as an evolutionary trade-off [13, 14]. In contrast, Nomocharis species are probably not limited by the need for protection from heavy rainfall, and may experience higher pollen transfer efficiency by virtue of their horizontally arranged, plate-shaped flowers [13, 14].

The profound effects of habitat specialization within the H-D Mountains and QTP regions on morphology is supported by evidence of convergent evolution among sympatric, distantly related Lilium-Nomocharis species. In particular, Nomocharis and N-N lilies share some morphological traits in common with species of the Lophophorum clade, despite their differences and with which they are sympatric in alpine areas of the QTP. Shared traits especially include inner perianth-segments that have crested or fringed glandular bases (e.g., L. nanum and L. lophophorum Additional file 6: Figure S6) and that are sometimes anthocyanin rich (e.g., L. henrici Additional file 6: Figure S6). These shared morphological traits appear to represent convergent evolution. Morphological convergence within QTP alpine plant genera has been noted in other plant genera including in Androsace (Primulaceae) [5], Pseudoeriocoryne (Asteraceae: Cardueae) [32], Rheum (Polygonaceae) [33] and the Ligularia-Cremanthodium-Parasenecio complex (Asteraceae) [2]. An alternative explanation for the shared morphology between Nomocharis and Lophophorum is hybridization. However, the monophyly of Lophophorum is supported by both ITS and plastid phylogenies (Figs. 3 and 4). Thus, convergence seems to better explain the morphological similarities and supports habitat specialization of Nomocharis and the N-N lilies within the H-D Mountains and QTP.

Detecting the environmental drivers of convergence remain beyond the scope of this study. However, it is noteworthy that many alpine plant groups exhibit floral traits that are well-adapted to the frequent but unpredictable rains experienced in alpine habitats [3436]. For example, the nodding flower orientation is thought to have evolved to avoid pollen damage and nectar dilution by rainfall [30, 31, 37, 38]. Floral orientation may also be strongly affected by niche features such as the presence and abundance of various types of pollinators. In particular, the horizontal orientation may increase the precision of pollen transfer in bilaterally symmetrical flowers (e.g. Lilium and Nomocharis) under some pollination syndromes [35, 36, 39]. However, morphological convergence among alpine plants may also be strongly affected by understudied environmental interactions, such as with the intense solar radiation experienced during the daytime in alpine areas or the cold night time temperatures [31]. Overall, morphological convergence within the QTP and H-D Mountains habitats is likely linked to the extreme morphological divergence between QTP and H-D Mountains endemics and their widespread relatives. Thus, morphological convergence among QTP and H-D Mountains species of Lilium-Nomocharis and within other plant groups merits more attention in future studies.


Our ITS and plastid gene trees reveal several signatures of possible hybridization. In particular, the gene trees exhibit incongruence. In the ITS phylogeny, Nomocharis and the N-N lilies form a clade in the ITS tree (Fig. 3) that is sister to Lilium sect. Liriotypus. This is in contrast to the plastid phylogeny, which shows poor resolution of Nomocharis and the N-N lilies and places them among species of sects. Sinomartagon, Martagon (Fig. 4). Incongruence between nuclear and plastid and nuclear gene trees is known to result from hybridization, but can also result from incomplete lineage sorting, which is common among vascular plants, and horizontal gene transfer, which is not [40, 41].

Another signature of hybridization may be the strong geographic clustering observed in the plastid phylogeny (Fig. 4) among clades, which are distantly related in the nuclear phylogeny (Fig. 3, Additional file 1: Figure S1). The sympatry of clades with closely related plastid genomes is consistent with secondary contact. Moreover, hybridization in Lilium-Nomocharis is most likely to occur among species that occur within reasonably close proximity due to the limited dispersability of seeds [42] and typically also of pollen via wind or pollinators [43].

If hybridization did occur between Nomocharis (including N-N lilies) and sympatric Lilium, it must have occurred following the evolution of the latter, ca. 12 Mya (Fig. 7). If the dates in the plastid phylogeny can be taken to represent the times of contact, then hybridization events occurred in Nomocharis 5.73 Mya with Sinomartagon and 4.85 Mya with Leucolirion species. These events seem to post-date late orogenies of the QTP ca. 7 Mya and pre-date uplift of the H-D Mountains, in the late Neogene (ca. 3.4 Mya, [25, 26]). However, 95 % CIs for the dates include the orogenic periods (Fig. 6) and may also be consistent with ecological expansion of some Nomocharis species into new elevational ranges (Fig. 7e).


Lilium-Nomocharis exhibits complex phylogenetic relationships typical of a pattern in which QTP and H-D Mountains endemic, morphologically and ecologically distinct vascular plant groups such as Nomocharis, are included within widespread ones, such as Lilium. Our phylogenetic results show that Nomocharis itself is paraphyletic and includes some species traditionally classified as Lilium; here, the N-N clade. Species of the N-N clade exhibit typical Lilium morphology, which distinguishes them from the Nomocharis species. Features characteristic of Nomocharis, such as horizon oriented and flat/open flowers are probably ancestral to the group, and evolved before the uplift of the QTP. However, such features may have enabled the invasion of the QTP and, later, the H-D Mountains by Nomocharis and should be the subject of future studies. Despite their differences, Nomocharis and the N-N clade have probably evolved some similarities due to differently timed expansions into diverse elevational habitats. Our phylogenetic results also show some circumstantial evidence for hybridization in among traditional Lilium and Nomocharis species, and that may help to explain the complex phylogenetic relationships within the Lilium-Nomocharis complex.


Plant materials

We reconstructed a molecular phylogeny of Lilium and Nomocharis using nuclear ITS and 294 total accessions, of which 67 were obtained from GenBank, 227 were collected with necessary permissions by the author, of which 30 were newly sequenced for this study (Table 2, Additional file 8: Table S1). Note that only 90 accessions used for our phylogenetic reconstruction have been sequenced for all plastid markers and ITS (Table 2, Additional file 8: Table S1). For molecular phylogenetic reconstructions of plastid DNA, we focused our sampling efforts on Nomocharis and its Lilium allies; namely Lilium species that are geographically and/or evolutionarily close to Nomocharis. Of particular note, we sampled L. henrici Franchet, L. xanthellum F. T. Wang & T. Tang, L. saccatum S. Y. Liang that are endemic to the H-D Mountains and have been sparsely sampled in previous studies. Among Nomocharis species, only N. synaptica Sealy, which is native to India, was not sampled. Additionally, we included representative species of Lilium from across the geographic and phylogenetic distribution of the genus. Altogether, for the plastid phylogeny we sampled 14 Nomocharis accessions representing seven of eight species, thirteen Lilium species for their geographic or evolutionary proximity to Nomocharis, and 29 additional Lilium species (Table 2). We selected representative accessions of other genera from within the Lilieae tribe as outgroups including two each of Notholirion, Cardiocrinum and Fritillaria (see [44]). Of the total 360 sequences that we used in this study, two hundred and sixty-five are new to our study, and these have collection, voucher, and Genbank accession information provided in Table 2. We have deposited downstream sequencing data, namely alignments and phylogenetic trees, in TreeBase (Submission number: 17567).

Table 2 Materials and GenBank accession numbers of five chloroplast makers and accession information

We surveyed the morphology of Nomocharis, its close allies, and major lineages throughout Lilium. In particular, we used photographs of specimens observed in the field, field collected materials, and greenhouse specimens to assess macromorphological traits of 14 species of Nomocharis and closely related species of Lilium. To evaluate the same characters more broadly in 10 major lineages of Lilium (based on our phylogenetic results) we examined preserved specimens available to us, utilized the Chinese Virtual Herbarium, and obtained data from the literature (e.g., Flora of China [20]).

DNA extraction, Polymerase Chain Reaction (PCR) and sequencing

We selected the nuclear marker ITS and the cpDNA regions trnL-F, rbcL, matK, rpl32-trnL and psbA-trnH to reconstruct the molecular phylogeny of Lilium-Nomocharis. We chose the five cpDNA makers because three of them have been proposed as DNA barcodes for their high resolution and amplification success [45], and the other two have shown suitable variation in preliminary analyses (data not shown). For PCR amplifications of nuclear and plastid markers, we used total DNA extractions from fresh or silica gel-dried leaf tissue using a modified cetyltrimethyl-ammonium bromide (CTAB) protocol by Doyle and Doyle [46] or the Plant Genomic DNA Kit (TIANGEN Biotech, Beijing, China). We amplified all six markers using the primers listed in Table 3. All PCR reactions were performed with 50 ng genomic DNA in 20 μl reactions in a GeneAmp PCR System 9700 (Applied Biosystems, USA). The ITS reactions were performed using the following thermocycler protocol: 94 °C denaturation for 2 min; 35 cycles of 94 °C denaturation for 30 s, 55 °C primer annealing for 30 s, and 72 °C extension for 60 s; and a final extension of 72 °C for 10 min. For the plastid markers, the amplification conditions were the same except that primer annealing was performed at 52 °C for 45 s each cycle. Our amplified PCR products were sent to Invitrogen Biotech Co. Ltd. (Shanghai, China) for purification and sequencing, which was done on an ABI-3730XL DNA sequencer. For each sequenced accession, forward and reverse sequencing reactions were performed for increased coverage. Sequencing of the psbA-trnH spacer failed in two species, Nomocharis basilissa and Lilium nepalense, due to homopolymers at ~200 bp from the 5’ end. Thus, all data for this marker for these two species was considered missing (i.e., '?’, [47]) in downstream phylogenetic analyses.

Table 3 Primers and sequences statistics of nuclear and chloroplast makers used in present study

Molecular analysis

We aligned our DNA sequences using ClustalX [48] and then by eye in MEGA4.0 [49] following the guidelines of Morrison [50]. We trimmed the sequences to the limits of the ITS and the plastid regions, respectively, by comparing with examples deposited in Genbank. We positioned gaps to minimize nucleotide mismatches. We combined the five cpDNA markers into a single dataset, and all six aligned, and curated datasets were used to calculate uncorrected pairwise nucleotide differences in PAUP* version 4.0b10 [51]. Our nuclear ITS dataset contained a total of 294 accessions, inclusive of our eight outgroups. The ITS matrix contained 673 characters of which 398 were variable and 271 were parsimony-informative. There were 90 accessions for which sequences of all chloroplast markers were available, including for six outgroups. Details of the five chloroplast makers are presented in Table 3. The combined cpDNA alignment was 3429 bp long and contained 336 variables sites, of which 218 (or 6.3 %) were parsimony informative.

For phylogenetic analyses, we combined all five plastid sequences, because chloroplast genes have shared evolutionary histories within the chloroplast genome and because they do not recombine. We treated the ITS dataset independently. Bayesian phylogenetic analyses of the combined chloroplast dataset and the ITS dataset were conducted using MrBayes version 3.1.2 [52] with the GTR+ G + I and GTR+ G models of nucleotide substitution, respectively. These models were selected under the Akaike information criterion (AIC) using MrModeltest version 2.2 [53]. For each of the two datasets, we performed two simultaneous Bayesian analyses that started from a random tree and ran for 10 million generations with sampling every 1000 generations. Within each simultaneous run, four independent MCMC chains were used and the temperature increment between chains was adjusted to 0.2 based on mixing observed in preliminary analyses. Variation in likelihood scores was examined graphically for each independent run using Tracer 1.4 [54] and was used to determine apparent stationarity. Based on observations in Tracer, the first 25 % (2500) of posterior trees were discarded from each run as “burn-in” and posterior probabilities (pp) of clades were calculated from the remaining trees. Following burnin, we selected the best tree from among the simultaneous analyses of the plastid and ITS dataset, independently, using maximum clade credibility.

Maximum parsimony (MP) analyses of the ITS and the combined chloroplast makers were carried out using PAUP* [51]. Characters were treated as unordered and unweighted. A heuristic search was performed with 1000 replicate analyses, random stepwise addition of taxa, tree-bisection-reconnection (TBR) branch swapping, and maximum trees set to 50,000. We summarized the resulting equally parsimonious topologies using majority-rule consensus and calculated bootstrap values from one million replicate analyses using fast stepwise addition of taxa. We retained the bootstrap values for clades consistent with the majority-rule consensus tree.

We carried out topological testing using Kishino-Hasegawa (KH) tests in PAUP*, because KH tests are known to exhibit very low type I error rates [55]. To perform the tests, we used a reduced dataset, which consisted of one sequence for each major evolutionary lineage that was mutually represented in the plastid and nuclear gene trees (Additional file 7: Figure S7). We confirmed that the selected samples produced the same arrangements of evolutionary lineages as the entire plastid and nuclear alignments by generating maximum likelihood (ML) trees using the GTR+ G + I and GTR+ G models, respectively (data not shown). Major lineages were manually organized into plastid and nuclear cladograms in Mesquite [56] (Additional file 7: Figure S7). The reduced alignments plus the cladograms were loaded into PAUP* for performing the KH tests. Specifically, we used the tests to determine if each tree represented a significantly better fit for the dataset from which it was reconstructed compared to tree resulting from the other dataset. We performed the KH tests under the GTR+ G + I and GTR+ G models for the plastid and nuclear datasets using a normal test distribution.

Statistical parsimony network

We expected that strictly bifurcating trees may not completely describe the evolutionary relationships within Lilium-Nomocharis, because hybridization in Lilium-Nomocharis has been postulated [13, 17, 57] and incomplete lineage sorting has been detected in many plant lineages [40]. Therefore, we used the statistical parsimony network approach implemented in TCS v.1.21 [58] to further evaluate evolutionary relationships within the Lilium-Nomocharis complex using the combined chloroplast sequences. We built the parsimony network using eighty-four accessions sequenced for all cpDNA markers except psbA-trnH, which was missing data for two taxa (see above). We tested whether removal of psbA-trnH would change relationships among species, by reconstructing a bifurcating plastid phylogeny without the marker, and it showed no differences compared to the tree constructed using whole dataset (results not shown). For the network analysis, we considered each indel as a single mutation event, and all indels were reduced to single characters (arbitrarily A or T) in a final alignment. The resulting plastid matrix was 3037 characters in length and contained 66 plastid haplotypes representing 84 accessions of Lilium-Nomocharis. We eliminated loops from the parsimony based on the principle that haplotypes with interior positions in the network are assumed to be ancestral [59].

Divergence estimation

Molecular dating in Liliales has been previously performed using distantly related fossils [60], calibrations from previous studies [44, 61], and single calibration points [17]. In particular, Bremer [60] dated nodes in the monocot phylogeny using fossils closely related to palms, aroids, grasses, and cattails and found that Liliales evolved approximately 112 Mya and began diversifying 82 ± 10 Mya. Deriving calibration points from Bremer [60], Patterson and Givnish [44] inferred the divergence time of the tribe Lilieae as 12 Mya and Vinnersten and Bremer [61] concluded that the monophyletic lineage comprised of Lilium, Nomocharis and Fritillaria diverged 6 ± 2.9 Ma. Gao et al. [17] provided a detailed review of Liliales fossils and performed dating using a single, reliable fossil of Smilax, Smilax wilcoxensis Berry [62], to calibrate the divergence between Liliaceae and Smilaceae. Their results showed that Lilieae evolved approximately 16mya. Despite these efforts, it has been widely discussed and shown that single calibration points and caibrations derived from prior studies lead to less reliable, and often younger, clade ages [6365].

We sought to more rigorously date events in Lilium-Nomocharis by applying two calibration points for dating analyses in BEAST (Additional file 2: Figure S2) [66, 67]. For one calibration, we constrained the divergence time of Liliaceae and Smilacaceae using Smilax wilcoxensis. In brief, Smilax wilcoxensis is known from the early Eocene (48.6–55.8 Mya) of the Tennessee Wilcox Formation [62, 68], which is assigned a relative age based on pollen [69, 70]. Specifically, we calibrated the Liliaceae-Smilacaceae node using a uniform prior with a lower bound (paleontologically upper) of 48.6 Mya and an upper bound of 131 Mya. Thus, we asserted our belief that Smilacaceae cannot be younger than Smilax wilcoxensis or older than the Barremian (i.e., 131 Mya), from which the oldest flowering plant fossil is known [71]. For the second calibration, we used Ripogonum tasmanicum Conran, et al. [72] to constrain the age of the ancestor of the monotypic Ripogonanceae and Philesiaceae (following Angiosperm Phylogeny Website, [73]). Ripogonum tasmanicum is reported from the Tasmanian Macquarie Harbour Formation [72], which is approximately 51–52 million years old based on a foraminiferal index [74]. Thus, we constrained the Ripogonanceae and Philesiaceae split using a uniform prior with a lower bound of 51 Mya and an upper bound of 131 Mya. The prior asserts our belief that Ripogonaceae cannot be younger than its fossil or older than the earliest known flowering plant.

The two fossils facilitated establishing calibration points that were well outside of the Nomocharis-Lilium complex. Therefore, we applied these two calibrations to infer the split between Lilium and Fritillaria using a dataset comprised of three cpDNA markers (aptF-H, matK and rbcL, see Additional file 9: Table S2, Additional file 2: Figure S2) that included 45 representative Liliales species and more than 3000 bp [75]. We applied the result mean and 95 % Highest Posterior Density (HPD) to constrain the Lilium and Fritillaria node using a normal prior distribution in an analysis of our plastid dataset. We take these results (Additional file 2: Figure S2) to be our best estimates of ages within Lilium-Nomocharis. More vetted fossils closer to Lilium may eliminate the need for the second dating step in the future.

Divergence time estimations were performed using BEAST ver. 1.5.3 [67] identically for the cpDNA and ITS datasets. The normal prior distribution on the age of the Lilium stem node (i.e., the split of Lilium and Fritillaria) was set using a mean of 14.92 Mya and a standard deviation of 2.5. The chosen standard deviation gave a 95 % HPD of 10.81-19.03 Ma, which was slightly narrower than the actual result of 6.32–25.71 Ma. A likelihood ratio test in PAUP 4.10b [51] rejected strict clocks for both datasets (P < 0.01), therefore we used an uncorrelated lognormal (UCLN), relaxed clock [76]. We used the GTR + G + I and GTR + G models of nucleotide substitution for combined plastid and nuclear ITS dataset, respectively. For the distribution of divergence times, a pure birth branching process (Yule model) was chosen as a prior. BEAST analyses were run on the Cyberinfrastructure for Phylogenetic Research (CIPRES) Science Gateway ( We ran two independent Markov chains, each for 50,000,000 generations, initiated with a random starting tree, and sampled every 1000 generations. The first 20 % of sampled trees from all runs were discarded as burn-in based on visual inspection in Tracer version 1.4 [54].

Ancestral Area Reconstructions (AAR)

We used the Bayesian Binary method (BBM) in Reconstruct Ancestral States in Phylogenies 2.1b (RASP 2.0) [7779] to reconstruct the biogeographic history of Lilium-Nomocharis on the ITS consensus phylogeny constructed from BEAST trees. Based on prior studies (e.g., [20, 80]) three areas of endemism were recognized: Qinghai-Tibetan Plateau (QTP, A), H-D Mountains (HDM, B), the geographic region now covered by Sino-Japanese Forest subkingdom (SJFS, C; A-C stand for each region in the RASP analyses, Table 2). We compared BBM results to results from Lagrange, which implements a likelihood method and the Dispersal-Extinction-Cladogenesis (DEC) model [81]. In Lagrange, we set migration probabilities among the three areas of endemism to 1.0 throughout time and did not limit the number of areas that a widespread taxon could occupy (Additional file 10: Table S3). We allowed Lagrange to estimate the extinction and dispersal parameters required for the DEC model.

Ancestral state reconstruction (ASR)

We reconstructed the ancestral states for four, variable macro morphological characters and the habitat characteristic, elevation, in Lilium-Nomocharis. We selected variable macromorphological characters with states that could be evaluated with confidence given the coarse availability of specimen data (see Taxon sampling above). Specifically, we performed reconstructions for corolla shape, flower orientation, the ratio of stigma versus stamen length, and leaf arrangement (Additional file 11: Table S4). We selected these characters from among other plausible ones, because they have previously been used to delimit species within Lilium and Nomocharis [19, 20, 80] but they have not been previously considered within a phylogenetic framework. For corolla shape, we coded species as having flat or open flowers, campaniform or bell shaped flowers, recurved, funnel or trumpet shaped, or bowl-shaped. Flower orientation states were coded as nodding, horizon, and up (i.e., upward facing). For stigma-stamen ratio, we coded states as being greater than 1.25, less than 0.75, or between 0.75 and 1.25. Using these ranges for stigma-stamen ratios enabled us to code species visually. Leaf arrangement was coded as being alternate or whorled. The whorled leaf character was assigned to species that have 3+ leaves arising from a single node and species with scattered leaves arising asynchronously [82]. For elevation, we acquired information from floras and specimen records on GBIF ( We treated elevation as categorical by using 1000 ft. increments for our discrete character states.

To reconstruct the ancestral character states we used BBM in RASP, which is not limited to historical biogeographic applications. We performed the reconstructions of ancestral morphological states across the dated ITS consensus tree resulting from the BEAST analysis and using the character matrices presented in Additional file 9: Table S2. We modified the BEAST consensus tree using TreeGraph 2.0 [83] by pruning outgroups and collapsing the major clades except Nomocharis. We did this to avoid confounding the issue with outgroups, which were not completely sampled or studied, and to simplify the reconstructions for less well sampled clades outside of Nomocharis. Branch length and divergence time information were preserved. The Bayesian analyses in RASP were carried out using default settings except that we ran the analyses for 1,000,000 MCMC generations and used the F81 + G model for changes between states.


  1. Myers N, Mittermeier RA, Mittermeier CG, da Fonseca GAB, Kent J. Biodiversity hotspots for conservation priorities. Nature. 2000;403:853–8.

    Article  CAS  PubMed  Google Scholar 

  2. Liu JQ, Wang YJ, Wang AL, Ohba H, Abbott RJ. Radiation and diversification within the Ligularia-Cremanthodium-Parasenecio complex (Asteraceae) triggered by uplift of the Qinghai-Tibetan Plateau. Mol Phylogenet Evol. 2006;38:31–49.

    Article  CAS  PubMed  Google Scholar 

  3. Yue JP, Sun H, Al-Shehbaz IA, Li JH. Support for an expanded Sloms-Laubachia (Brassicaceae): evidence from sequence of chloroplast and nuclear genes. Ann Mo Bot Gard. 2006;93:402–11.

    Article  Google Scholar 

  4. Peterson A, Levichev IG, Peterson J. Systematics of Gagea and Lloydia (Liliaceae) and infrageneric classification of Gagea based on molecular and morphological data. Mol Phylogenet Evol. 2008;46(2):446–65.

    Article  CAS  PubMed  Google Scholar 

  5. Wang YJ, Li XX, Hao G, Li JQ. Molecular phylogeny and biogeography of Androsace (Primulaceae) and the convergent evolution of cushion morphology. Acta Phytotaxonomica Sinica. 2004;42(6):481–99.

    Google Scholar 

  6. Liu JQ, Chen ZD, Lu AM. A preliminary study of the phylogeny of the Swertiinae based on ITS data (Gentianaceae). Irsael J Plant Sci. 2001;49:345–9.

    Google Scholar 

  7. Friesen N, Fritsch RM, Pollner S, Blattner F. Molecular and morphological evidence for an origin of the aberrant genus Milula within Himalayan species of Allium (Alliaceae). Mol Phylogenet Evol. 2000;17:209–18.

    Article  CAS  PubMed  Google Scholar 

  8. Tang HG, Meng LH, Ao SM, Liu JQ. Origin of the Qinghai-Tibetan Plateau endemic Milula (Liliaceae): further insights from karyological comparisons with Allium. Caryologia. 2005;58(4):320–31.

    Article  Google Scholar 

  9. Liu JQ, Tian B. Origin, evolution, and systematics of Himalaya endemic genera. Newslett Himalayan Botany. 2007;40:20–7.

    Google Scholar 

  10. Raudnitschka D, Hensen I, Oberprieler C. Introgressive hybridization of Senecio hercynicus and S. ovatus (Compositae, Senecioneae) along an altitudinal gradient in Harz National Park (Germany). Syst Biodivers. 2007;5(3):333–44.

    Article  Google Scholar 

  11. Carling MD, Thomassen HA. The Role of Environmental Heterogeneity in Maintaining Reproductive Isolation between Hybridizing Passerina (Aves: Cardinalidae) Buntings. Int J Ecol. 2012;2012:295463.

    Article  Google Scholar 

  12. Balfour B. The Genus Nomocharis. Bot J Scotl. 1918;27(3):273–300.

    Google Scholar 

  13. Sealy JR. Nomocharis and Lilium. Kew Bull. 1950;5(2):273–97.

    Article  Google Scholar 

  14. Sealy JR. A revision of the genus Nomocharis Franchet. Bot J Linn Soc. 1983;87:285–323.

    Article  Google Scholar 

  15. Liang SY. Studies on the genus Nomocharis (Liliaceae). Bull Botanical Res. 1984;4(3):163–78.

    CAS  Google Scholar 

  16. Gao YD, Hohenegger M, Harris A, Zhou SD, He XJ, Wan J. A new species in the genus Nomocharis Franchet (Liliaceae): evidence that brings the genus Nomocharis into Lilium. Plant Syst Evol. 2012;298:69–85.

    Article  Google Scholar 

  17. Gao YD, Harris A, Zhou SD, He XJ. Evolutionary events in Lilium (including Nomocharis, Liliaceae) are temporally correlated with orogenies of the Q-T plateau and the Hengduan Mountains. Mol Phylogenet Evol. 2013;68(3):443–60.

    Article  PubMed  Google Scholar 

  18. Hayashi K, Kawano S. Molecular systematics of Lilium and allied genera (Liliaceae): phylogenetic relationships among Lilium and related genera based on the rbcL and matK gene sequence data. Plant Species Biol. 2000;15:73–93.

    Article  Google Scholar 

  19. Liang SY. The genus Lilium L. In: Wang FZ, Tang J, editors. Flora Reipublicae Popularis Sinicae, Anagiospermae, Monocotyledoneae Liliaceae (I), vol. 14. Beijing: Science Press; 1980. p. 116–57.

    Google Scholar 

  20. Liang SY, Tamura M. Lilium L. In: Wu ZY, Raven PH, editors. Flora of China, vol. 24. Beijing; St. Louis: Science Press; Missouri Botanical Garden Press; 2000. p. 135–59.

    Google Scholar 

  21. Gao YD, Zhou SD, He XJ. Lilium yapingense (Liliaceae), a new species from Yunnan, China, and its systematic significance relative to Nomocharis. Ann Bot Fenn. 2013;50:187–94.

    Article  Google Scholar 

  22. Wu ZY, Wu SG. A proposal for a new floristic kingdom (realm) - the E. Asiatic kingdom, its delimitation and characteristics. In: Zhang AL, Wu SG, editors. Proceedings of the First International Symposium on Floristic Characteristics and Diversity of East Asian Plants. Beijing and Berlin: Higher Education Press, Springer-Verlag Heidelberg; 1996. p. 3–42.

    Google Scholar 

  23. Comber HF. A new classification of genus Lilium. Royal Horticult Soc Liliy Year Book. 1949;13:85.

    Google Scholar 

  24. Haw SG. The Lilies of China. Portland: Timber Press; 1986.

    Google Scholar 

  25. Chen FB. H-D event: an important tectonic event of the late Cenozoic in Eastern Asia. Mt Res. 1992;10(4):195–202.

    Google Scholar 

  26. Chen FB. Second discussion on the H-D movement. Mt Res. 1996;17(3–4):14–22.

    Google Scholar 

  27. Harrison TM, Copeland P, Kidd WSF, Yin A. Raising Tibet. Science. 1992;255:1663–70.

    Article  CAS  PubMed  Google Scholar 

  28. Li JJ, Shi YF, Li BY. Uplift of the Qinghai-Xizang (Tibet) Plateau and Global Change. Lanzhou, China: Lanzhou University Press; 1995.

    Google Scholar 

  29. Chapman MA, Hiscock SJ, Filatov DA. Genomic divergence during speciation driven by adaptation to altitude. Mol Biol Evol. 2013;30(12):2553–67.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  30. Mao YY, Huang SQ. Pollen resistance to water in 80 angiosperm species: flower structures protect rain-susceptible pollen. New Phytol. 2009;183(3):892–9.

    Article  PubMed  Google Scholar 

  31. Wang Y, Meng LL, Yang YP, Duan YW. Change infloral orientation in Anisodus luridus (Solanaceae) protects pollen grains and facilitates development of fertilized ovules. Am J Bot. 2010;97:1618–24.

    Article  PubMed  Google Scholar 

  32. Wang YJ, Liu JQ. Phylogenetic analyses of Saussurea sect. Pseudoeriocoryne (Asteraceae: Cardueae) based on chloroplast DNA trnL-F sequences. Biochem Sys Ecol. 2004;32:1009–21.

    Article  CAS  Google Scholar 

  33. Wang AL, Yang MY, Liu JQ. Molecular phylogeny, recent radiation and evolution of gross morphology of the Rhubarb genus Rheum (Polygonaceae) inferred from chloroplast DNA trnL-F sequences. Ann Bot. 2005;96:489–98.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  34. Köener C. Alpine plant life: functional plant ecology of high mountain ecosystems. Berlin: Springer; 2003.

    Google Scholar 

  35. Ushimaru A, Dohzono I, Takami Y, Hyodo F. Flower orientation enhances pollen transfer in bilaterally symmetrical flowers. Oecologia. 2009;160:667–74.

    Article  PubMed  Google Scholar 

  36. Ushimaru A, Hyodo F. Why do bilaterally symmetrical flowers orient vertically? Flower orientation in fluences pollinator landing behaviour. Evol Ecol Res. 2005;7:151–60.

    Google Scholar 

  37. Huang SQ, Takahashi Y, Dafni A. Why does the flower stalk of Pulsatilla cernua (Ranunculaceae) bend during anthesis? Am J Bot. 2002;89:1599–603.

    Article  PubMed  Google Scholar 

  38. Sun JF, Gong YB, Renner SS, Huang SQ. Multifunctional bracts in the Dove tree Davidia involucrate (Nyssaceae: Cornales): rain protection and pollinator attraction. Am Nat. 2008;117:119–24.

    Article  Google Scholar 

  39. Fenster CB, Armbruster WS, Dudash MR. Specialization of flowers: is floral orientation an overlooked step? New Phytol. 2009;183(3):502–6.

    Article  PubMed  Google Scholar 

  40. Maddison W, Knowles LL. Inferring phylogeny despite incomplete lineage sorting. Syst Biol. 2006;55(1):21–30.

    Article  PubMed  Google Scholar 

  41. Richardson AO, Palmer JD. Horizontal gene transfer in plants. J Exp Bot. 2007;58(1):1–9.

    Article  CAS  PubMed  Google Scholar 

  42. Minami S, Azuma A. Various flying modes of wind-dispersal seeds. J Theor Biol. 2003;225(1):1–14.

    Article  PubMed  Google Scholar 

  43. Hiramatsu M, Ii K, Okubo H, Huang KL, Huang CW. Biogeography and origin of Lilium longiflorum and L. formosanum (Liliaceae) endemic to the Ryukyu Archipelago and Taiwan as determined by allozyme diversity. Am J Bot. 2001;88(7):1230–9.

    Article  CAS  PubMed  Google Scholar 

  44. Patterson TB, Givnish TJ. Phylogeny, concerted convergence, and phylogenetic niche conservatism in the core Liliales: Insights from rbcL and ndhF sequence data. Evolution. 2002;56(2):233–52.

    Article  CAS  PubMed  Google Scholar 

  45. Shaw J, Lickey EB, Schilling EE, Small RL. Comparison of whole chloroplast genome sequences to choose noncoding regions for phylogenetic studies in angiosperms: the tortoise and the hare III. Am J Bot. 2007;94(3):275–88.

    Article  CAS  PubMed  Google Scholar 

  46. Doyle JJ, Doyle JL. A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem Bull. 1987;19:11–5.

    Google Scholar 

  47. Wiens JJ, Morrill MC. Missing data in phylogenetic analysis: reconciling results from simulations and empirical data. Syst Biol. 2011;60(5):719–31.

    Article  PubMed  Google Scholar 

  48. Thompson JD, Gibson TJ, Plewniak F, Jeanmougin F, Higgins DG. The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 1997;25:4876–82.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  49. Tamura K, Dudley J, Nei M, Kumar S. MEGA4: Molecular Evolutionary Genetics Analysis (MEGA) software version 4.0. Mol Biol Evol. 2004;24:1596–9.

    Article  Google Scholar 

  50. Morrison DA. A framework for phylogenetic sequence alignment. Plant Syst Evol. 2009;282:127–49.

    Article  Google Scholar 

  51. Swofford DL. PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). 4th ed. Sunderland, Massachusetts: Sinauer Associates; 2003.

    Google Scholar 

  52. Ronquist F, Huelsenbeck JP. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003;19(12):1572–4.

    Article  CAS  PubMed  Google Scholar 

  53. Nylander JAA. MrModeltest 2.0. In: Department of Systematic Zoology. 20th ed. Uppsala: EBC, Uppsala University; 2004.

    Google Scholar 

  54. Rambaut A, Drummond AJ. Tracer v1.4. 2007.

    Google Scholar 

  55. Susko E. Tests for two trees using likelihood methods. Molecular Biology and Evolution. 2014;31(4):1029–39.

    Article  CAS  PubMed  Google Scholar 

  56. Maddison WP, Maddison DR. Mesquite: a modular system for evolutionary analysis. In: 2.75 edn; 2011.

  57. Douglas NA, Wall WA, Xiang QY, Hoffman WA, Wentworth TR, Gray JB, et al. Recent vicariance and the origin of the rare, edaphically specialized Sandhills lily, Lilium pyrophilum (Liliaceae): evidence from phylogenetic and coalescent analyses. Mol Ecol. 2011;20:2901–15.

  58. Clement M, Posada D, Crandall KA. TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000;9:1657–60.

    Article  CAS  PubMed  Google Scholar 

  59. Crandall KA, Templeton AR. Empirical tests of some predictions from coalescent theory with applications to intraspecific phylogeny reconstruction. Genetics. 1993;134:959–69.

    CAS  PubMed Central  PubMed  Google Scholar 

  60. Bremer K. Early Cretaceous lineages of monocot flowering plants. Proc Natl Acad Sci U S A. 2000;97(9):4707–11.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  61. Vinnersten A, Bremer K. Age and biogeography of major clades in Liliales. Am J Bot. 2001;88:1695–703.

    Article  CAS  PubMed  Google Scholar 

  62. Berry EW. Revision of the Lower Eocene Wilcox flora of the southeastern states: With descriptions of new species, chiefly from Tennessee and Kentucky, vol. 156. Washington, DC: US Government Printing Office; 1930.

  63. Sauquet H, Ho SYW, Gandolfo MA, Jordan GJ, Wilf P, Cantrill DJ, et al. Testing the impact of calibration on molecular divergence times using a Fossil-Rich Group: The Case of Nothofagus (Fagales). Syst Biol. 2012;61(2):289–313.

  64. Rambaut A, Bromham L. Estimating divergence dates from molecular sequences. Mol Biol Evol. 1998;15(4):442–8.

    Article  CAS  PubMed  Google Scholar 

  65. Graur D, Martin W. Reading the entrails of chickens: molecular timescales of evolution and the illusion of precision. Trends Genet. 2004;20(2):80–6.

    Article  CAS  PubMed  Google Scholar 

  66. Drummond AJ, Rambaut A. BEAST: Bayesian Evolutionary Analysis Sampling Trees v1.3. 2003, 2012.

  67. Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:214.

    Article  PubMed Central  PubMed  Google Scholar 

  68. Dilcher DL, Lott TA. A middle Eocene fossil plant assemblage (Powers Clay Pit) from western Tennessee. Florida Mus Nat Hist Bull. 2005;45:1–43.

    Google Scholar 

  69. Potter F, Dilcher D. Biostratigraphic analysis of Eocene clay deposits in Henry County, Tennessee. Biostratigraph Fossil Plants. 1980;211:225.

    Google Scholar 

  70. Tschudy RH. Stratigraphic distribution of significant Eocene palynomorphs of the Mississippi embayment. Washington, DC: US Government Printing Office; 1973.

  71. Sun G, Dilcher DL, Wang H, Chen Z. A eudicot from the Early Cretaceous of China. Nature. 2011;471(7340):625–8.

    Article  CAS  PubMed  Google Scholar 

  72. Conran JG, Carpenter RJ, Jordan GJ. Early Eocene Ripogonum (Liliales: Ripogonaceae) leaf macrofossils from southern Australia. Aust Syst Bot. 2009;22(3):219–28.

    Article  Google Scholar 

  73. Stevens PF. Angiosperm Phylogeny Website. Version 12 []

  74. Carpenter RJ, Jordan GJ, Hill RS. A toothed Lauraceae leaf from the Early Eocene of Tasmania. Austr Int J Plant Sci. 2007;168(8):1191–8.

    Article  Google Scholar 

  75. Kim JS, Hong JK, Chase MW, Fay MF, Kim JH. Familial relationships of the monocot order Liliales based on a molecular phylogenetic analysis using four plastid loci: matK, rbcL, atpB and atpF-H. Bot J Linn Soc. 2013;172(1):5–21.

    Article  Google Scholar 

  76. Drummond AJ, Ho SYW, Phillips MJ, Rambaut A. Relaxed phylogenetics and dating with confidence. PLoS Biol. 2006;4(5):e88.

    Article  PubMed Central  PubMed  Google Scholar 

  77. Yu Y, Harris A, He XJ. S-DIVA (statistical dispersal-vicariance analysis): a tool for inferring biogeographic histories. Mol Phylogenet Evol. 2010;56:848–50.

    Article  PubMed  Google Scholar 

  78. Yu Y, Harris A, He XJ. RASP (Reconstruct Ancestral State in Phylogenies) 2.0 beta. In:; 2012.

  79. Yu Y, Harris AJ, He XJ. A novel Bayesian method for reconstructing geographic ranges and ancestral states on phylogenies. 2011.

    Google Scholar 

  80. McRae EA. Lilies: a guide for growers and collectors. Portland: Timber Press; 1998.

    Google Scholar 

  81. Ree RH, Smith SA. Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis. Syst Biol. 2008;57(1):4–14.

    Article  PubMed  Google Scholar 

  82. Rutishauser R. Polymerous leaf whorls in vascular plants: developmental morphology and fuzziness of organ identities. Int J Plant Sci. 1999;160(S6):S81–103.

    Article  PubMed  Google Scholar 

  83. Stöver BC, Müller KF. TreeGraph 2: combining and visualizing evidence from different phylogenetic analyses. BMC Bioinformatics. 2010;11:7.

    Article  PubMed Central  PubMed  Google Scholar 

  84. Ho TN, Pringle JS. Gentianaceae. In: Wu ZY, Raven PH, editors. Flora of China, vol. 16. Beijing and St. Louis: Science Press and Missouri Botanical Garden Press; 1995. p. 1–140.

    Google Scholar 

  85. Levan A. Cytological studies in Allium, II Chromosome morphological contributions. Hereditas. 1932;16:257–94.

    Article  Google Scholar 

  86. Tian XM, Luo J, Wang AL, Mao KS, Liu JQ. On the origin of the woody buckwheat Fagopyrum tibeticum (=Parapteropyrum tibeticum) in the Qinghai-Tibetan Plateau. Mol Phylogenet Evol. 2011;61(2):515–20.

    Article  PubMed  Google Scholar 

  87. Shi Z, Chen YL, Chen YS, Lin YR, Liu SW, Ge XJ, et al. Asteraceae (Compositae). In: Wu ZY, Raven PH, editors. Flora of China, vol. 20–21. Beijing and St. Louis: Science Press and Missouri Botanical Garden Press; 2011. p. 1–8.

  88. Hu CM, Kelso S. Lysimachia Linnaeus. In: Wu ZY, Raven PH, editors. Flora of China, vol. 15. Beijing: Science Press; St. Louis: Missouri Botanical Garden Press; 1996. p. 39–78.

    Google Scholar 

  89. Schneeweiss GM, Schonswetter P, Kelso S, Niklfeld H. Complex biogeographic patterns in Androsace (Primulaceae) and related genera: evidence from phylogenetic analyses of nuclear internal transcribed spacer and plastid trnL-F sequences. Syst Biol. 2004;53(6):856–76.

    Article  PubMed  Google Scholar 

  90. Cheo TY, Lu L, Yang G, Al-Shenbaz I, Dorofeev V. Brassicaceae. In: Wu ZY, Raven PH, editors. Flora of China, vol. 8. Beijing, St. Louis: Science Press, Missouri Botanical Garden Press; 2001. p. 1–193.

    Google Scholar 

  91. Yue JP, Gu ZJ, Al-Shehbaz IA, Sun H. Cytological studies on the Sino-Himalayan endemic Solms-laubachia (Brassicaceae) and two related genera. Bot J Linn Soc. 2004;145(1):77–86.

    Article  Google Scholar 

  92. White TJ, Bruns T, Lee S, Taylor J. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. In: Innis MA, Gelfand DH, Shinsky JJ, White TJ, editors. PCR protocols: a guide to methods and applications. San Diego: Academic Press; 1990. p. 315–22.

    Google Scholar 

  93. Fay MF, Swensen SM, Chase MW. Taxonomic affinities of Medusagyne oppositifolia (Medusagynaceae). Kew Bull. 1997;52(1):111–20.

    Article  Google Scholar 

  94. Cuenoud P, Savolainen V, Chatrou LW, Powell M, Grayer RJ, Chase MW. Molecular phylogenetics of Caryophyllales based on nuclear 18S rDNA and plastid rbcL, atpB, and matK DNA sequences. Am J Bot. 2002;89(1):132–44.

    Article  CAS  PubMed  Google Scholar 

  95. Taberlet PGL, Pautou G, Bouvet J. Universal primers for amplification of three noncoding regions of chloroplast DNA. Plant Mol Biol. 1991;17:1105–9.

    Article  CAS  PubMed  Google Scholar 

  96. Hamilton MB. Four primer pairs for the amplification of chloroplast intergenic regions with intraspecific variation. Mol Ecol. 1999;8(3):521–3.

    CAS  PubMed  Google Scholar 

Download references


We thank Dr. You-Sheng Chen and Xiao-Hua Jin from Institution of Botany, Chinese Academy of Sciences for the help in field work and material collection. We also thank Qin-Qin Li, Cheng-Yang Liao, Chang-Bao Wang, Qiang Wang and Xiang-Guang Ma with the help in the field work. This work was supported by the National Natural Science Foundation of China (Grant Nos. 31270241, 31470009), and the Specimen Platform of China, Teaching Specimen’s sub-platform (

Author information

Authors and Affiliations


Corresponding authors

Correspondence to Yun-Dong Gao or Xing-Jin He.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

Conceived, designed, and performed the laboratory experiments: YDG. Performed analyses the data: YDG AJ-H. Contributed reagents/materials/analysis tools: YDG AJ-H XJH. Wrote the paper: YDG AJ-H. All authors read and approved the final manuscript.

Additional files

Additional file 1: Figure S1.

Reconstructed phylogenetic relationship of whole Lilium-Nomocharis based on Bayesian inferences of ITS dataset. Names of terminal clades based on Comber [23] and Liang [19]. The Sinomartagon I clade is highlighted.

Additional file 2: Figure S2.

Divergence dating of major clades of Liliales using two fossil calibrations (1 and 2).

Additional file 3: Figure S3.

Pictures from western China showing: a-c Lilium henrici var. henrici; d-f L. lophophorum; g-i L. saccatum; j-l L. yapingense.

Additional file 4: Figure S4.

Pictures from western China showing: a-e Lilium xanthellum with variations on tepal morphology within a same locality; g-i, flower of L. souilei, L. nanum and L. nepalense; j-l, habit of L. souilei, L. nanum and L. nepalense.

Additional file 5: Figure S5.

Pictures from western China showing Nomocharis: a-c, N. pardanthina; d-f, N. saluenensis; g-i, N. pardanthina f. punctulata.

Additional file 6: Figure S6.

Outer and inner tepals comparison in a, Lilium henrici; b, L. lophophorum; c-d, two types of L. xanthellum; e, L. yapingense; f, L. saccatum; g, Nomocharis saluenensis; h, N. pardanthina f. punctulata; i-j, two types of N. aperta (Zhongdian and Fugong, respectively); k, N. basilissa; l, N. gongshanensis; m, N. pardanthina; n, N. meleagrina.

Additional file 7: Figure S7.

Results of KH tests for ITS and combined plastid datasets.

Additional file 8: Table S1.

Sources of ITS sequence data.

Additional file 9: Table S2.

Genbank accessions used in diversification dating of major clades of Liliales.

Additional file 10: Table S3.

The matrix of model used in AAR analysis of LARANGE.

Additional file 11: Table S4.

Morphological character states used in ancestral state reconstruction.

Rights and permissions

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Gao, YD., Harris, A. & He, XJ. Morphological and ecological divergence of Lilium and Nomocharis within the Hengduan Mountains and Qinghai-Tibetan Plateau may result from habitat specialization and hybridization. BMC Evol Biol 15, 147 (2015).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: