Morphological and ecological divergence of Lilium and Nomocharis within the Hengduan Mountains and Qinghai-Tibetan Plateau may result from habitat specialization and hybridization
Gao et al. BMC Evolutionary Biology
Morphological and ecological divergence of Lilium and Nomocharis within the Hengduan Mountains and Qinghai-Tibetan Plateau may result from habitat specialization and hybridization
Yun-Dong Gao 0 2
AJ Harris 1
Xing-Jin He 0
0 Key Laboratory of Bio-Resources and Eco-Environment of Ministry of Education, College of Life Science, Sichuan University , Chengdu , China
1 Department of Botany, Oklahoma State University , 301 Physical Sciences, Stillwater, OK 74078-3013 , USA
2 Key Laboratory of Mountain Ecological Restoration and Bioresource Utilization & Ecological Restoration Biodiversity Conservation Key Laboratory of Sichuan Province, Chengdu Institute of Biology, Chinese Academy of Sciences , Chengdu 610041 , China
Background: 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. Results: 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. (Continued on next page)
(Continued from previous page)
Conclusions: 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 . 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, [2–8]). 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 , 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.
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) [12–15]. Currently,
there are eight recognized species of Nomocharis, of
Table 1 Morphologically distinctive plant species that are endemic to the QTP but phylogenetically indistinct (i.e., nested within)
Gentianaceae [9, 84]
Brassicaceae [3, 9, 90, 91]
QTP and adjacent regions Polygonaceae 
to the south and east
eastern North America and
temperate South America
and arctic areas in eastern
and central Asia and
Morphology of allies
Lomatogoniopsis Lomatogonium A. 2n = 12; Petals bearing one 2n = 18; Petals bearing
T. N. Ho & S. W. Liu Braun nectary each; Nectaries two nectarines each;
appendaged, not in pits Nectaries not
appendaged, in pits
n = 10; inflorescence n = 16 or multiples; Throughout the Northern
spicate; sepals fused over inflorescence umbellate; Hemisphere and in Africa,
1/3 or more of length sepals free or fused and Central and South
only at base America
woody; flowers bisexual
Parasenecio W. W.
Smith & J. Small
capitula discoid; roots
fruit capsule operculate
unique suite of characters unique suite of
fruit capsule opening
along longitudinal slits
Fig. 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
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 . 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 [18–20].
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.
A large ITS dataset confirmed the phylogentic position
of Nomocharis within Lilium and showed no major
Fig. 2 Pictures from western China showing Nomocharis: (a-c), N. basilissa; (d-f), N. farreri; (g-i), N. gongshanensis; (j-l), N. meleagrina
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
. 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).
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 . Cluster II (PP = 1.0, BS = 90 %)
contained Nomocharis and species of Lilium that occur
within the H-D Mountains and adjacent Himalayas.
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
Fig. 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 
Fig. 4 Maximum credibility tree resulting from a Bayesian analysis of combined plastid DNA. Clade names based on Comber  and Liang .
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
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).
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
Fig. 5 Parsimony network conducted by TCS  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
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).
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
Fig. 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
within the H-D Mountains and the adjacent
SinoJapanese 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
Fig. 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
postulated for the H-D Mountain uplift [25, 26].
LiliumNomocharis 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 nod
ding 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 .
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 . 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) , Pseudoeriocoryne (Asteraceae: Cardueae)
, Rheum (Polygonaceae)  and the
LigulariaCremanthodium-Parasenecio complex (Asteraceae) .
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 [34–36]. 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 . 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 
and typically also of pollen via wind or pollinators .
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
predate 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.
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 ). 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).
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 ).
DNA extraction, Polymerase Chain Reaction (PCR) and
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 , 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
cetyltrimethylammonium bromide (CTAB) protocol by Doyle and Doyle
 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., '?’, ) in downstream
