\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Effects of Mountain Uplift and Climatic Oscillations on Phylogeography
and Species Divergence of Notholirion (Liliaceae)}



\author[1]{Rui-Yu Cheng}%
\author[1]{Juan Li}%
\author[1]{Deng-Feng Xie}%
\author[1]{renxiu zhou}%
\author[2]{Qing Li}%
\author[2]{Yanglina Yu}%
\author[1]{Xingjin He}%
\author[1]{Songdong Zhou}%
\affil[1]{Sichuan University}%
\affil[2]{Chengdu Branch of Giant Panda National Park}%


\vspace{-1em}



  
  \date{May 21, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Exploring the geological events and climate change in the
Himalaya-Hengduan Mountains (HHM) region and the Qinghai-Tibet Plateau
(QTP) is crucial for understanding the impact of environmental change on
biogeographic distribution and biological evolution. To delve deeper
into these mechanisms, we reconstructed the evolutionary history of
three Notholirion species that span these regions. Here, we examined a
total of 254 individuals from 31 populations of these three species,
utilizing five chloroplast DNA (cpDNA) (matK, ndhA, ndhG-ndhI,
petB-petD, and petL-petG), along with one nuclear DNA region (internal
transcribed spacer, ITS). We identified 14 haplotypes from cpDNA and 27
haplotypes from the ITS, each specific to corresponding species. Robust
haplotype trees were detected, and obvious discrepancies were found
between the cpDNA and ITS trees. A total of 147 chloroplast genomes were
used for divergence time estimation, of which 10 chloroplast genomes
from distinct populations of Notholirion species provided a
comprehensive representation of the genus. The divergence time
estimation results suggested that species of Notholirion genus
originated in the southern Himalayan region during the Late Oligocene
period (25.05 Ma), and the three Notholirion species then diverged
during the Late Pliocene period (7.43 Ma). Our maximum model forecasts
that the overall distribution range of Notholirion over four different
periods remains relatively stable,from LIG to the future. The origin of
the genus Notholirion was triggered by sustained climate fluctuations
during the Late Oligocene, with the uplift of the Himalayas and the
subsequent orogenic movements intensifying climate changes and further
promoting the species divergence of Notholirion. A long evolutionary
history, coupled with sexual reproduction and habitat fragmentation,
likely contributed to the higher genetic diversity of Notholirion. The
higher genetic differentiation among Notholirion populations may be
attributed to drastic changes in the external environment within their
range, as well as their constrained capacity for seed production and
dispersal.%
\end{abstract}%



\sloppy


\textbf{1 Introduction}

The Qinghai-Tibet Plateau (QTP), as the ''roof of the world'', has
fostered a variety of plateau organisms, and the Himalayan-Hengduan
Mountains (HHM) on \emph{ITS} margin represent two distinct biodiversity
hotspots (Deng, Ding, \& Deng, 2015; Favre et al., 2015; H. H. Meng et
al., 2017). As key features of the biodiversity hotspots of East and
South Asia, Himalaya Mountains define the southern edge of the QTP,
whereas the Hengduan Mountains of Southwest China form the southeastern
frontier of the plateau(Marchese, 2015; Y. L. Zhang, Li, \& Zheng,
2002). Both regions have impressive animal and plant diversity, with
alpine plant diversity being an important contributor to the hotspot (D.
Z. Li, 2008; J. Li et al., 2020; C. Xie et al., 2019). Previous studies
have suggested that the Himalayas were formed by the collision of the
Indo-Australian plate with Eurasian plate, accompanied by continuous and
prolonged uplifts(Aitchison \& Ali, 2008; Dupont-Nivet, Lippert,
Hinsbergen, Meijers, \& Kapp, 2010; H. B. Zheng, Powell, An, Zhou, \&
Dong, 2000). The collision of the Indian plate and Eurasia resulted in
the consequent rise of the Himalayas and QTP. This collision also
created folded mountains and faulted basins that formed the present-day
HHM(Chatterjee, Goswami, \& Scotese, 2013; Dupont-Nivet et al., 2010).
Most studies suggest that the collision began in the Eocene, around
55-50 Mya, followed by a continuous uplift of the Himalayas. parts of
the region, such as Thakkhola, Gyirong and Zhada, reached an average
elevation of about 4000 to 6000 m during the Miocene (Chatterjee et al.,
2013; Dupont-Nivet et al., 2010; Saylor et al., 2009; Q. Zhang \&
Willems, 2012). During the Late Miocene and the Pliocene, the QTP
experienced further uplift. The uplift of the HHM on the southeastern
edge of the QTP was particularly pronounced and reached \emph{ITS} peak
elevation rapidly before the Late Pliocene(Mulch \& Chamberlain, 2006;
Xiong et al., 2022; H. B. Zheng et al., 2000).

The formation and uplift of the QTP have had a profound impact on the
local area and \emph{ITS} surrounding regions. This has resulted in a
series of topographic and climatic changes, which has especially shaped
the unusually rugged and complex topography of the HHM and the special
East Asian monsoon climate of this region(Cao, Shi, Zhang, \& Wang,
2009; Huntington, Blythe, \& Hodges, 2006; X. Sun \& Wang, 2005; P. Z.
Zhang et al., 2004). Since the Quaternary period, climatic fluctuations
have had significant effects on the distribution and history of plants
populations in the Northern Hemisphere, causing their range-wide
migration or extinction and possibly driving local adaptation(J. M.
Chen, Zhi-Yuan, Yuan, \& Wang, 2014; G. M. Hewitt, 2004; Hickerson et
al., 2010). Previous studies have shown that there were four major ice
ages and three interglacials that occurred during the Pleistocene, with
each ice period being separated by interglacial period. The repeated
alternation of ice ages and interglacials caused oscillations of warm
and cold climates, which profoundly affected the genetic diversity and
spatial distribution patterns of biological taxa found on the QTP and
adjacent areas(Borns, 1994; G. Hewitt, 2000; Keller \& Seehausen, 2012;
Z. Zhang \& Sun, 2011). And the repeated alternating cycles of
population contraction during the ice age and expansion after the ice
age became the dominant pattern for many plants in the Northern
Hemisphere, enabling their survival throughout the Quaternary ice
age(Godfrey M. Hewitt, 1996). On the QTP and \emph{ITS} adjacent areas,
the HHM region in the southeastern part of the QTP may have served as an
important refuge for endemic species during the ice age, thanks
to\emph{ITS} distinctive mountainous topography(Q. B. Gao et al., 2012;
J. Q. Liu, Duan, Hao, XuemgUn, \& Sun, 2014; Qiu, Fu, \& Comes, 2011).
This region is characterized by north-south peaks separated by valley
bottoms(Song et al., 2016). From valleys situated approximately 1,000 m
above sea level to the peakof Minya Konka at 7,556 m, the huge drop in
altitude can form ``sky islands''(Kai \& Jiang, 2014; Qian, Ricklefs, \&
Thuiller, 2021). During the Quaternary, these low elevation valley areas
on the edge of the ice age impact may have provided shelter for certain
local plants to survive the ice age(Q. B. Gao et al., 2012; J. Q. Liu et
al., 2014). In addition, many studies have suggested that during the Ice
Age, plant species from the QTP region may have survived in invisible
refuges at high altitudes(J. Q. Liu, Sun, Ge, Gao, \& Qiu, 2012; Luo,
Yue, \& Sun, 2016). These different habitat zones and ecosystems, such
as the refuges found in different elevation zones and sky islands, have
also resulted in speciation and radiation(J. Q. Liu et al., 2014; Qiu et
al., 2011; H. Wang et al., 2010; Xu et al., 2010). Thus, the climatic
and specific topographic changes caused by the uplift of the QTP have
played a key role in the origin, speciation and evolution of several
plant taxa throughout the HHM and adjacent areas.

The HHM region has a special climatic environment, complex geological
landscape, and bioecological diversity, making it an ideal and
captivating location for studying species diversification and adaptive
evolution. Previous studies have mostly focused on how species in the
HHM adapt to climatic oscillations brought by uplift of the plateau(Luo
et al., 2016; Shahzad, Jia, Chen, Zeb, \& Li, 2017; Xu et al., 2010; J.
Yang, Zhou, Huang, \& He, 2018). The refuge pattern of species in the
QTP region during the ice age, as well as their migration out of the
plateau, is also a topic of current research (Q. Li et al., 2020; H. R.
Liu et al., 2022; M. L. Liu, He, L\selectlanguage{ngerman}ópez-Pujol, Jia, \& Li, 2019; Shahzad
et al., 2017). Available phylogenetic and biogeographic studies have
suggested that population differentiation and demography within species
have also been strongly influenced by the uplift of the QTP and the
consequent climate oscillations during the Quaternary ice age in the
region(H. R. Liu et al., 2022; Y. Sun et al., 2014; H. Wang et al.,
2010; F. S. Yang, Qin, Li, \& Wang, 2012). In recent years, while most
studies continue to focus on population differentiation and biogeography
within species, very few researchers have focused on an entire endemic
genus in the HHM region(H. Y. Zheng, Guo, Price, He, \& Zhou, 2021).
Therefore, we expect to investigate the origins, diversification, and
population response of species to climatic oscillations since the
Pleistocene by studying the distribution and population history of plant
taxa that are mainly distributed in the HHM and QTP region.

\emph{Notholirion} is a genus of Liliaceae and is mainly distributed in
the HHM and QTP region, where it grows in alpine meadows and
shrublands(M. Yang, Zhou, Xingjin, \& Peng, 2016).\emph{Notholirion}
species are distinguished from other genera of Liliaceae by the presence
of racemes and numerous small bulbs surrounding their bulbs. According
to the records of \emph{Notholirion}in the Flora of China and Li et al.
(2021), four morphologically distinct species are accepted, and three
species occur in the HHM and QTP region, except for \emph{N. koeiei} ,
which are found in Iran and Iraq in the West Asia region(J. Li et al.,
2022; J. Q. Liu et al., 2014). Particularly, the two species with
distribution in China (\emph{N. macrophyllum} and \emph{N. bulbuliferum}
) both live at altitudes above 2800 m and are typical plants of the
QTP(M. Yang et al., 2016). \emph{N. bulbuliferum} exhib\emph{ITS} the
widest distribution in China, ranging from the Qinling Mountains in
Shaanxi to Western Sichuan, northwestern Yunnan, and eastern Tibet.
\emph{N. macrophyllum} is restricted to specific areas such as Jilong in
Tibet and Daocheng in Sichuan(J. Li et al., 2022). Additionally, we were
fortunate to acquire a population of \emph{N. thomsonianum} from Nepal.
\emph{N. thomsonianum} does not occur naturally in China, but to
maximize the representation of \emph{Notholirion} species, we included
the only available population from Nepal in our analysis. As for
\emph{N. koeiei} , there are only five recorded specimens from Iraq and
Iran in western Asia, and unfortunately, we did not have access to any
samples of this species. Based on the descriptions in the Flora of China
and our own field observations, we have noted several distinct
characteristics among the \emph{Notholirion} species. Mature plants of
\emph{N. macrophyllum} typically have a height of less than 40 cm, and
their racemes typically bear 3-7 flowers. On the other hand, \emph{N.
bulbuliferum} and \emph{N. thomsonianum} generally grow taller than 40
cm, and their racemes typically bear a larger number of flowers.
Furthermore, \emph{N. bulbuliferum} lacks any stripes at the base of the
tepals and typically flowers from July to August, whereas \emph{N.
thomsonianum} exhib\emph{ITS} longitudinal dark stripes consisting of
spots at the base of the tepals and flowers from March to May(S. C.
Chen, Liang, Xu, \& Minoru, 2000; J. Li et al., 2022). The recent
phylogenetic study indicated that \emph{N. bulbuliferum} and \emph{N.
macrophyllum} are sister taxa, and they form a sister clade together
with \emph{N. thomsonianum} . The three \emph{Notholirion} species form
a distinct monophyletic clade. Therefore, this genus is an ideal taxon
to investigate the evolution and diversification of alpine plants in the
HHM and QTP region(J. Li et al., 2022).

In this study, we reconstructed the phylogenetic tree
of\emph{Notholirion} based on data from 31 populations of three species
that distributed in the HHM and QTP regions. We used biogeographical and
population genetics approaches to identify the mechanism of origin,
genetic structure, and lineage evolutionary history in the
genus\emph{Notholirion} . Our aims were to: (i) reconstruct the
phylogenetic relationships within the genus \emph{Notholirion} using
materials from 31 populations; (ii) investigate the origin and species
diversification of \emph{Notholirion} in the HHM and QTP regions; (iii)
explore the geological and climatic changes that influenced distribution
shifts and genetic diversity of \emph{Notholirion} . We hope that this
study will provide some clues to explore the evolutionary patterns of
species in the HHM region.

\textbf{2 Materials and Methods}

\textbf{2.1 Sample collection}

To observe as much variation of the species as possible in the wild, a
total of 254 individuals belonging to three species of\emph{Notholirion}
were collected based on our field explorations and past herbarium
records. These individuals were from 31 populations at elevations above
2800 m in Nepal, Tibet (Xizang), Sichuan and Yunnan provinces
\textbf{(Figure 1 andSupplementary Table 1)} . Sampling localities
covered almost the entire distribution of the three species and 6 to 10
individuals were randomly sampled from each population with at least 20
m between the collected samples. Fresh leaves were collected and dried
immediately with silica gel. The voucher specimens were deposited at
Sichuan University Herbarium (SZ).

\textbf{2.2 DNA Extraction, PCR Amplification, and Sequencing}

Total genomic DNA was extracted from dried leaf tissue using a modified
cetyltrimethylammonium bromide (CTAB) protocol(E. \& C., 1980). Raw data
of the 10 newly sequenced whole plastid genomes from
\emph{Notholirion}species were generated by Illumina platform,
generating 150~bp paired-end reads at Novogene (Beijing, China). The
reads were assembled using the program NOVOPlasty v2.6.2 (Dierckxsens,
Mardulyn, \& Smits, 2017) and the genome annotation was performed using
Plastid Genome Annotator (PGA)(Qu, Moore, Li, \& Yi, 2019). Manual
adjustment compared with related species' plastomes was conducted in
Geneious v9.0.2(Kearse et al., 2012). For the plastid DNA (cpDNA)
fragment, we designed 15 primer pairs using the whole plastid genome of
\emph{Notholirion}(GenBank accessions NC046464, MH011354) and finally
got five chloroplast fragments (\emph{matK} , \emph{ndhA} ,
\emph{ndhG-I} , \emph{petB-D} and\emph{petL-G} ) for analysis. We
amplified and sequenced \emph{ITS}(\emph{ITS1--5.8sRNA--ITS2} )(White,
Bruns, Lee, \& Taylor, 1990) and the five chloroplast fragments above at
Sangon Biotech Co. Ltd. (Chengdu, China) and all DNA fragments were
aligned and manually adjusted using the software MEGA X(Kumar, Stecher,
Li, Knyaz, \& Tamura, 2018).

\textbf{2.3 Phylogenetic analyses}

The phylogeny reconstruction based on two databases of haplotypes
(\emph{ITS} and cpDNA sequences) were performed using Bayesian inference
(BI) and maximum likelihood (ML). Gaps were treated as missing data
during the tree searches. The two programs MrBayes version 3.2(Ronquist
et al., 2012) and RAxML v.8(Stamatakis, 2014) were used to conduct BI
and ML analyses. The best model for \emph{ITS} and cpDNA were both GTR+G
(BI and ML) inferred by ModelFinder(Kalyaanamoorthy, Minh, Wong, von
Haeseler, \& Jermiin, 2017) based on the Akaike information criterion
(AIC). Markov chain Monte Carlo (MCMC) method was used for BI and set to
run 10 million generations, sample every 1,000 generations, discarde
25\% of the trees as burn-in. The convergence of MCMC inference was
evaluated using Tracer v.1.7.11 to ensure sufficient samples (ESS
\textgreater{} 200) by checking the ESS (Effective Sample Size) value.
Besides, we ran RAxML on\emph{ITS} and cpDNA with 1000 bootstrap
replicates. Finally, BI posterior probabilities (BI PP) and ML bootstrap
support (ML BS) were presented at the nodes separately.

\textbf{2.4 Divergence Time Estimate}

There are currently no well-documented fossils in Liliaceae, and thus
fossil constraints were limited to Liliales. We used 10 populations
of\emph{Notholirion} plastomes and 137 other plastomes to calibrate the
divergence time of \emph{Notholirion} \textbf{(Supplementary Table 2)} .
Referring to previous studies(H. T. Li et al., 2019; J. Li et al., 2022;
D. F. Xie et al., 2020) and to minimize the effects of missing data, we
only used combined single-copy CDS genes data set derived from 147
plastomes for the estimation of divergence time. The common CDS of 147
plastomes were extracted, respectively aligned, and then concatenated as
the plastome CDS dataset, using PhyloSuite v1.2.2(D. Zhang et al.,
2020). Estimations of divergence time were performed using an
uncorrelated lognormal relaxed molecular clock method implemented in the
BEAST 1.10.4 program(Suchard et al., 2018). BEAUti was used to set
criteria for analyses where we used a GTR + G substitution model
selected by ModelFinder and a Yule process(Gernhard, 2008). An optimal
partitioning scheme was determined using the PartitionFinder 2
program(Lanfear, Frandsen, Wright, Senfeld, \& Calcott, 2017) and six
calibration points were used to calibrate time \textbf{(Supplementary
Table 3)} . MCMC analyses were run for 100 million generations with
parameters sampled every 10,000 generations after discarding the first
10\% of generations as burnin. The convergence of the stationary
distribution was accessed by ESS values (\textgreater{}200) using the
Tracer v.1.7.11. Maximum clade credibility (MCC) trees were produced
using TreeAnnotator v1.8.4(Helfrich, Rieb, Abrami, Lücking, \& Mehler,
2018) and visualized in FigTree v1.4.2(Rambaut, 2014).

\textbf{2.5 Ancestral Area Reconstruction}

To explore the biogeographic history of \emph{Notholirion} plants, we
used three regions based on the paleogeographic and climatic evidence,
geographic distribution and floral composition throughout the
Himalayas-Hengduan Mountains (HHM) for ancestral area
reconstruction(Buerki et al., 2012; Tu, Volis, Dillon, Sun, \& Wen,
2010; H. Y. Zheng et al., 2021): (A) southern Himalayas, (B) eastern
Himalayas, (C) Hengduan Mountains. The Statistical Dispersal-Vicariance
(S-DIVA) analysis as implemented in the RASP v4(Y. Yu, Blair, \& He,
2020) was used for the reconstruction. We used the BI tree based on the
cpDNA datasets (obtained by phylogenetic analyses) for the S-DIVA
analyses. The uncertainty in the root areas of an outgroup can lead to
biased inferences for the crown node of the ingroups(Y. Yu, Harris,
Blair, \& He, 2015), outgroups were removed before ancestral-state
reconstruction. To explore the effects of area constraints, the maximum
number of areas at each node was set to two.

\textbf{2.6 Population Genetics, Phylogeographic Analyses, and
Demographic History}

All the DNA sequences were edited by SeqMan (DNAstar package; DNAStar
Inc., Madison, WI, United States) to obtain consensus sequences. The
program MAFFT v.7.369b(Katoh \& Standley, 2013) was used to align for
subsequent manual adjustments. Haplotype diversity
(H\textsubscript{d})(Nei \& Tajima, 1981) and nucleotide diversity
(\selectlanguage{greek}π\selectlanguage{english})(Nei \& Li, 1979) for each population were calculated using DNAsp
v6.0(Rozas et al., 2017) to verify the degrees and patterns of
diversity. The haplotypes of the two gene fragment datasets and the
distribution of all haplotypes on the map were then plotted separately
in the software PopART using the Minimum spanning method(Bandelt,
Forster, \& R\selectlanguage{ngerman}öhl, 1999). PERMUT was used to access the total diversity
(H\textsubscript{T}), within-population diversity (H\textsubscript{S})
and population differentiation indices (G\textsubscript{ST} and
N\textsubscript{ST})(Pons \& Petit, 1996). And a U-statistic was used to
test the phylogeographic structure by comparing G\textsubscript{ST} and
N\textsubscript{ST}. In addition, analyses of molecular variance (AMOVA)
was performed using ARLEQUIN v3.5(Excoffier \& Lischer, 2010) to assess
the genetic differentiation within and between populations, and 1000
random permutations were conducted to test the significance of
partitioning. To detect whether the populations in
\emph{Notholirion}experienced a recent population expansion, a mismatch
distribution analysis (MDA)(Schneider \& Excoffier, 1999) was carried
out using the ARLEQUIN v3.5(Excoffier \& Lischer, 2010). Additionally,
Neutrality test (Fu's Fs; Tajima's D)(Tajima, 1989) was also conducted
to test whether there was potential population expansion in
\emph{Notholirion} using DnaSP v6.0 program. The smoothness of observed
mismatch distribution was detected by calculating the Harpending's
raggedness index (Hrag)(Schneider \& Excoffier, 1999) and the sum of
squared deviations (SSD) between observed and expected mismatch
distributions.

\textbf{2.7 Species Distribution Modeling}

MAXENT 3.3.3e(Phillips, Anderson, \& Schapire, 2006) was used to predict
the distribution of \emph{Notholirion} in China during 4 time periods:
Last Interglacial (LIG), Last Glacial Maximum (LGM), and the present and
future. The samples collected in this study covered the distribution of
all 30 populations of the genus \emph{Notholirion} throughout China for
which specimens were recorded. Nineteen bioclimatic environment
variables(Alemu \& Wimberly, 2020) from four time periods above were
downloaded from the WorldClim database(Fick \& Hijmans, 2017) and
employed in modeling analyses. In MAXENT, due to the small sample size
of our data (15--79 localities), only the linear + quadratic + hinge
functional forms were used. Here, models were built using linear +
quadratic + hinge features, random training data was set to 25\%, repeat
run type was set to Subsample, maximum iterations were set to 1000, and
other parameters were set by default. We also calculated the area under
the Receiver Operating Characteristic (ROC)(Fawcett, 2005) and Operating
Characteric Curve (AUC)(Ellis, Yahr, Belinchón, \& Coppins, 2014) to
observe the accuracy of each model prediction. Usually, good model
performance was assessed by AUC values above 0.7(Elith \& Leathwick,
2009; Fielding, 1997; Peterson, Pape, \& Soberón., 2008). The spatial
distribution range map of the final modeling results was plotted by the
software DIVA-GIS v 7.5(Hijmans, Cameron, Parra, Jones, \& Jarvis,
2010).

\textbf{3 Results}

\textbf{3.1 Genetic Diversity and Structure}

\emph{ITS} and five chloroplast fragments (\emph{matK} , \emph{ndhA}
,\emph{ndhG-I} , \emph{petB-D} , \emph{petL-G} ) were used to analyze
254 individuals from 31 populations of the nine \emph{Notholirion}
species. The \emph{ITS} and cpDNA haplotype frequencies of each
population are listed in \textbf{Supplementary Table 4} , and the
geographical distributions of haplotypes are shown in \textbf{Figures 2,
4} . The length of \emph{ITS} sequence we used in this study was 671 bp
containing 60 polymorphic sites. We detected 27 \emph{ITS} haplotypes
(N1--N27), which were species-specific as none were shared by any two
species \textbf{(Figure 3)} . Most of the populations have only one or
two haplotypes and populations with more than three haplotypes only
occurred in \emph{N. bulbuliferum} . Besides, the total length of the
aligned sequences of cpDNA was 2753 bp containing 39 polymorphic sites
were recovered with 14 haplotypes. Like \emph{ITS} haplotypes, all cpDNA
haplotypes were spec 09ies-specific in \emph{Notholirion}\textbf{(Figure
5)} . All the populations have only one or two haplotypes, with only
22.58\% of the populations having two haplotypes.

The haplotype diversity (Hd) of \emph{ITS} ranged from 0.000 to 0.867,
and nucleotide diversity (\selectlanguage{greek}π\selectlanguage{english}) ranged from 0.000 to 0.011 across all
species \textbf{(Supplementary Table 4)} . Total gene diversity
(H\textsubscript{T}) value was all higher than the average gene
diversity within populations (H\textsubscript{S}) at the level of
species and genus. Additionally, the number of substitution types
(N\textsubscript{ST}) was higher than interpopulation differentiation
(G\textsubscript{ST}), which indicated that a significant
phylogeographic structure existed in \emph{Notholirion} (Table 1). The
haplotype diversity (H\textsubscript{d}) of cpDNA ranged from 0.000 to
0.800, and nucleotide diversity (\selectlanguage{greek}π\selectlanguage{english}) ranged from 0.000 to 0.004 across
all species \textbf{(Supplementary Table 4)} . Like
\emph{ITS}haplotypes, the total gene diversity (H\textsubscript{T})
value was higher than the average gene diversity within populations
(H\textsubscript{S}), whether at the level of species or genus. Total
NST was significantly higher than GST, indicating significant
phylogeographic structure in \emph{Notholirion} \textbf{(Table 1)} .

\emph{ITS} AMOVA revealed that the primary genetic variation (74.86\% -
98.68\%) in \emph{N. bulbuliferum} and \emph{N. macrophyllum} occurred
among populations (\textbf{Table 2} ). Furthermore, the AMOVA analysis
of cpDNA data indicated a genetic differentiation pattern similar to
that of \emph{ITS} , with 93.73\% - 98.85\% of the total variation
distributed among populations (\textbf{Table 2} ).

\textbf{Table} \textbf{1 Genetic diversity and genetic differentiation
of three species within}

\emph{Notholirion} \textbf{based on cpDNA and \emph{ITS}.}\selectlanguage{english}
\begin{longtable}[]{@{}lllll@{}}
\toprule
\textbf{Species} & \textbf{H\textsubscript{S}} &
\textbf{H\textsubscript{T}} & \textbf{G\textsubscript{ST}} &
\textbf{N\textsubscript{ST}}\tabularnewline
\midrule
\endhead
\textbf{nrDNA} & & & &\tabularnewline
\emph{N. bulbuliferum} & 0.199 (0.0611) & 0.750 (0.0711) & 0.705
(0.0741) & 0.716 (0.0781)\tabularnewline
\emph{N. macrophyllum} & 0.156 (0.1556) & 0.667 (0.3156) & 0.767 (NC) &
0.986 (NC)\tabularnewline
\emph{N. thomsonianum} & - & - & - & -\tabularnewline
Total & 0.227 (0.0617) & 0.796 (0.0623) & 0.715 (0.0705) & 0.767
(0.0899)\tabularnewline
\textbf{cpDNA} & & & &\tabularnewline
\emph{N. bulbuliferum} & 0.107 (0.0395) & 0.766 (0.0529) & 0.860
(0.0476) & 0.907 (0.0369)\tabularnewline
\emph{N. macrophyllum} & 0.119 (0.1185) & 0.667 (0.2933) & 0.822 (NC) &
0.983 (NC)\tabularnewline
\emph{N. thomsonianum} & - & - & - & -\tabularnewline
Total & 0.677 (0.0600) & 0.931 (0.0080) & 0.273 (0.0640) & 0.402
(0.1710)\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table2 Analysis of molecular variance (AMOVA) based on the cpDNA
and nrDNA}\selectlanguage{english}
\begin{longtable}[]{@{}lllllll@{}}
\toprule
\begin{minipage}[b]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[b]{0.12\columnwidth}\raggedright\strut
Source of variation\strut
\end{minipage} & \begin{minipage}[b]{0.12\columnwidth}\raggedright\strut
d.f.\strut
\end{minipage} & \begin{minipage}[b]{0.12\columnwidth}\raggedright\strut
SS\strut
\end{minipage} & \begin{minipage}[b]{0.12\columnwidth}\raggedright\strut
VC\strut
\end{minipage} & \begin{minipage}[b]{0.12\columnwidth}\raggedright\strut
PV (\%)\strut
\end{minipage} & \begin{minipage}[b]{0.12\columnwidth}\raggedright\strut
Fixation Indices\strut
\end{minipage}\tabularnewline
\midrule
\endhead
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{ITS}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{ITS}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{ITS}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{ITS}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{ITS}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{ITS}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{ITS}\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
All species\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
Among species Among populations within species Within populations\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
1\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
195.447\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
5.12027\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
70.37\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
FSC: 0.79869\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
27\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
383.922\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
1.72216\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
23.67\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
FST: 0.94035\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
205\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
88.986\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
0.43408\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
5.97\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
FCT: 0.70367\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\emph{N.} \emph{bulbuliferum}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Among populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
26\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
296.338\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
1.37585\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
74.86\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
FST: 0.74856\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Within populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
188\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
86.886\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
0.46216\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
25.14\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\emph{N. macrophyllum}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Among populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
1\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
87.584\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
9.23196\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
98.68\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
FST: 0.98680\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Within populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
17\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
2.100\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
0.12353\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
1.32\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\emph{N. thomsonianum}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Among populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\emph{Within populations}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\emph{-}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\emph{-}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\emph{-}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\emph{-}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\emph{-}\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
All samples\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Among populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
30\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
1060.215\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
4.27072\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
91.45\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
FST: 0.91455\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Within populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
223\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
88.986\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
0.39904\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
8.55\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{cpDNA}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{cpDNA}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{cpDNA}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{cpDNA}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{cpDNA}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{cpDNA}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\textbf{cpDNA}\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
All species\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
Among species Among populations within species Within populations\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
1\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
166.421\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
3.18145\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
76.19\strut
\end{minipage} & \begin{minipage}[t]{0.14\columnwidth}\raggedright\strut
FSC: 0.95113\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
28\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
214.406\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
0.94567\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
22.65\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
FST: 0.98836\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
213\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
10.350\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
0.04859\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
1.16\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
FCT: 0.76189\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\emph{N.} \emph{bulbuliferum}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Among populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
26\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
156.485\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
0.75095\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
93.73\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
FST: 0.93726\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Within populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
188\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
9.450\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
0.05027\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
6.27\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\emph{N. macrophyllum}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Among populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
2\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
57.921\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
3.10303\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
98.85\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
FST: 0.98853\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Within populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
25\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
0.900\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
0.03600\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
1.15\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\emph{N. thomsonianum}\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Among populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Within populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
-\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
All samples\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Among populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
30\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
569.516\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
2.32391\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
98.03\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
FST: 0.98033\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
Within populations\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
222\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
10.350\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
0.04662\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
1.97\strut
\end{minipage} & \begin{minipage}[t]{0.12\columnwidth}\raggedright\strut
\strut
\end{minipage}\tabularnewline
\bottomrule
\end{longtable}

Notes: FCT = differentiation among groups; FST = differentiation among
populations; FSC =

differentiation among populations within groups. ** P \textless{} 0.001.

\textbf{3.2 Phylogeny Reconstruction}

The ML and BI topologies of \emph{ITS} and cpDNA haplotype phylogenetic
trees were consistent on branches with strong support \textbf{(Figure
6)} . Partial inconsistencies between the \emph{ITS} and cpDNA trees
occurred regarding the phylogenetic relationships of \emph{N.
bulbuliferum} and \emph{N. macrophyllum} . In the \emph{ITS} haplotype
phylogenetic trees, \emph{N. bulbuliferum} and \emph{N.
macrophyllum}each had a haplotype (N7 and N26) embedded in the branch of
the other\textbf{(Figure 6a)} . In the cpDNA haplotype phylogenetic
trees, the haplotypes of all species are well clustered into a single
branch, forming a good monophyletic clade \textbf{(Figure 6b)} .
Overall, no haplotypes were shared among the three \emph{Notholirion}
species, and each species in general formed \emph{ITS} own individual
branch in the\emph{ITS} and cpDNA trees. The \emph{ITS} haplotype
network map and the\emph{ITS} haplotype phylogenetic tree constructed by
TCS network were consistent, with \emph{N. bulbuliferum} and \emph{N.
macrophyllum} each having a haplotype (N7 and N26) embedded in the
branch of the other\textbf{(Figure 3)} . while in cpDNA, the haplotype
network map and haplotype phylogenetic tree were consistent, with each
of the three species forming a separate monophyletic \textbf{(Figure 5)}
.