We aligned our DNA sequences using ClustalX  and
then by eye in MEGA4.0  following the guidelines of
Morrison . 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 . 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  with
the GTR+ G + I and GTR+ G models of nucleotide
Table 2 Materials and GenBank accession numbers of five chloroplast makers and accession information
Genbank accession numbers (bold indicated contributed by this study)
matK rbcL trnL-trnF rpl32-trnL psbA-trnH
Table 3 Primers and sequences statistics of nuclear and chloroplast makers used in present study
ACTGCCTTGATCCACTTGGC CGAAGCTCCATCTACAAATGG   
substitution, respectively. These models were selected
under the Akaike information criterion (AIC) using
MrModeltest version 2.2 . 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  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* . 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
KishinoHasegawa (KH) tests in PAUP*, because KH tests are
known to exhibit very low type I error rates . 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  (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
LiliumNomocharis has been postulated [13, 17, 57] and
incomplete lineage sorting has been detected in many plant
lineages . Therefore, we used the statistical
parsimony network approach implemented in TCS v.1.21
 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
Molecular dating in Liliales has been previously
performed using distantly related fossils , calibrations
from previous studies [44, 61], and single calibration
points . In particular, Bremer  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 ,
Patterson and Givnish  inferred the divergence time of
the tribe Lilieae as 12 Mya and Vinnersten and Bremer
 concluded that the monophyletic lineage comprised
of Lilium, Nomocharis and Fritillaria diverged 6 ± 2.9 Ma.
Gao et al.  provided a detailed review of Liliales fossils
and performed dating using a single, reliable fossil of
Smilax, Smilax wilcoxensis Berry , 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 [63–65].
We sought to more rigorously date events in
LiliumNomocharis 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
LiliaceaeSmilacaceae 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 . For the second
calibration, we used Ripogonum tasmanicum Conran, et al.
 to constrain the age of the ancestor of the monotypic
Ripogonanceae and Philesiaceae (following Angiosperm
Phylogeny Website, ). Ripogonum tasmanicum is
reported from the Tasmanian Macquarie Harbour Formation
, which is approximately 51–52 million years old based
on a foraminiferal index . 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
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 . 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  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 
rejected strict clocks for both datasets (P < 0.01), therefore
we used an uncorrelated lognormal (UCLN), relaxed
clock . 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 (http://www.phylo.org/portal2).
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 .
Ancestral Area Reconstructions (AAR)
We used the Bayesian Binary method (BBM) in
Reconstruct Ancestral States in Phylogenies 2.1b (RASP 2.0)
[77–79] 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
SinoJapanese 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-ExtinctionCladogenesis (DEC) model . 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
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 . For elevation, we
acquired information from floras and specimen records on
GBIF (http://www.gbif.org/). We treated elevation as
categorical by using 1000 ft. increments for our discrete
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  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.
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  and Liang . 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.
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
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
Additional file 11: Table S4. Morphological character states used in
ancestral state reconstruction.
The authors declare that they have no competing interests.
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
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
1. Myers N , Mittermeier RA , Mittermeier CG , da Fonseca GAB , Kent J. Biodiversity hotspots for conservation priorities . Nature . 2000 ; 403 : 853 - 8 .
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 .
3. Yue JP , Sun H , Al-Shehbaz IA , Li JH . Support for an expanded SlomsLaubachia (Brassicaceae): evidence from sequence of chloroplast and nuclear genes . Ann Mo Bot Gard . 2006 ; 93 : 402 - 11 .
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 .
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 .
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 .
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 .
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 .
9. Liu JQ , Tian B. Origin , evolution, and systematics of Himalaya endemic genera . Newslett Himalayan Botany . 2007 ; 40 : 20 - 7 .
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 .
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 .
12. Balfour B. The Genus Nomocharis . Bot J Scotl . 1918 ; 27 ( 3 ): 273 - 300 .
13. Sealy JR . Nomocharis and Lilium. Kew Bull . 1950 ; 5 ( 2 ): 273 - 97 .
14. Sealy JR . A revision of the genus Nomocharis Franchet . Bot J Linn Soc . 1983 ; 87 : 285 - 323 .
15. Liang SY . Studies on the genus Nomocharis (Liliaceae). Bull Botanical Res . 1984 ; 4 ( 3 ): 163 - 78 .
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 .
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 .
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 .
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 .
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 .
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 .
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 .
23. Comber HF . A new classification of genus Lilium . Royal Horticult Soc Liliy Year Book . 1949 ; 13 : 85 .
24. Haw SG . The Lilies of China . Portland: Timber Press; 1986 .
25. Chen FB . H-D event: an important tectonic event of the late Cenozoic in Eastern Asia . Mt Res . 1992 ; 10 ( 4 ): 195 - 202 .
26. Chen FB . Second discussion on the H-D movement. Mt Res . 1996 ; 17 ( 3-4 ): 14 - 22 .
27. Harrison TM , Copeland P , Kidd WSF , Yin A . Raising Tibet. Science . 1992 ; 255 : 1663 - 70 .
28. Li JJ , Shi YF , Li BY . Uplift of the Qinghai-Xizang (Tibet) Plateau and Global Change . Lanzhou, China: Lanzhou University Press ; 1995 .
29. Chapman MA , Hiscock SJ , Filatov DA . Genomic divergence during speciation driven by adaptation to altitude . Mol Biol Evol . 2013 ; 30 ( 12 ): 2553 - 67 .
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 .
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 .
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 .
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 .
34. Köener C. Alpine plant life: functional plant ecology of high mountain ecosystems . Berlin: Springer; 2003 .
35. Ushimaru A , Dohzono I , Takami Y , Hyodo F. Flower orientation enhances pollen transfer in bilaterally symmetrical flowers . Oecologia . 2009 ; 160 : 667 - 74 .
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 .
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 .
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 .
39. Fenster CB , Armbruster WS , Dudash MR . Specialization of flowers: is floral orientation an overlooked step? New Phytol . 2009 ; 183 ( 3 ): 502 - 6 .
40. Maddison W , Knowles LL . Inferring phylogeny despite incomplete lineage sorting . Syst Biol . 2006 ; 55 ( 1 ): 21 - 30 .
41. Richardson AO , Palmer JD . Horizontal gene transfer in plants . J Exp Bot . 2007 ; 58 ( 1 ): 1 - 9 .
42. Minami S , Azuma A. Various flying modes of wind-dispersal seeds . J Theor Biol . 2003 ; 225 ( 1 ): 1 - 14 .
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 .
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 .
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 .
46. Doyle JJ , Doyle JL. A rapid DNA isolation procedure for small quantities of fresh leaf tissue . Phytochem Bull . 1987 ; 19 : 11 - 5 .
47. Wiens JJ , Morrill MC . Missing data in phylogenetic analysis: reconciling results from simulations and empirical data . Syst Biol . 2011 ; 60 ( 5 ): 719 - 31 .
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 .
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 .
50. Morrison DA . A framework for phylogenetic sequence alignment . Plant Syst Evol . 2009 ; 282 : 127 - 49 .
51. Swofford DL . PAUP* . Phylogenetic Analysis Using Parsimony (*and Other Methods) . 4th ed. Sunderland, Massachusetts: Sinauer Associates; 2003 .
52. Ronquist F , Huelsenbeck JP . MrBayes 3: Bayesian phylogenetic inference under mixed models . Bioinformatics . 2003 ; 19 ( 12 ): 1572 - 4 .
53. Nylander JAA . MrModeltest 2.0. In: Department of Systematic Zoology. 20th ed . Uppsala: EBC, Uppsala University; 2004 .
54. Rambaut A , Drummond AJ . Tracer v1 .4. 2007 . http://beast.bio.ed.ac.uk/Tracer.
55. Susko E. Tests for two trees using likelihood methods . Molecular Biology and Evolution . 2014 ; 31 ( 4 ): 1029 - 39 .
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 .
59. Crandall KA , Templeton AR . Empirical tests of some predictions from coalescent theory with applications to intraspecific phylogeny reconstruction . Genetics . 1993 ; 134 : 959 - 69 .
60. Bremer K. Early Cretaceous lineages of monocot flowering plants . Proc Natl Acad Sci U S A . 2000 ; 97 ( 9 ): 4707 - 11 .
61. Vinnersten A , Bremer K. Age and biogeography of major clades in Liliales . Am J Bot . 2001 ; 88 : 1695 - 703 .
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 .
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 .
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 .
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 .
69. Potter F , Dilcher D. Biostratigraphic analysis of Eocene clay deposits in Henry County , Tennessee. Biostratigraph Fossil Plants . 1980 ; 211 : 225 .
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 .
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 .
73. Stevens PF . Angiosperm Phylogeny Website. Version 12 [http://www.mobot.org/MOBOT/research/APweb/]
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 .
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 .
76. Drummond AJ , Ho SYW , Phillips MJ , Rambaut A. Relaxed phylogenetics and dating with confidence . PLoS Biol . 2006 ; 4 ( 5 ): e88 .
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 .
78. Yu Y , Harris A , He XJ . RASP (Reconstruct Ancestral State in Phylogenies) 2.0 beta . In: http://www.mnh.scu.edu.cn/soft/blog/RASP; 2012 .
79. Yu Y , Harris AJ , He XJ . A novel Bayesian method for reconstructing geographic ranges and ancestral states on phylogenies . 2011 .
80. McRae EA . Lilies: a guide for growers and collectors . Portland: Timber Press ; 1998 .
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 .
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 .
83. Stöver BC , Müller KF. TreeGraph 2: combining and visualizing evidence from different phylogenetic analyses . BMC Bioinformatics . 2010 ; 11 :7.
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 .
85. Levan A. Cytological studies in Allium, II Chromosome morphological contributions . Hereditas . 1932 ; 16 : 257 - 94 .
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 .
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 .
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 .
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 .
91. Yue JP , Gu ZJ , Al-Shehbaz IA , Sun H. Cytological studies on the SinoHimalayan endemic Solms-laubachia (Brassicaceae) and two related genera . Bot J Linn Soc . 2004 ; 145 ( 1 ): 77 - 86 .
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 .
93. Fay MF , Swensen SM , Chase MW . Taxonomic affinities of Medusagyne oppositifolia (Medusagynaceae) . Kew Bull . 1997 ; 52 ( 1 ): 111 - 20 .
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 .
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 .
96. Hamilton MB . Four primer pairs for the amplification of chloroplast intergenic regions with intraspecific variation . Mol Ecol . 1999 ; 8 ( 3 ): 521 - 3 .