\textbf{3.3 Divergence Time Estimation}

As we all know, the molecular markers of chloroplast are closer to
neutral and have matrilineal genetic stability and widely used in the
process of divergence time estimation and biogeographic deduction. We
finally selected the results obtained based on common CDS for the
subsequent elaboration of the differentiation history
of\emph{Notholirion} . According to the beast-derived age estimate based
on the common CDS of 147 plastid genomes and 4 fossil calibration sites,
the crown age of \emph{Notholirion} was dated to be 25.05 Ma (95\%
HPD:35.56-16.62 Ma). The divergence times of \emph{N. thomsonianum} was
estimated at 7.43 Ma (95\% HPD:14.86-2.53 Ma). And The divergence
of\emph{N. macrophyllum} and \emph{N. bulbuliferum} occurred at 2.46 Ma
(95\% HPD:4.74-0.81 Ma). The results were shown in \textbf{Figures 10,
11} .

\textbf{3.4 Ancestral Area Reconstruction}

As the results of ancestral area reconstruction based on cpDNA phylogeny
showed, the most recent common ancestor of the \emph{Notholirion}species
probably originated in the Southern Himalayas\textbf{(region A, blue in
the Figure 9, Nodes I)} . Meanwhile, the most recent common ancestor
of\emph{N. macrophyllum} and \emph{N. bulbuliferum} most likely
originated in the Hengduan Mountains \textbf{(region C, brown in the
Figure 9, Nodes II)} . The origin of the genus \emph{Notholirion} is
hypothesized to have begun with initial divergence in the Southern
Himalayas, followed by migration of a portion of the ancestral taxa to
the Eastern Himalayas. From there, some populations continued to migrate
and spread eastward to the Hengduan Mountains, where they underwent
further divergence.

\textbf{3.5 Demographic History and Species Distribution Modeling}

Based on \emph{ITS} and cpDNA data, we performed mismatch distribution
analysis and neutrality test to determine the demographic history
of\emph{Notholirion} species. The mismatch distributions graphs for our
study species were multimodal and/or very ragged instead of a smooth
single-peaked curve, which indicates populations are stable and not
shrinking or expanding \textbf{(Figure 10)} . In addition, the values of
Tajima's D and Fu's Fs were non-significant for all
\emph{Notholirion}species in the neutrality test, which also means
populations have not experienced expansion or sustained growth
\textbf{(Table 3)} .

Considering \emph{Notholirion} and climate scenarios (the LIG, the LGM,
the present and the future), the area under the AUC value for the
potential climatically suitable areas of \emph{Notholirion} species was
high (\textgreater{}0.95), indicating the highest predictive capacity.
And the distribution ranges predicted for \emph{Notholirion} species
were consistent with the actual geographic distributions. Based on the
modeling results, there are slight variations in both the overall
distribution range and the optimal distribution range
of\emph{Notholirion} across four different periods. The overall
distribution range throughout these periods predominantly encompasses
the HHM region, while the optimal distribution range (highlighted in
red) is mainly concentrated in the Hengduan Mountains and along the
margin of the QTP. The suitable distribution area of
\emph{Notholirion}on the eastern part of the Tibetan Plateau contracted
significantly toward \emph{ITS} eastern edge during the LGM period, and
the main best distribution area shifted from the northern part of the
Hengduan Mountains to the southern region. And from the LGM period to
the present, the distribution area of \emph{Notholirion} re-expands to
the eastern part of the Tibetan Plateau. Predictions indicate that the
future overall distribution range and optimal distribution area will not
significantly differ from the present distribution \textbf{(Figure 11)}
.

\textbf{Table3 Results of the mismatch distribution analysis and
neutrality tests for\emph{Notholirion}}\selectlanguage{english}
\begin{longtable}[]{@{}lllll@{}}
\toprule
Floras & Tajima's \emph{D} (p) & Fu's \emph{Fs} (p) & SSD &
\emph{Hrag}\tabularnewline
\midrule
\endhead
\textbf{ITS} & \textbf{ITS} & \textbf{ITS} & \textbf{ITS} &
\textbf{ITS}\tabularnewline
\emph{N.} \emph{bulbuliferum} & 0.07080(0.52800) & -5.26370(0.08000) &
0.04602(0.26000) & 0.06807(0.41000)\tabularnewline
\emph{N. macrophyllum} & 2.86519(0.99500) & 12.88345(0.99999) &
0.20769(0.01000) & 0.34121(0.91000)\tabularnewline
\emph{N. thomsonianum} & - & - & - & -\tabularnewline
All & -0.49517(0.32500) & 0.55610(0.63000) & 0.03180(0.69000) &
0.04341(0.43000)\tabularnewline
\textbf{cpDNA} & \textbf{cpDNA} & \textbf{cpDNA} & \textbf{cpDNA} &
\textbf{cpDNA}\tabularnewline
\emph{N.} \emph{bulbuliferum} & 0.10895(0.54300) & -1.60875(0.33900) &
0.01034(0.16000) & 0.04323(0.25000)\tabularnewline
\emph{N. macrophyllum} & 2.23738(0.95100) & 8.15095(0.99999) &
0.27275(0.00001) & 0.55839(0.99999)\tabularnewline
\emph{N. thomsonianum} & - & - & - & -\tabularnewline
All & -0.74429(0.22900) & 2.67736(0.79700) & 0.03655 (0.08000) &
0.07168(0.06000)\tabularnewline
\bottomrule
\end{longtable}

\textbf{4 Discussion}

\textbf{4.1 Genetic Diversity and Structure of \emph{Notholirion}}

Biogenetic diversity is the result of the long-term evolution of species
and serves as an essential foundation for organisms to adapt to their
environment and evolve. It is also a key factor in maintaining the
stability and function of ecosystems. The higher the genetic diversity
of a species, the better it can adapt to environmental changes.
Conversely, when genetic diversity is low, a species is often vulnerable
to environmental changes and may lose \emph{ITS} distribution or even
become extinct(Soltis \& Soltis, 1991; Vranckx, Jacquemyn, Muys, \&
Honnay, 2012). Previous studies have reported that the average genetic
diversity of angiosperms, as inferred from plastid gene data, is
0.67(Petit et al., 2005). Our analysis of both plastid and
\emph{ITS}data sets revealed that the overall genetic diversity
(\emph{H\textsubscript{T}} value) of the \emph{Notholirion} was
significantly higher than this average (cpDNA: \emph{H\textsubscript{T}}
= 0.931; \emph{ITS} : \emph{H\textsubscript{T}} = 0.796), indicating a
remarkably high level of genetic diversity of this group (Table 1).
Through a series of analyses of the evolutionary aspects of the origin
of the \emph{Notholirion} , we suggested that two factors may have
contributed the most to the higher genetic diversity
of\emph{Notholirion} : (1) \emph{Notholirion} is the earliest divergent
clade of the tribe Lilieae, during the late Oligocene. During
\emph{ITS}long evolutionary history, it experienced numerous geological
and climatic events such as orogenic movements, monsoonal climatic
events, and drought events, which led to the accumulation of a large
amount of genetic variation to adapt to environmental changes;
(2)\emph{Notholirion} grows at high altitudes with an average altitude
of about 3000 meters, mainly in the Himalayan-Hengduan Mountains. The
area is influenced by frequent mountain-building movements in the
Miocene, resulting in an intricate and complex topography with
interlocking rows of high mountain valleys that form the ''sky
islands''. The unique ''sky island'' environment fragments habitat,
leading to alpine plant populations that form relatively isolated and
special ecosystems subject to different environmental selection and
ecological niche differentiation. This increases the opportunity for
genetic drift within populations, which results in a wealth of genetic
variation. Thus, the plant taxa distributed in this environment often
exhibit a high degree of genetic diversity and unique genotypic
assemblages, adapted to the challenges of these extreme environments.
Similar results have been found in other plants distributed in the HHM
region, such as\emph{Chamaesium} (cpDNA:\emph{H\textsubscript{T}} =
0.794; \emph{ITS} : \emph{H\textsubscript{T}} = 0.718),
\emph{Notopterygium incisum} (cpDNA:\emph{H\textsubscript{T}} = 0.939;
\emph{ITS} : \emph{H\textsubscript{T}}= 0.725) and \emph{Allium} section
\emph{Sikkimensia} (cpDNA:\emph{H\textsubscript{T}} = 0.974; \emph{ITS}
: \emph{H\textsubscript{T}} = 0.988)(Shahzad et al., 2017; C. Xie et
al., 2019; H. Y. Zheng et al., 2021).

In addition, AMOVA analyses were performed based on both \emph{ITS} and
cpDNA datasets. The results showed (Table2) that if \emph{Notholirion}is
considered as a whole, the genetic variation among the 31 populations
was large (\emph{ITS} : 91.45\%; cpDNA: 98.85\%) and within the
populations was relatively small (\emph{ITS} : 8.55\%; cpDNA: 1.15\%),
which may be related to the following reasons: (1) field surveys
revealed that although \emph{Notholirion} species are distributed in
four provinces of southwestern China, they have very specific habitat
requirements. Apart from altitude requirements, this genus also
necessitates specific heights of companion plants. For instance,\emph{N.
macrophyllum} and \emph{N. thomsonianum} prefer to grow in scrub or tall
grasses at around 3000m above sea level, where the plants are roughly
the same height as their own height. Hence, almost all populations
collected in the field have an extremely limited local distribution, and
communication between flora and populations is rare(J. Li et al., 2022);
(2) Through our field investigations, we also have observed that
\emph{Notholirion} species have very limited seed dispersal ability, and
their seed production cycle is quite lengthy. It takes approximately 5
years from seedling to flowering, and after flowering, the primary bulbs
of the roots wither, making it impossible for seed production to occur
in the following year. Mature seed production is only attainable through
newborn seeds and asexual reproduction of small bulbs, which take an
additional 5 years to reach the mature stage. This situation with seeds
also makes interflora and interpopulation communication very difficult;
(3) The genetic differentiation of many plant taxa in the HHM region has
been observed to be primarily influenced by the complex geological
activities and climate changes that occur in the area. For example,
orogenic movements can lead to fragmented habitats, which can reduce
gene flow among populations and consequently increase genetic variation
among them. Therefore, it is likely that the high genetic
differentiation observed among \emph{Notholirion} flora is due to the
significant fluctuations in the external environment in the HHM region
and \emph{ITS} own specificity. Similar results were found in
\emph{Chamaesium} ,\emph{Allium} section \emph{Sikkimensia} and
\emph{Polygonatum} (Xia et al., 2022; C. Xie et al., 2019; H. Y. Zheng
et al., 2021). However, AMOVA results of \emph{N. bulbuliferum} based on
\emph{ITS} data showed a lower percentage of genetic variation among
\emph{ITS} populations (74.86\%) compared to the results based on cpDNA
data (93.73). The reason for this inconsistent genetic structure is
hypothesized to be as follows: according to the \emph{ITS} data,
multiple haplotypes were detected in 27 \emph{N. bulbuliferum}
populations, of which N1 was a haplotype common to 15 populations, thus
reducing genetic variation among populations.

Both cpDNA and \emph{ITS} datasets showed high levels of haplotype
diversity in the haplotype analysis (cpDNA: Hd = 0.801; \emph{ITS} : Hd
= 0.788). There were no shared haplotypes among the
three\emph{Notholirion} species, and most haplotypes were restricted to
a single population or between several geographically adjacent
populations, which is similar to the results found in\emph{Chamaesium}
(H. Y. Zheng et al., 2021). We speculate that\emph{Notholirion} , as a
high-altitude distributed taxon, has different populations that remain
isolated and lack communication with each other during ice ages and
interglacial periods. This situation may be due to several environmental
factors and the species \emph{ITS} elf, including the inability of
high-altitude plants to adapt to low-altitude environments, the
fragmentation of habitats in the HHM region, and the formation of unique
habitats such as sky islands that make communication between populations
difficult, and the long seed production cycle and difficulty of
dissemination of \emph{Notholirion} species themselves that limit
communication between populations.

\textbf{4.2 Phylogenetic relationships of \emph{Notholirion}}

The phylogenetic analysis using both nuclear gene fragments and plastid
fragments confirmed the monophyletic nature of \emph{Notholirion} ,
which is consistent with previous studies(Huang et al., 2018; Kim \&
Kim, 2018; J. Li et al., 2022). Within the genus \emph{Notholirion} ,
the species \emph{N. thomsonianum} was found to be a sister taxon to the
branches of \emph{N. macrophyllum} and \emph{N. bulbuliferum} ,
representing the earliest divergence event as supported by the BEAST
divergence time estimates. Notably, some populations showed
inconsistencies in their nuclear genetic makeup. The cpDNA tree revealed
that all populations of the three \emph{Notholirion} species formed
distinct monophyletic branches. However, in the \emph{ITS} tree, two
populations of \emph{N. macrophyllum} (MNC and MZT) clustered together
with most populations of \emph{N. bulbuliferum} , while \emph{N.
bulbuliferum} exhibited a population (BDY) clustered into a single
branch with the remaining two populations of \emph{N. macrophyllum} .

The phenomenon of inconsistency between nuclear gene trees and
chloroplast gene trees is very common in plant phylogenetic analyses and
is typically attributed to convergent evolution, incomplete lineage
sorting, hybridization, or gene introgression(Y. Gao, Harris, Li, \&
Gao, 2020; Y. D. Gao, Harris, \& He, 2015; Y. D. Gao, Harris, Zhou, \&
He, 2013; Huang et al., 2018; X. Liu, Wang, Shao, Ye, \& Zhang, 2016).
The selected gene fragments in this study are mostly evolutionarily
conserved and less influenced by the environment, reducing the
likelihood of convergent evolution. Therefore, the main factors
considered are incomplete lineage sorting and hybridization. And in many
studies on the phylogenetic relationships of plants in the family
Liliaceae, both incomplete lineage sorting and hybridization have been
mentioned as factors contributing to conflicts in the phylogenetic
relationships(Y. Gao et al., 2020; Y. D. Gao et al., 2015; Y. D. Gao et
al., 2013; Huang et al., 2018). Studies have shown that incomplete
lineage sorting is more likely to occur when species rapidly expand or
have large population sizes(Maddison, 1997). Due to the stochastic
nature of the coalescent process, incomplete lineage sorting can lead to
random patterns of gene trees among taxa, which may result in
inconsistencies between different gene trees(Buckley, Cordeiro,
Marshall, \& Simon, 2006). Furthermore, lineage sorting in plants
generally takes hundreds of millions of years(X. Liu et al., 2016;
Pelser et al., 2010). Both \emph{N. macrophyllum} and \emph{N.
bulbuliferum} are relatively young species that diverged in the last few
million years, which is far from reaching the timescale of hundreds of
millions of years. Therefore, we speculate that incomplete lineage
sorting may be one of the reasons for the observed incongruence between
nuclear and chloroplast genomes in \emph{Notholirion} . We know that
hybridization plays a significant role in biodiversity and speciation by
either reinforcing or disrupting reproductive barriers(Les et al.,
2015). Due to the maternal inheritance of organelle genes, gene flow
resulting from biparental hybridization is usually manifested only in
nuclear genes, leading to inconsistencies between the phylogenetic trees
constructed based on nuclear genes and organellar genes(X. Liu et al.,
2016). Hybridization is widely observed in the family Liliaceae, as
evidenced by numerous reports in wild \emph{Lilium} species
and\emph{ITS} frequent use in horticulture for plant breeding
purposes(Y. Gao et al., 2020; Y. D. Gao et al., 2015; Y. D. Gao et al.,
2013; Pelkonen, Niittyvuopio, Pirttil\selectlanguage{ngerman}ä, Laine, \& Hohtola, 2007).
Furthermore, some researchs have suggested that hybridization may be a
common phenomenon in Liliaceae plants(Douglas et al., 2011; Muratovi\selectlanguage{polish}ć,
\selectlanguage{english}Robin, Boguni\selectlanguage{polish}ć, Š\selectlanguage{english}oljan, \& Siljak-Yakovlev, 2010). In conclusion, we
speculate that incomplete lineage sorting and hybridization may be the
primary factors leading to the conflicting phylogenetic trees observed
between the nuclear genes and chloroplast genes in \emph{Notholirion} .

\textbf{4.3 Origin and Diversification of \emph{Notholirion}}

In the vast timeline of Earth's evolution, every geological event and
climate change has had unprecedented effects on organisms, making it a
highly researched topic in the present day. Based on the complete
plastid genomic dataset, we estimated the origin of the
genus\emph{Notholirion} and the divergence times of the three species
within the genus \textbf{(Figure 7, 8)} . Additionally, we reconstructed
their ancestral distribution using the chloroplast fragments
\textbf{(Figure 9)} . The results suggest that the ancestral lineage of
the\emph{Notholirion} genus originated in the southern Himalayas during
the end of the Tertiary period in the Cenozoic era (25.05 Ma, 95\% HPD:
35.56-16.62 Ma). The prevailing theory regarding this period suggests
that from 38.0-23.03 Ma, the Earth transitioned from the end of the
Eocene epoch through the Oligocene epoch and into the early Miocene
epoch. During this time, the climate gradually cooled, resulting in a
moderate icehouse climate, and plants and animals evolved and
diversified rapidly(O'Brien et al., 2020; Straume, Nummelin, Gaina, \&
Nisancioglu, 2022). It is during this period that the evolution and
spread of modern types of flowering plants primarily took place(Nge,
Biffin, Thiele, \& Waycott, 2020). Furthermore, during this period, as a
result of the gradual collision between the Indian subcontinent and the
Eurasian Plate, the Qinghai-Tibetan Plateau continued to uplift both to
the south and north. The ongoing uplift of the Qinghai-Tibet Plateau has
caused dramatic changes in the topography of the region and has directly
contributed to the formation of the monsoon climate(Favre et al., 2015).
The significant changes in habitat have also facilitated the formation
and differentiation of populations, as seen in \emph{Chamaesium}
,\emph{Cardiocrinum} , and \emph{Sinopteris} (J. Li et al., 2022; L.
Wang, Yang, Zhang, Zhang, \& Zhang, 2023; H. Y. Zheng et al., 2021).
Therefore, we speculate that the origin of the genus
\emph{Notholirion}is closely related to the uplift events that occurred
on Tibetan Plateau.

During the early to middle Miocene (23.03-10 Ma), the Earth's climate
experienced significant cooling and entered a period of massive
glaciation. The uplift of the Himalayas and surrounding mountains caused
aridification of the local climate and a decrease in temperature,
limiting plant growth and reproduction and having a profound impact on
the evolution of local species, such as \emph{Gentiana crassicaulis}
,\emph{Polygonatum} and \emph{Allium} section \emph{Daghestanica} (M. J.
Li, Yu, Guo, \& He, 2021; Ni, Li, Zhao, Gaawe, \& Liu, 2022; Xia et al.,
2022). By the late Miocene (10-5.33 Ma), the Qinghai-Tibetan Plateau
experienced a more significant uplift and expansion, which affected
numerous mountain ranges in the plateau's eastern margin region. Among
them, the Himalayas and the Gangdise Mountains run parallel to each
other in a northwest-southeast direction, forming a natural barrier.
This allows the warm and humid airflow from the Indian Ocean to enter
China exclusively through the Hengduan Mountains. This has brought
abundant rainfall to the southeastern Tibetan Plateau, significantly
impacting species in that region, such as \emph{Cardiocrinum}
,\emph{Chamaesium} \emph{thalictrifolium} , \emph{Allium
sectionSikkimensia} , and\emph{Rabdosia} (J. Li et al., 2022; C. Xie et
al., 2019; X. Q. Yu et al., 2014; H. Y. Zheng et al., 2021)\emph{.} Our
molecular dating analysis suggested that \emph{N. thomsonianum} diverged
during the late Miocene (7.43 Ma, 95\% HPD: 14.86-2.53 Ma). Therefore,
we infer that climate change and complex geological activity play
important roles in the differentiation and dispersal
of\emph{Notholirion} species during this period.

During the Pleistocene period (2.588-0.126 Ma), the most far-reaching
geological event was the Last Ice Age, a repeated cycle of alternating
ice ages and interglacials that caused widespread climatic upheaval with
unprecedented effects on the biology of the period(G. M. Hewitt, 2004).
During the Ice Age, a continuous ice cap did not form on the Tibetan
Plateau and surrounding areas, especially in the Hengduan Mountains
located on \emph{ITS} eastern edge. Intense orogenic movements in the
Hengduan Mountains during and after the late Miocene caused the
mountains to stretch north and south, resulting in high mountain valleys
and diverse topography characterized by folded mountains and faulted
basins. While the ice age caused most of the region to become cold and
dry, there were still relatively warm and moist areas in the Hengduan
Mountains region where organisms could survive. The harsh conditions of
the ice age forced organisms to seek refuge in these areas, resulting in
natural isolation that may have led to the divergence and formation of
new species(Q. Meng et al., 2022). Through BEAST analysis, we discovered
that the divergence of\emph{N. macrophyllum} and \emph{N. bulbuliferum}
occurred during the Last Glacial Maximum (2.46 Ma, 95\% HPD: 4.74-0.81
Ma), indicating that the significant climate fluctuations and natural
isolation between different refuges may play a crucial role in the
formation and distribution of these species.

\textbf{4.4 Population dynamics history within \emph{Notholirion}}

The Quaternary Ice Age has been a highly significant and widely studied
event in the field of biogeography. The climatic upheavals associated
with this ice age, along with interglacial cycles, have had a profound
impact on species ranges, often leading to population migration and even
extinction(G. M. Hewitt, 2004). Based on our previous research on
divergence timing, we have found that among the three species in
the\emph{Notholirion} , \emph{N. thomsonianum} diverged as the early
branch during the late Miocene, whereas the other two species did not
fully diverge until the LGM period. Therefore, we speculate that the
climatic fluctuations during the Quaternary Ice Age, particularly the
significant temperature drop during the LGM, played a crucial role in
shaping the range of the ancestral taxa of the \emph{Notholirion} genus.
These changes likely in turn contributed to the divergence of \emph{N.
bulbuliferum} and \emph{N. macrophyllum} , leading to the formation of
distinct species. According to the results predicted by our distribution
models for the three periods (LIG, LGM and present), the suitable
distribution area of \emph{Notholirion} in the eastern part of the
Tibetan Plateau contracted significantly toward \emph{ITS} eastern edge
during the Last Glacial Maximum (LGM), and then re-expanded toward the
eastern part of the Plateau to \emph{ITS} present distribution pattern
after the LGM. Additionally, our mismatch distribution analysis based on
both \emph{ITS} and cpDNA data exhibited non-smooth bimodal or
multi-peak curves (Figure). This, coupled with the non-star-like
phylogeographical structure of haplotypes and non-significant results of
the neutrality test, suggests that \emph{Notholirion} had not recently
undergone rapid expansion events (Figure, Table). We propose that the
stability of \emph{Notholirion} in \emph{ITS} overall distribution
range, as well as the observed contraction and re-expansion in the
eastern Tibetan Plateau, can be attributed to \emph{ITS} distinctive
characteristics as an autochthonous high-altitude taxon, along with the
combined influence of external environmental factors. These findings
also indicate that the response of plant taxa across the eastern Tibetan
Plateau and the Hengduan Mountains to Quaternary ice age climate change
is more intricate and nuanced than previously assumed.

Numerous studies have provided evidence that during the Quaternary
period, there was no uniform ice cap covering both the Tibetan Plateau
and the Hengduan Mountains(G. M. Hewitt, 2004; Taberlet, Fumagalli,
Wust-Saucy, \& Cosson, 1998). Additionally, \emph{Notholirion} , as a
high-altitude plant taxon primarily found at elevations around 3000 m,
exhib\emph{ITS} enhanced cold and frost tolerance and a greater capacity
to adapt to the environmental fluctuations of the ice age. Consequently,
it is plausible that certain widely distributed and cold-tolerant plant
taxa, such as \emph{Chamaesium} , \emph{Allium}
section\emph{Sikkimensia} , and \emph{Cardiocrinum} ,
including\emph{Notholirion} , survived the ice age by persisting in
multiple refugia and subsequently expanded their distribution following
the ice age(J. Li et al., 2022; C. Xie et al., 2019; H. Y. Zheng et al.,
2021). Furthermore, we observed notable levels of genetic diversity
(\emph{ITS} : HD\textgreater{}0.7) and unique haplotypes specific in
populations from BDY, BLQ, and BZD. In the context of genealogical
geographic studies, regions exhibiting elevated genetic diversity or the
presence of ancient and endemic haplotypes are often regarded as
potential ice age refuges for plants(Tzedakis, Lawson, Frogley, Hewitt,
\& Preece, 2002). Based on these findings, we propose that the Hengduan
Mountains and their neighboring areas could have served as refuges
during the ice age for \emph{Notholirion} species inhabiting both the
eastern Tibetan Plateau and the Hengduan Mountains. Notably, this region
has been recognized as a significant refuge for numerous plant taxa on
the Tibetan Plateau, including \emph{Potentilla glabra} , \emph{Rhodiola
alsia} , \emph{Eriophyton wallichii} , and \emph{Chamaesium} (Q. B. Gao
et al., 2012; L. Y. Wang, Ikeda, Liu, Wang, \& Liu, 2009; X. X. Wang,
Yue, Sun, \& Li, 2011; H. Y. Zheng et al., 2021). These findings suggest
that the Quaternary ice age may have had less impact on the survival and
widespread distribution of these high-elevation cold-tolerant plants in
the subtropical region of the Tibetan Plateau than we thought. They
survived on multiple ice-age refuges and rapidly regained their
pre-glacial ranges after the ice age.

\textbf{5 Conclusion}

In conclusion, our study on the genetic diversity and structure
of\emph{Notholirion} sheds light on the intricate evolutionary history
and adaptive mechanisms of this genus. The observed high genetic
diversity of \emph{Notholirion} , surpassing the average of angiosperms,
is attributed to \emph{ITS} early divergence in the late Oligocene
and\emph{ITS} adaptation to diverse climatic and geological events
during\emph{ITS} evolution. The unique ''sky island'' environments,
resulting from complex topography in the Himalayan-Hengduan Mountains,
contribute to the isolation of populations and, consequently, high
genetic differentiation. The AMOVA analyses revealed substantial genetic
variation among \emph{Notholirion} populations, indicating limited
communication between populations due to specific habitat requirements,
restricted seed dispersal, and environmental factors influencing gene
flow. This genetic differentiation is consistent with the impact of
complex geological activities and climate changes in the
Himalayan-Hengduan Mountains region. Phylogenetic analysis confirmed the
monophyletic nature of \emph{Notholirion} , with \emph{N.
thomsonianum}identified as the earliest divergent taxon. However,
incongruences between nuclear and chloroplast gene trees suggest the
influence of incomplete lineage sorting and hybridization, common
phenomena in plant phylogenetics. The origin and diversification of
\emph{Notholirion} are closely tied to the uplift events on the Tibetan
Plateau. The late Miocene and Pleistocene periods played pivotal roles
in the evolution and distribution of \emph{Notholirion} species, with
significant climate fluctuations contributing to species divergence. The
Last Glacial Maximum, in particular, led to the formation of distinct
species, emphasizing the impact of Quaternary Ice Age climate changes on
plant ranges and population dynamics. Our findings propose
that\emph{Notholirion} , as a high-altitude taxon, survived the ice age
through persistence in multiple refugia, particularly in the Hengduan
Mountains region. The study challenges previous assumptions about the
uniformity of ice cap coverage in the Tibetan Plateau and Hengduan
Mountains during the Quaternary period, suggesting a more nuanced
response of cold-tolerant plant taxa to climate changes. Overall, this
research contributes valuable insights into the evolutionary and
ecological dynamics of \emph{Notholirion} , highlighting the importance
of considering multiple factors, including geological events, climatic
changes, and specific ecological, in understanding the genetic diversity
and adaptive strategies of high-altitude plant taxa.

\textbf{Acknowledgments}

This work was supported by the National Natural Science Foundation of
China (Grant Nos. 32100180, 32070221, 32170209), and Survey on the
Background Resources of Chengdu Area of Giant Panda National Park
(Project No. 510101202200376).

\textbf{Author contribution statemen}

Rui-Yu Cheng: Conceptualization, Formal analysis, Investigation,
Methodology, Writing--original draft. Juan Li: Conceptualization,
Investigation, Data curation, Formal analysis. Deng-Feng Xie: Formal
analysis, Investigation, Methodology. Ren-Xiu Zhou: Formal analysis,
Investigation. Qing Li: Funding acquisition. Yanglina Yu: Funding
acquisition. Xing-Jin He: Investigation, Validation, Writing--review \&
editing. Song-Dong Zhou: Investigation, Validation, Writing--review \&
editing.

\textbf{Declaration of Interest Statement}

The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to
influence the work reported in this paper.

\textbf{Data availability}

The datasets presented in this study can be found in online
repositories. The names of the repositories and accession number(s) can
be found below: https://www.ncbi.nlm.nih.gov/genbank/, (chloroplast
genomes) OR899315 - OR899324; (\emph{ITS} ) OR898935 - OR899150 and
OR899157 - OR899194, (\emph{matK} ) OR914272 - OR914525; (\emph{ndhA} )
OR897910 - OR898163; (\emph{ndhG-I} ) OR914680 - OR914933;
(\emph{petB-D} ) OR917093 - OR917346; (\emph{petL-G} ) OR917454 -
OR917707.

\textbf{Reference} :

Aitchison, J. C., \& Ali, J. R. (2008). When and where did India and
Asia collide? \emph{Geol Bull China} (9).

Alemu, W. G., \& Wimberly, M. C. (2020). Evaluation of Remotely Sensed
and Interpolated Environmental Datasets for Vector-Borne Disease
Monitoring Using In Situ Observations Over the Amhara Region,
Ethiopia.\emph{Sensors (Basel), 20} (5). doi:10.3390/s20051316

Bandelt, H. J., Forster, P., \& R\selectlanguage{ngerman}öhl, A. (1999). Median-joining networks
for inferring intraspecific phylogenies. \emph{Mol Biol Evol, 16} (1),
37-48. doi:10.1093/oxfordjournals.molbev.a026036

Borns, H. W. (1994). \emph{The ice age world : an introduction to
quaternary history and research with emphasis on North America and
Northern Europe during the last 2.5 million years} : Scandinavian
University Press.

Buckley, T. R., Cordeiro, M., Marshall, D. C., \& Simon, C. (2006).
Differentiating between hypotheses of lineage sorting and introgression
in New Zealand alpine cicadas (Maoricicada Dugdale). \emph{Syst Biol,
55} (3), 411-425. doi:10.1080/10635150600697283

Buerki, S., Jose, S., Yadav, S. R., Goldblatt, P., Manning, J. C., \&
Forest, F. (2012). Contrasting biogeographic and diversification
patterns in two Mediterranean-type ecosystems. \emph{PLoS One, 7} (6),
e39377. doi:10.1371/journal.pone.0039377

Cao, J. L., Shi, Y. L., Zhang, H., \& Wang, H. (2009). Numerical
simulation of GPS observed clockwise rotation around the eastern
Himalayan syntax in the Tibetan Plateau. \emph{Chin Sci Bull} .

Chatterjee, S., Goswami, A., \& Scotese, C. R. (2013). The longest
voyage: Tectonic, magmatic, and paleoclimatic evolution of the Indian
plate during its northward flight from Gondwana to Asia. \emph{Gondwana
Res, 23} (1), 238-267.

Chen, J. M., Zhi-Yuan, D. U., Yuan, Y. Y., \& Wang, Q. F. (2014).
Phylogeography of an alpine aquatic herb Ranunculus bungei
(Ranunculaceae) on the Qinghai--Tibet Plateau. \emph{J Syst Evol, 52}
(003), 313-325.

Chen, S. C., Liang, S. J., Xu, J. M., \& Minoru, N. T. (2000). Liliaceae
A. L. Jussieu. \emph{Flora of China, 14} , 73-263.

Deng, T. D., Tao., Ding, L. D., Lin., \& Deng, T. (2015). Paleoaltimetry
reconstructions of the Tibetan Plateau: progress and
contradictions.\emph{Natl Sci Rev, 2} (4), 417-437.

Dierckxsens, N., Mardulyn, P., \& Smits, G. (2017). NOVOPlasty: de novo
assembly of organelle genomes from whole genome data. \emph{Nucleic
Acids Res, 45} (4), e18. doi:10.1093/nar/gkw955

Douglas, N. A., Wall, W. A., Xiang, Q. Y., Hoffmann, W. A., Wentworth,
T. R., Gray, J. B., \& Hohmann, M. G. (2011). Recent vicariance and the
origin of the rare, edaphically specialized Sandhills lily, Lilium
pyrophilum (Liliaceae): evidence from phylogenetic and coalescent
analyses. \emph{Mol Ecol, 20} (14), 2901-2915.
doi:10.1111/j.1365-294X.2011.05151.x

Dupont-Nivet, G., Lippert, P. C., Hinsbergen, D., Meijers, M., \& Kapp,
P. (2010). Paleolatitude and age of the Indo-Asia collision:
Paleomagnetic constraints. \emph{Geophys J Int, 182} (3), 1189-1198.

E., P., \& C., G. (1980). A rapid DNA isolation procedure for small
quantities of fresh leaf tissue. \emph{Phytochemistry} (19), 11-13.

Elith, J., \& Leathwick, J. R. (2009). Species Distribution Models:
Ecological Explanation and Prediction Across Space and Time. \emph{Annu
Rev Ecol Evol Syst, 40} (1), 677-697.

Ellis, C. J., Yahr, R., Belinchón, R., \& Coppins, B. J. (2014).
Archaeobotanical evidence for climate as a driver of ecological
community change across the anthropocene boundary. \emph{Glob Chang
Biol, 20} (7), 2211-2220. doi:10.1111/gcb.12548

Excoffier, L., \& Lischer, H. E. (2010). Arlequin suite ver 3.5: a new
series of programs to perform population genetics analyses under Linux
and Windows. \emph{Mol Ecol Resour, 10} (3), 564-567.
doi:10.1111/j.1755-0998.2010.02847.x

Favre, A., Päckert, M., Pauls, S. U., Jähnig, S. C., Uhl, D., Michalak,
I., \& Muellner-Riehl, A. N. (2015). The role of the uplift of the
Qinghai-Tibetan Plateau for the evolution of Tibetan biotas. \emph{Biol
Rev Camb Philos Soc, 90} (1), 236-253. doi:10.1111/brv.12107

Fawcett, T. (2005). An introduction to ROC analysis. \emph{Pattern
Recognit Lett, 27} (8), 861-874.

Fick, S. E., \& Hijmans, R. J. (2017). WorldClim 2: New 1-km spatial
resolution climate surfaces for global land areas. \emph{Int J Climatol,
37} (12).

Fielding, A. (1997). A review of methods for the assessment of
prediction errors in conservation presence/absence models. \emph{Environ
Conser, 24} (1), 38-49.

Gao, Q. B., Zhang, D. J., Duan, Y. Z., Zhang, F. Q., Yinhu, L. I.,
Pengcheng, F. U., \& Chen, S. L. (2012). Intraspecific divergences of
Rhodiola alsia (Crassulaceae) based on plastid DNA and internal
transcribed spacer fragments. \emph{BIOL J LINN SOC} .

Gao, Y., Harris, A. J., Li, H., \& Gao, X. (2020). Hybrid Speciation and
Introgression Both Underlie the Genetic Structures and Evolutionary
Relationships of Three Morphologically Distinct Species of Lilium
(Liliaceae) Forming a Hybrid Zone Along an Elevational
Gradient.\emph{Front Plant Sci, 11} , 576407.
doi:10.3389/fpls.2020.576407

Gao, Y. D., Harris, A. J., \& He, X. J. (2015). Morphological and
ecological divergence of Lilium and Nomocharis within the Hengduan
Mountains and Qinghai-Tibetan Plateau may result from habitat
specialization and hybridization. \emph{BMC Evol Biol, 15} , 147.
doi:10.1186/s12862-015-0405-2

Gao, Y. D., Harris, A. J., Zhou, S. D., \& He, X. J. (2013).
Evolutionary events in Lilium (including Nomocharis, Liliaceae) are
temporally correlated with orogenies of the Q-T plateau and the Hengduan
Mountains. \emph{Mol Phylogenet Evol, 68} (3), 443-460.
doi:10.1016/j.ympev.2013.04.026

Gernhard, T. (2008). The conditioned reconstructed process. \emph{J
Theor Biol, 253} (4), 769-778. doi:10.1016/j.jtbi.2008.04.005

Helfrich, P., Rieb, E., Abrami, G., Lücking, A., \& Mehler, A.
(2018).\emph{TreeAnnotator: versatile visual annotation of hierarchical
text relations.} Paper presented at the Proceedings of the eleventh
international conference on language resources and evaluation (LREC
2018).

Hewitt, G. (2000). The genetic legacy of the Quaternary ice
ages.\emph{Nature, 405} (6789), 907-913. doi:10.1038/35016000

Hewitt, G. M. (1996). Some genetic consequences of ice ages, and their
role in divergence and speciation. \emph{BIOL J LINN SOC} (3), 247-276.

Hewitt, G. M. (2004). Genetic consequences of climatic oscillations in
the Quaternary. \emph{Philos Trans R Soc Lond B Biol Sci, 359} (1442),
183-195; discussion 195. doi:10.1098/rstb.2003.1388

Hickerson, M. J., Carstens, B. C., Cavender-Bares, J., Crandall, K. A.,
Graham, C. H., Johnson, J. B., . . . Yoder, A. D. (2010).
Phylogeography's past, present, and future: 10 years after Avise,
2000.\emph{Mol Phylogenet Evol, 54} (1), 291-301.
doi:10.1016/j.ympev.2009.09.016

Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., \& Jarvis,
A. (2010). Very high resolution interpolated climate surfaces for global
land areas. \emph{Int J Climatol, 25} (15), 1965-1978.

Huang, J., Yang, L. Q., Yu, Y., Liu, Y. M., Xie, D. F., Li, J., . . .
Zhou, S. D. (2018). Molecular phylogenetics and historical biogeography
of the tribe Lilieae (Liliaceae): bi-directional dispersal between
biodiversity hotspots in Eurasia. \emph{Ann Bot, 122} (7), 1245-1262.
doi:10.1093/aob/mcy138

Huntington, K. W., Blythe, A. E., \& Hodges, K. V. (2006). Climate
change and Late Pliocene acceleration of erosion in the
Himalaya.\emph{Earth Planet Sci Lett, 252} (1-2), 107-118.

Kai, H., \& Jiang, X. (2014). Sky islands of southwest China. I: An
overview of phylogeographic patterns. \emph{Chin Sci Bull, 59} (7),
1-13.

Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., von Haeseler, A., \&
Jermiin, L. S. (2017). ModelFinder: fast model selection for accurate
phylogenetic estimates. \emph{Nat Methods, 14} (6), 587-589.
doi:10.1038/nmeth.4285

Katoh, K., \& Standley, D. M. (2013). MAFFT multiple sequence alignment
software version 7: improvements in performance and usability. \emph{Mol
Biol Evol, 30} (4), 772-780. doi:10.1093/molbev/mst010

Kearse, M., Moir, R., Wilson, A., Stones-Havas, S., Cheung, M.,
Sturrock, S., . . . Drummond, A. (2012). Geneious Basic: an integrated
and extendable desktop software platform for the organization and
analysis of sequence data. \emph{Bioinformatics, 28} (12), 1647-1649.
doi:10.1093/bioinformatics/bts199

Keller, I., \& Seehausen, O. (2012). Thermal adaptation and ecological
speciation. \emph{Mol Ecol, 21} (4), 782-799.
doi:10.1111/j.1365-294X.2011.05397.x

Kim, J. S., \& Kim, J. H. (2018). Updated molecular phylogenetic
analysis, dating and biogeographical history of the lily family
(Liliaceae: Liliales). \emph{Bot. J. Linn. Soc., 187} (4), 579-593.

Kumar, S., Stecher, G., Li, M., Knyaz, C., \& Tamura, K. (2018). MEGA X:
Molecular Evolutionary Genetics Analysis across Computing
Platforms.\emph{Mol Biol Evol, 35} (6), 1547-1549.
doi:10.1093/molbev/msy096

Lanfear, R., Frandsen, P. B., Wright, A. M., Senfeld, T., \& Calcott, B.
(2017). PartitionFinder 2: New Methods for Selecting Partitioned Models
of Evolution for Molecular and Morphological Phylogenetic
Analyses.\emph{Mol Biol Evol, 34} (3), 772-773.
doi:10.1093/molbev/msw260

Les, D. H., Peredo, E. L., King, U. M., Benoit, L. K., Tippery, N. P.,
Ball, C. J., \& Shannon, R. K. (2015). Through thick and thin: cryptic
sympatric speciation in the submersed genus Najas
(Hydrocharitaceae).\emph{Mol Phylogenet Evol, 82 Pt A} , 15-30.
doi:10.1016/j.ympev.2014.09.022

Li, D. Z. (2008). Floristics and plant biogeography in China. \emph{J
Integr Plant Biol, 50} (7), 771-777.
doi:10.1111/j.1744-7909.2008.00711.x

Li, H. T., Yi, T. S., Gao, L. M., Ma, P. F., Zhang, T., Yang, J. B., . .
. Li, D. Z. (2019). Origin of angiosperms and the puzzle of the Jurassic
gap. \emph{Nat Plants, 5} (5), 461-470. doi:10.1038/s41477-019-0421-0

Li, J., Cai, J., Qin, H. H., Price, M., Zhang, Z., Yu, Y., . . . Gao, X.
F. (2022). Phylogeny, Age, and Evolution of Tribe Lilieae (Liliaceae)
Based on Whole Plastid Genomes. \emph{Front Plant Sci, 12} , 699226.
doi:10.3389/fpls.2021.699226

Li, J., Milne, R. I., Ru, D., Miao, J., Tao, W., Zhang, L., . . . Mao,
K. (2020). Allopatric divergence and hybridization within Cupressus
chengiana (Cupressaceae), a threatened conifer in the northern Hengduan
Mountains of western China. \emph{Mol Ecol, 29} (7), 1250-1266.
doi:10.1111/mec.15407

Li, M. J., Yu, H. X., Guo, X. L., \& He, X. J. (2021). Out of the
Qinghai-Tibetan Plateau and rapid radiation across Eurasia for Allium
section Daghestanica (Amaryllidaceae). \emph{AoB Plants, 13} (3),
plab017. doi:10.1093/aobpla/plab017

Li, Q., Guo, X., Niu, J., Duojie, D., Li, X., Opgenoorth, L., \& Zou, J.
(2020). Molecular Phylogeography and Evolutionary History of the Endemic
Species Corydalis hendersonii (Papaveraceae) on the Tibetan Plateau
Inferred From Chloroplast DNA and ITS Sequence Variation. \emph{Front
Plant Sci, 11} , 436. doi:10.3389/fpls.2020.00436

Liu, H. R., Khan, G., Gao, Q., Zhang, F., Liu, W., Wang, Y., . . .
Afridi, S. G. (2022). Dispersal into the Qinghai-Tibet plateau: evidence
from the genetic structure and demography of the alpine plant Triosteum
pinnatifidum. \emph{PeerJ, 10} , e12754. doi:10.7717/peerj.12754

Liu, J. Q., Duan, Y., Hao, G., Xue\selectlanguage{english}mgUn, G., \& Sun, H. (2014).
Evolutionary history and underlying adaptation of alpine plants on the
Qinghai-Tibet Plateau. \emph{J Syst Evol, 52} (3), 241-249.

Liu, J. Q., Sun, Y. S., Ge, X. J., Gao, L. M., \& Qiu, Y. X. (2012).
Phylogeographic studies of plants in China: Advances in the past and
directions in the future. \emph{J Biogeogr, 50} (3), 267-275.

Liu, M. L., He, Y. L., Lopez-Pujol, J., Jia, Y., \& Li, Z. H. (2019).
Complex population evolutionary history of four cold-tolerant
Notopterygium herb species in the Qinghai-Tibetan Plateau and adjacent
areas. \emph{Heredity (Edinb), 123} (2), 242-263.
doi:10.1038/s41437-019-0186-2

Liu, X., Wang, Z., Shao, W., Ye, Z., \& Zhang, J. (2016). Phylogenetic
and Taxonomic Status Analyses of the Abaso Section from Multiple Nuclear
Genes and Plastid Fragments Reveal New Insights into the North America
Origin of Populus (Salicaceae). \emph{Front Plant Sci, 7} , 2022.
doi:10.3389/fpls.2016.02022

Luo, D., Yue, J. P., \& Sun, W. G. (2016). Evolutionary history of the
subnival flora of the Himalaya-Hengduan Mountains: first insights from
comparative phylogeography of four perennial herbs. \emph{J Biogeogr} .

Maddison, W. P. (1997). Gene trees in species trees. \emph{Syst. Biol.,
46} (3), 523-536.

Marchese, C. (2015). Biodiversity hotspots: A shortcut for a more
complicated concept. \emph{Glob Ecol Conserv, 3} (74), 297-309.

Meng, H. H., Su, T., Gao, X. Y., Li, J., Jiang, X. L., Sun, H., \& Zhou,
Z. K. (2017). Warm-cold colonization: response of oaks to uplift of the
Himalaya-Hengduan Mountains. \emph{Mol Ecol, 26} (12), 3276-3294.
doi:10.1111/mec.14092

Meng, Q., Xie, Z., Xu, H., Guo, J., Tang, Y., Ma, T., . . . Gao, T.
(2022). Out of the Qinghai-Tibetan plateau: Origin, evolution and
historical biogeography of Morchella (both Elata and Esculenta
clades).\emph{Front Microbiol, 13} , 1078663.
doi:10.3389/fmicb.2022.1078663

Mulch, A., \& Chamberlain, C. P. (2006). Earth science: the rise and
growth of Tibet. \emph{Nature, 439} (7077), 670-671. doi:10.1038/439670a

Muratovi\selectlanguage{polish}ć, \selectlanguage{english}E., Robin, O., Boguni\selectlanguage{polish}ć, \selectlanguage{english}F., \selectlanguage{polish}Š\selectlanguage{english}oljan, D., \& Siljak-Yakovlev,
S. (2010). Karyotype evolution and speciation of European lilies from
Lilium sect. Liriotypus. \emph{Taxon, 59} (1), 165-175.

Nei, M., \& Li, W. H. (1979). Mathematical model for studying genetic
variation in terms of restriction endonucleases. \emph{Proc Natl Acad
Sci U S A, 76} (10), 5269-5273. doi:10.1073/pnas.76.10.5269

Nei, M., \& Tajima, F. (1981). DNA polymorphism detectable by
restriction endonucleases. \emph{Genetics, 97} (1), 145-163.
doi:10.1093/genetics/97.1.145

Nge, F. J., Biffin, E., Thiele, K. R., \& Waycott, M. (2020). Extinction
pulse at Eocene-Oligocene boundary drives diversification dynamics of
two Australian temperate floras. \emph{Proc Biol Sci, 287} (1919),
20192546. doi:10.1098/rspb.2019.2546

Ni, L., Li, W., Zhao, Z., Gaawe, D., \& Liu, T. (2022). Migration
patterns of Gentiana crassicaulis, an alpine gentian endemic to the
Himalaya-Hengduan Mountains. \emph{Ecol Evol, 12} (3), e8703.
doi:10.1002/ece3.8703

O'Brien, C. L., Huber, M., Thomas, E., Pagani, M., Super, J. R., Elder,
L. E., \& Hull, P. M. (2020). The enigma of Oligocene climate and global
surface temperature evolution. \emph{Proc Natl Acad Sci U S A, 117}
(41), 25302-25309. doi:10.1073/pnas.2003914117

Pelkonen, V. P., Niittyvuopio, A., Pirttil\selectlanguage{ngerman}ä, A. M., Laine, K., \&
Hohtola, A. (2007). Phylogenetic background of Orange lily (Lilium
bulbiferum s.l.) cultivars from a genetically isolated
environment.\emph{Plant Biol (Stuttg), 9} (4), 534-540.
doi:10.1055/s-2007-965042

Pelser, P. B., Kennedy, A. H., Tepe, E. J., Shidler, J. B., Nordenstam,
B., Kadereit, J. W., \& Watson, L. E. (2010). Patterns and causes of
incongruence between plastid and nuclear Senecioneae (Asteraceae)
phylogenies. \emph{Am J Bot, 97} (5), 856-873. doi:10.3732/ajb.0900287

Peterson, A. T., Pape, M., \& Soberón., J. (2008). Rethinking receiver
operating characteristic analysis applications in ecological niche
modeling. \emph{Ecol Modell} (1).

Petit, R. J., Duminil, J., Fineschi, S., Hampe, A., Salvini, D., \&
Vendramin, G. G. (2005). Comparative organization of chloroplast,
mitochondrial and nuclear diversity in plant populations. \emph{Mol
Ecol, 14} (3), 689-701. doi:10.1111/j.1365-294X.2004.02410.x

Phillips, S. J., Anderson, R. P., \& Schapire, R. E. (2006). Maximum
entropy modeling of species geographic distributions. \emph{Ecol Modell,
190} (3-4), 231-259.

Pons, O., \& Petit, R. J. (1996). Measuring and testing genetic
differentiation with ordered versus unordered alleles. \emph{Genetics,
144} (3), 1237-1245. doi:10.1093/genetics/144.3.1237

Qian, H., Ricklefs, R. E., \& Thuiller, W. (2021). Evolutionary assembly
of flowering plants into sky islands. \emph{Nat Ecol Evol, 5} (5),
640-646. doi:10.1038/s41559-021-01423-1

Qiu, Y. X., Fu, C. X., \& Comes, H. P. (2011). Plant molecular
phylogeography in China and adjacent regions: Tracing the genetic
imprints of Quaternary climate and environmental change in the world's
most diverse temperate flora. \emph{Mol Phylogenet Evol, 59} (1),
225-244. doi:10.1016/j.ympev.2011.01.012

Qu, X. J., Moore, M. J., Li, D. Z., \& Yi, T. S. (2019). PGA: a software
package for rapid, accurate, and flexible batch annotation of
plastomes.\emph{Plant Methods, 15} (1), 1-12.
doi:10.1186/s13007-019-0435-7

Rambaut, A. (2014). FigTree: a graphical viewer of phylogenetic trees,
version 1.4. 2. \emph{Institute of Evolutionary Biology, Univ.
Edinburgh} .

Ronquist, F., Teslenko, M., van der Mark, P., Ayres, D. L., Darling, A.,
Höhna, S., . . . Huelsenbeck, J. P. (2012). MrBayes 3.2: efficient
Bayesian phylogenetic inference and model choice across a large model
space. \emph{Syst Biol, 61} (3), 539-542. doi:10.1093/sysbio/sys029

Rozas, J., Ferrer-Mata, A., Sánchez-DelBarrio, J. C., Guirao-Rico, S.,
Librado, P., Ramos-Onsins, S. E., \& Sánchez-Gracia, A. (2017). DnaSP 6:
DNA Sequence Polymorphism Analysis of Large Data Sets. \emph{Mol Biol
Evol, 34} (12), 3299-3302. doi:10.1093/molbev/msx248

Saylor, J. E., Quade, J., Dettman, D. L., Decelles, P. G., Kapp, P. A.,
\& Ding, L. (2009). The late Miocene through present paleoelevation
history of southwestern Tibet. \emph{Am J Sci, 309} (1), 1-42.

Schneider, S., \& Excoffier, L. (1999). Estimation of past demographic
parameters from the distribution of pairwise differences when the
mutation rates vary among sites: application to human mitochondrial
DNA.\emph{Genetics, 152} (3), 1079-1089. doi:10.1093/genetics/152.3.1079

Shahzad, K., Jia, Y., Chen, F. L., Zeb, U., \& Li, Z. H. (2017). Effects
of Mountain Uplift and Climatic Oscillations on Phylogeography and
Species Divergence in Four Endangered Notopterygium Herbs. \emph{Front
Plant Sci, 8} , 1929. doi:10.3389/fpls.2017.01929

Soltis, P. S., \& Soltis, D. E. (1991). Genetic Variation in Endemic and
Widespread Plant Species. \emph{Aliso, 13} (1).

Song, G., Zhang, R. Y., Qu, Y. H., Wang, Z. H., Dong, L., Kristin, A., .
. . Lei, F. M. (2016). A zoogeographical boundary between the
Palaearctic and Sino-Japanese realms documented by consistent
north/south phylogeographical divergences in three woodland birds in
eastern China. \emph{J Biogeogr, 43} (11).

Stamatakis, A. (2014). RAxML version 8: a tool for phylogenetic analysis
and post-analysis of large phylogenies. \emph{Bioinformatics, 30} (9),
1312-1313. doi:10.1093/bioinformatics/btu033

Straume, E. O., Nummelin, A., Gaina, C., \& Nisancioglu, K. H. (2022).
Climate transition at the Eocene-Oligocene influenced by bathymetric
changes to the Atlantic-Arctic oceanic gateways. \emph{Proc Natl Acad
Sci U S A, 119} (17), e2115346119. doi:10.1073/pnas.2115346119

Suchard, M. A., Lemey, P., Baele, G., Ayres, D. L., Drummond, A. J., \&
Rambaut, A. (2018). Bayesian phylogenetic and phylodynamic data
integration using BEAST 1.10. \emph{Virus Evol, 4} (1), vey016.
doi:10.1093/ve/vey016

Sun, X., \& Wang, P. (2005). How old is the Asian monsoon system?
Paleobotanical records from China. \emph{Paleogeogr Paleoclimatol
Paleoecol, 222} (3-4), 181-222.

Sun, Y., Abbott, R. J., Li, L., Li, L., Zou, J., \& Liu, J. (2014).
Evolutionary history of Purple cone spruce (Picea purpurea) in the
Qinghai-Tibet Plateau: homoploid hybrid origin and Pleistocene
expansion. \emph{Mol Ecol, 23} (2), 343-359. doi:10.1111/mec.12599

Taberlet, P., Fumagalli, L., Wust-Saucy, A. G., \& Cosson, J. F. (1998).
Comparative phylogeography and postglacial colonization routes in
Europe. \emph{Mol. Ecol, 7} (4), 453-464.

Tajima, F. (1989). Statistical method for testing the neutral mutation
hypothesis by DNA polymorphism. \emph{Genetics, 123} (3), 585-595.
doi:10.1093/genetics/123.3.585

Tu, T., Volis, S., Dillon, M. O., Sun, H., \& Wen, J. (2010). Dispersals
of Hyoscyameae and Mandragoreae (Solanaceae) from the New World to
Eurasia in the early Miocene and their biogeographic diversification
within Eurasia. \emph{Mol Phylogenet Evol, 57} (3), 1226-1237.
doi:10.1016/j.ympev.2010.09.007

Tzedakis, P. C., Lawson, I. T., Frogley, M. R., Hewitt, G. M., \&
Preece, R. C. (2002). Buffered tree population changes in a Quaternary
refugium: evolutionary implications. \emph{Science, 297} (5589),
2044-2047.

Vranckx, G., Jacquemyn, H., Muys, B., \& Honnay, O. (2012).
Meta-analysis of susceptibility of woody plants to loss of genetic
diversity through habitat fragmentation. \emph{Conserv Biol, 26} (2),
228-237. doi:10.1111/j.1523-1739.2011.01778.x

Wang, H., Laqiong, Sun, K., Lu, F., Wang, Y., Song, Z., . . . Zhang, W.
(2010). Phylogeographic structure of Hippophae tibetana (Elaeagnaceae)
highlights the highest microrefugia and the rapid uplift of the
Qinghai-Tibetan Plateau. \emph{Mol Ecol, 19} (14), 2964-2979.
doi:10.1111/j.1365-294X.2010.04729.x

Wang, L., Yang, W., Zhang, X., Zhang, D., \& Zhang, G. (2023). Insights
into the evolutionary history and taxonomic status of Sinopteris
(Pteridaceae). \emph{Mol Phylogenet Evol, 180} , 107672.
doi:10.1016/j.ympev.2022.107672

Wang, L. Y., Ikeda, H., Liu, T. L., Wang, Y. J., \& Liu, J. Q. (2009).
Repeated range expansion and glacial endurance of Potentilla glabra
(Rosaceae) in the Qinghai-Tibetan plateau. \emph{J Integr Plant Biol,
51} (7), 698-706. doi:10.1111/j.1744-7909.2009.00818.x

Wang, X. X., Yue, J. P., Sun, H., \& Li, Z. M. (2011). Phylogeographical
study on Eriophyton wallichii (Labiatae) from alpine scree of
Qinghai-Tibetan Plateau. \emph{Plant Divers. Resour., 33} (6), 605-614.

White, T. J., Bruns, S., Lee, S., \& Taylor, J. (1990). Amplification
and direct sequencing of fungal ribosomal RNA genes for
phylogenetics.\emph{Trends Biochem Sci} (1), 315-322. doi:
10.1016/B978-0-12-372180-8.50042-1

Xia, M., Liu, Y., Liu, J., Chen, D., Shi, Y., Chen, Z., . . . Qiu, Y.
(2022). Out of the Himalaya-Hengduan Mountains: Phylogenomics,
biogeography and diversification of Polygonatum Mill. (Asparagaceae) in
the Northern Hemisphere. \emph{Mol Phylogenet Evol, 169} , 107431.
doi:10.1016/j.ympev.2022.107431

Xie, C., Xie, D. F., Zhong, Y., Guo, X. L., Liu, Q., Zhou, S. D., \& He,
X. J. (2019). The effect of Hengduan Mountains Region (HMR) uplift to
environmental changes in the HMR and its eastern adjacent area: Tracing
the evolutionary history of \emph{Allium} section
\emph{Sikkimensia}(Amaryllidaceae). \emph{Mol Phylogenet Evol, 130} ,
380-396. doi:10.1016/j.ympev.2018.09.011

Xie, D. F., Tan, J. B., Yu, Y., Gui, L. J., Su, D. M., Zhou, S. D., \&
He, X. J. (2020). Insights into phylogeny, age and evolution
of\emph{Allium} (Amaryllidaceae) based on the whole plastome
sequences.\emph{Ann Bot, 125} (7), 1039-1055. doi:10.1093/aob/mcaa024

Xiong, Z., Liu, X., Ding, L., Farnsworth, A., Spicer, R. A., Xu, Q., . .
. Yue, Y. (2022). The rise and demise of the Paleogene Central Tibetan
Valley. \emph{Sci Adv, 8} (6), eabj0944. doi:10.1126/sciadv.abj0944

Xu, T., Abbott, R. J., Milne, R. I., Mao, K., Du, F. K., Wu, G., . . .
Liu, J. (2010). Phylogeography and allopatric divergence of cypress
species (Cupressus L.) in the Qinghai-Tibetan Plateau and adjacent
regions. \emph{BMC Evol Biol, 10} , 194. doi:10.1186/1471-2148-10-194

Yang, F. S., Qin, A. L., Li, Y. F., \& Wang, X. Q. (2012). Great genetic
differentiation among populations of Meconopsis integrifolia and its
implication for plant speciation in the Qinghai-Tibetan
Plateau.\emph{PLoS One, 7} (5), e37196. doi:10.1371/journal.pone.0037196

Yang, J., Zhou, S., Huang, D., \& He, X. (2018). Phylogeography of two
closely related species of Allium endemic to East Asia: Population
evolution in response to climate oscillations. \emph{Ecol Evol, 8} (16),
7986-7999. doi:10.1002/ece3.4338

Yang, M., Zhou, S., Xingjin, H. E., \& Peng, C. (2016). Genetic
Relationships of \emph{Notholirion} and Genetic Diversity of
N.bulbuliferum in China Revealed by ISSR Markers. \emph{Acta Botanica
Boreali-Occidentalia Sinica} .

Yu, X. Q., Maki, M., Drew, B. T., Paton, A. J., Li, H. W., Zhao, J. L.,
. . . Li, J. (2014). Phylogeny and historical biogeography of Isodon
(Lamiaceae): rapid radiation in south-west China and Miocene overland
dispersal into Africa. \emph{Mol Phylogenet Evol, 77} , 183-194.
doi:10.1016/j.ympev.2014.04.017

Yu, Y., Blair, C., \& He, X. (2020). RASP 4: Ancestral State
Reconstruction Tool for Multiple Genes and Characters. \emph{Mol Biol
Evol, 37} (2), 604-606. doi:10.1093/molbev/msz257

Yu, Y., Harris, A. J., Blair, C., \& He, X. (2015). RASP (Reconstruct
Ancestral State in Phylogenies): a tool for historical
biogeography.\emph{Mol Phylogenet Evol, 87} , 46-49.
doi:10.1016/j.ympev.2015.03.008

Zhang, D., Gao, F., Jakovli\selectlanguage{polish}ć, \selectlanguage{english}I., Zou, H., Zhang, J., Li, W. X., \&
Wang, G. T. (2020). PhyloSuite: An integrated and scalable desktop
platform for streamlined molecular sequence data management and
evolutionary phylogenetics studies. \emph{Mol Ecol Resour, 20} (1),
348-355. doi:10.1111/1755-0998.13096

Zhang, P. Z., Shen, Z., Min, W., Gan, W., B\selectlanguage{ngerman}ürgmann, R., Molnar, P., . .
. Wu, J. (2004). Continuous deformation of the Tibetan Plateau from
global positioning system data. \emph{Geology, 32} (9), 809-812.

Zhang, Q., \& Willems, H. (2012). Initial India-Asia Continental
Collision and Foreland Basin Evolution in the Tethyan Himalaya of Tibet:
Evidence from Stratigraphy and Paleontology. \emph{J Geol, 120} (3),
347-347.

Zhang, Y. L., Li, B. Y., \& Zheng, D. (2002). A discussion on the
boundary and area of the Tibetan Plateau in China. \emph{Geogr Res, 21}
(001), 1-8.

Zhang, Z., \& Sun, J. (2011). Palynological evidence for Neogene
environmental change in the foreland basin of the southern Tianshan
range, northwestern China. \emph{Glob Planet Change, 75} (1-2), 56-66.

Zheng, H. B., Powell, C. M., An, Z. S., Zhou, J., \& Dong, G. R. (2000).
Pliocene uplift of the northern Tibetan Plateau. \emph{Geology, 28} (8),
715.

Zheng, H. Y., Guo, X. L., Price, M., He, X. J., \& Zhou, S. D. (2021).
Effects of Mountain Uplift and Climatic Oscillations on Phylogeography
and Species Divergence of Chamaesium (Apiaceae). \emph{Front Plant Sci,
12} , 673200. doi:10.3389/fpls.2021.673200

\textbf{Figure legends}

\textbf{Fig. 1 Sampling locations of \emph{Notholirion}.}

\textbf{Fig. 2 Geographic distribution of \emph{ITS} haplotypes
for\emph{Notholirion}.} Circles represent populations and the colored
outlines of the distinguish different species. Each color represents a
haplotype.

\textbf{Fig. 3 Haplotypes network based on \emph{ITS}
within\emph{Notholirion}.} Each color represents a specie. The size of
circles in the network corresponds to the frequency of each haplotype.
The horizontal bars on the horizontal lines between different haplotypes
represent the number of variant sites. Small solid black circles denoted
hypothetic unsampled or extinct ancestral haplotypes.

\textbf{Fig. 4 Geographic distribution of cpDNA haplotypes
for\emph{Notholirion}.} Each circle represents a

population and the colored outlines of the distinguish different
species. Each color represents a

haplotype.

\textbf{Fig. 5 Haplotypes network based on cpDNA
within\emph{Notholirion}.} Each color represents a specie.

The size of circles in the network corresponds to the frequency of each
haplotype. The horizontal bars on the horizontal lines between different
haplotypes represent the number of variant sites. Small solid black
circles denoted hypothetic unsampled or extinct ancestral haplotypes.

\textbf{Fig. 6 Phylogenetic relationship of \emph{Notholirion} recovered
from the (a) \emph{ITS} haplotypes and (b) chloroplast haplotypes.} Bule
clades: \emph{N. bulbuliferum} , yellow clades: \emph{N. macrophyllum} ,
red purple clades: \emph{N. thomsonianum} . Maximum-likelihood bootstrap
support/Bayesian posterior probability values are shown near
corresponding nodes (`\selectlanguage{english}-' indicates support values less than 50\%;
`\textbf{*} ' represents 100\%/1 support).

\textbf{Fig. 7 Chronogram of \emph{Notholirion} from the phylogenetic
relationship of Liliales based on}

\textbf{plastomes.} The red boxes show species of the\emph{Notholirion}
.

\textbf{Fig. 8 Divergence time estimation based on plastomes
of\emph{Notholirion}.} Blue bars and numbers above represent 95\%
highest posterior density (95\% HPDs) for each node.

\textbf{Fig. 9 Ancestral area reconstructions based on the cpDNA
phylogeny within \emph{Notholirion}.} Pie charts show proportions of the
ancestral ranges. The map shows the three geographical distribution
areas of \emph{Notholirion} : \textbf{(A)} Southern
Himalayas,\textbf{(B)} Eastern Himalayas, \textbf{(C)} Hengduan
Mountains. The divergence time (in Ma) based on cpDNA. Numbers in the
brackets show the 95\% HPD of divergence time (in Ma) of the main nodes.

\textbf{Fig. 10 Mismacth distribution analysis for \emph{ITS} data (a)
and cpDNA data (b).} The line represents the distributions of an
expected population expansion, the dashed line show observed (Obs)
values.

\textbf{Fig. 11 Potential distribution ranges of different periods
for\emph{Notholirion} were simulated by Species Distribution Models
using Bioclimatic variables.} Colors represent bioclimatic suitability,
from most suitable (red) to unsuitable (gray). LIG: Last Interglacial,
LGM: Last Glacial Maximum.

\textbf{Supporting information}

\textbf{Table S1} The information of sample collection and geographical
area for \emph{Notholirion} .

\textbf{Table S2} The genbank accessions for species used in
differentiation time estimates.

\textbf{Table S3} Priors and mean ages used to calibrate the Liliales
phylogenomic tree.

\textbf{Table S4} The haplotype distribution, haplotype diversity (Hd)
and nucleotide diversity (Pi) in each population
within\emph{Notholirion} .\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Figure 1: Sampling locations of \emph{Notholirion}}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Figure 2: Geographic distribution of ITS haplotypes
for\emph{Notholirion}.}

Circles represent populations and the colored outlines of the
distinguish different species. Each color represents a haplotype.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Figure 3: Haplotypes network based on ITS
within\emph{Notholirion}.}

Each color represents a specie. The size of circles in the network
corresponds to the frequency of each haplotype. The horizontal bars on
the horizontal lines between different haplotypes represent the number
of variant sites. Small solid black circles denoted hypothetic unsampled
or extinct ancestral haplotypes.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{Figure 4: Geographic distribution of cpDNA haplotypes
for\emph{Notholirion}.}

Each circle represents a population and the colored outlines of the
distinguish different species. Each color represents a haplotype.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\textbf{Figure 5: Haplotypes network based on cpDNA
within\emph{Notholirion}.}

Each color represents a specie. The size of circles in the network
corresponds to the frequency of each haplotype. The horizontal bars on
the horizontal lines between different haplotypes represent the number
of variant sites. Small solid black circles denoted hypothetic unsampled
or extinct ancestral haplotypes.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

\textbf{Figure 6: Phylogenetic relationship of
\emph{Notholirion}recovered from the (A) ITS haplotypes and (B)
chloroplast haplotypes.}

Bule clades: \emph{N. bulbuliferum} , yellow clades: \emph{N.
macrophyllum} , red purple clades: \emph{N. thomsonianum} .
Maximum-likelihood bootstrap support/Bayesian posterior probability
values are shown near corresponding nodes (`-' indicates support values
less than 50\%; `*' represents 100\%/1 support).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

\textbf{Figure 7: Chronogram of \emph{Notholirion} from the phylogenetic
relationship of Liliales based on}

\textbf{plastomes.}

The red boxes show species of the \emph{Notholirion} .\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image8/image8}

\end{center}
\end{figure}

\textbf{Figure 8: Divergence time estimation based on plastomes
of\emph{Notholirion}.}

Blue bars and numbers above represent 95\% highest posterior density
(95\% HPDs) for each node.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image9/image9}

\end{center}
\end{figure}

\textbf{Figure 9: Ancestral area reconstructions based on the cpDNA
phylogeny within \emph{Notholirion}.}

Pie charts show proportions of the ancestral ranges. The map shows the
three geographical distribution areas of \emph{Notholirion} : (A)
Southern Himalayas, (B) Eastern Himalayas, (C) Hengduan Mountains. The
divergence time (in Ma) based on cpDNA. Numbers in the brackets show the
95\% HPD of divergence time (in Ma) of the main nodes.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image10/image10}

\end{center}
\end{figure}

\textbf{Figure 10: Mismacth distribution analysis for \emph{ITS} data
(A) and cpDNA data (B).}

The line represents the distributions of an expected population
expansion, the dashed line show observed (Obs) values.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image11/image11}

\end{center}
\end{figure}

\textbf{Figure 11: Potential distribution ranges of different periods
for \emph{Notholirion} were simulated by Species Distribution Models
using Bioclimatic variables.}

Colors represent bioclimatic suitability, from most suitable (red) to
unsuitable (gray). LIG: Last Interglacial, LGM: Last Glacial Maximum.

\textbf{Hosted file}

\verb`Supplementary table1.docx`  available at \url{https://authorea.com/users/753582/articles/723781-effects-of-mountain-uplift-and-climatic-oscillations-on-phylogeography-and-species-divergence-of-notholirion-liliaceae}

\textbf{Hosted file}

\verb`Supplementary table2.docx`  available at \url{https://authorea.com/users/753582/articles/723781-effects-of-mountain-uplift-and-climatic-oscillations-on-phylogeography-and-species-divergence-of-notholirion-liliaceae}

\textbf{Hosted file}

\verb`Supplementary table3.docx`  available at \url{https://authorea.com/users/753582/articles/723781-effects-of-mountain-uplift-and-climatic-oscillations-on-phylogeography-and-species-divergence-of-notholirion-liliaceae}

\textbf{Hosted file}

\verb`Supplementary table4.docx`  available at \url{https://authorea.com/users/753582/articles/723781-effects-of-mountain-uplift-and-climatic-oscillations-on-phylogeography-and-species-divergence-of-notholirion-liliaceae}

\selectlanguage{english}
\FloatBarrier
\end{document}

