Phylogenetics and Biogeography of Lilium ledebourii from the Hyrcanian Forest

: Lilium ledebourii (Baker) Boiss is one of the most endangered lilies, restricted to only a few small and fragmented areas in the Hyrcanian forest. This study aimed at evaluating the taxonomy of this unique Iranian lily and reconstructing divergence time from other species of the genus Lilium to address the role of this region in its diversiﬁcation. Phylogenetic trees based on nuclear ITS and chloroplastic matK strongly supported the monophyly of the genus Lilium and division into subclades hardly matching prior morphological classiﬁcations. Biogeographic analyses using S-DIVA revealed East Asia as the ancestral range from where Lilium presented a multidirectional expansion towards North America, West-Central Asia, North Asia, and Europe. Diverging from ancestral Lilium during the beginning of Eocene (50 Ma; 95% HDP: 68.8–36.8). Speciﬁc members of Lilium colonized Iran (Western Asia) separated from the Clade IV (West-Central Asia and Europe lineage), and then yielded the Iranian L. ledebourri . Accordingly, the north of Iran appears to have promoted both long-term persistence and migration of Lily species from Asia to the Europe.


Introduction
Lilium L. (Liliaceae) is one of the most fascinating genera among the plant kingdom, and it adds to its importance on horticulture, medicine and food [1][2][3].Lilium is a perennial herbs with subterranean bulbs and mostly spring-flowering that grow in steppes and mountain meadows [4].This genus occurs in Eurasia and North America and has three main ranges including the Qinghai-Tibetan Plateau in East Asia, North America, and the Caucasus [5,6].The genus of Lilium, with approximately 100 species, has been classified into five to eleven sections [7][8][9].Due to the frequent gene flow among sections, the major phylogenetic clades of Lilium are still controversial but have been clear [2,3,8,10].Traditionally, seven sections were recognized based on morphological taxonomy and molecular phylogenetic methods [2,8,10].
The Hyrcanian area in the southern shores of the Caspian Sea presents remnants of natural deciduous forests [11] and forms a unique vegetation belt from the Talish region in Azerbaijan to Golestan National Park in Iran (between 48 • -56 • E and 38 • 55 -35 • 05 N) [12].Decreasing water levels during the Quaternary, coupled with the rising of the Earth's crust in the northern Caucasus, led to the separation of the Black Sea and the Caspian Sea, promoting the isolation of the Hyrcanian region [13].The presence of many Arcto-Tertiary relict species in the area was generally viewed as coherent with the Hyrcanian forest having acted as major refugia sheltering the long-term persistence of plant species [14,15].To confirm this hypothesis, additional data are necessary.
Lilium ledebourii (Baker) Boiss (Iranian lily) is a rare and endemic lily species distributed only in a few small areas of the Hyrcanian forest [16].Its habitat was recorded as a natural monument in Iran's protected areas in 1976 [17].This species is attractive and widely used in ornamental breeding programs, and some biochemical features were also studied [18,19].Limited information is available regarding its phylogenetic position; however, there was one recent study based on the nuclear rDNA-internal transcribed spacers (ITS) and another based on the 5.8 S ribosomal DNA sequence [19,20].However, due to the incomplete sampling on the distribution range of L. ledebourii and single sequence, its phylogenetic position is still unclear.Furthermore, there are several biogeographic studies of the genus Lilium missing this species.As the famous species is restricted to the biodiversity hotspots, L. ledebourii presents a less known biogeographic synthesis, and deserves further attention [21,22].
Despite the advent of genomics, DNA barcoding based on nucleotide sequences from variable chloroplastic and nuclear loci is still facilitated to monitor the phylogenetic position accurately and economically [23][24][25].Here, we thus address the phylogenetic position of L. ledebourii concerning worldwide Lilium species and shed light on its role for the biogeography and diversification of the genus.The first purpose of this study is to present the phylogenetic relationships of the Iranian Lily based on nuclear and plastid DNA markers; the second goal is to estimate the divergence time of the Iranian Lily and discuss the biogeography of Lilium after adding the Iranian Lily.

Sampling and DNA Extraction
Leaf samples of L. ledebourii were collected from four natural populations (all reported sites for this species) of the Hyrcanian forest (Figure 1).Total genomic DNA was extracted from at least 10 individuals from each population following Murray and Thompson (1980), with modifications according to Janfaza et al. [26].

PCR Amplification and Sequencing
For amplification and sequencing, four candidate DNA barcoding regions were selected: the internal transcribed spacer regions (ITS) from nuclear genomes and three plastid regions (matK, trnL-F and psbA-trnH).PCR amplifications were performed in 20 µL using the AccuPower HotStart PCR PreMix kit (Bioneer, Daejeon, Korea), with either universal ITS-1 and ITS-4 primers from White et al. [27] or matK primers from Johnson et al. [28] or primers for the trnL-F intergenic spacer from Taberlet et al. [29] or primers for the psbA-trnH intergenic spacer from Sang et al. [30].For ITS, psbA-trnH and trnL-F regions, PCR cycles consisted of an initial denaturation for 6 min at 95 • C, followed by 32 cycles of 60 s at 95 • C, 45 s at 56 • C, 90 s at 72 • C, and a final extension of 5 min 72 • C. The matK region was amplified following Hayashi and Kawano (2000) [31], including 35 cycles of 60 s at 94 • C, 120 s at 50 • C, 180 s at 72 • C and a final extension at for 7 min at 72 • C).After a check through electrophoresis, PCR products were sent to Bioneer service (Bioneer, Daejeon, Korea) for sequencing, and the resulting electropherograms were manually inspected using Chromas version 2.3 (www.technelysium.com.au(15 February 2020) to yield accurate sequence data.

PCR Amplification and Sequencing
For amplification and sequencing, four candidate DNA barcoding regions were selected: the internal transcribed spacer regions (ITS) from nuclear genomes and three plastid regions (matK, trnL-F and psbA-trnH).PCR amplifications were performed in 20 μL using the AccuPower HotStart PCR PreMix kit (Bioneer, Daejeon, Korea), with either universal ITS-1 and ITS-4 primers from White et al. [27] or matK primers from Johnson et al. [28] or primers for the trnL-F intergenic spacer from Taberlet et al. [29] or primers for the psbA-trnH intergenic spacer from Sang et al. [30].For ITS, psbA-trnH and trnL-F regions, PCR cycles consisted of an initial denaturation for 6 min at 95 °C, followed by 32 cycles of

Phylogenetic Analysis and ITS2 Secondary Structure
There are eight individuals of L. ledebourii were used in the phylogenetic analysis in this study.The analysis based on ITS sequences used 35 Lilium species and took five Tulipa species as outgroups.The analysis based on matK sequences used 21 Lilium species and taken two Tulipa species as outgroups.
Alignments of sequences were performed for each locus using MUSCLE and nucleotide composition, the number of variables, parsimony and conserved sites among studied taxa were estimated in MEGA 6 [32].A median-joining (MJ) network was inferred from all variable characters of the complete alignment of concatenated sequences using NETWORK 3.1.1.1 [33].The phylogenetic tree was constructed based on maximum likelihood (ML) and Bayesian inference (BI) for all taxa.
Bayesian inference tree was constructed by MrBayes v. 3.2.6 with a cold chain and three incrementally heated chains (T = 0.2), running for 10,000,000 generations with a sampling frequency of 1000.The first 25% of the trees were discarded (burnin) and the remaining trees were used to build a consensus tree and estimate Bayesian posterior probabilities.The MEGA (version 6) [32] software was used to choose the best model and create phylogenetic trees using maximum likelihood (ML) methods, [34].Bootstrap analysis with 1000 replications was used to analyze the reliability of each branch.ITS2 secondary structure for each taxon was inferred by minimizing free energy and homology search across the ITS2 database (http://its2.bioapps.biozentrum.uni-wuerzburg.de/(27 February 2020).

Divergence Time Estimate and Biogeographic Analysis
To estimate the divergence time of Lilium and reconstruct the ancestral geographical range of this genus, we used 31 Lilium species covering its entire distribution area.Partition homogeneity or incongruence length difference (ILD) [35] test was conducted in PAUP v4.0b10 [36] following a heuristic search approach of 1000 replicates with 100 random stepwise additions and tree bisection reconstruction (TBR) branch swapping supported congruity of ITS and matK datasets (p = 0.68).Similarly, the likelihood ratio test (LRT) carried out in PAUP v4.0b10 [36] supported a relaxed molecular clock approach for our combined dataset (p = 0.00001).
Divergence times were estimated by calibrating the stem node of Ripogonaceae, Luzuriaga, according to Kim and Kim (2018) [37], and the crown node of Smilax (the outgroup) to a mean age of 55 million years ago (Ma), 24 Ma and 46 Ma respectively, using the log-normal distribution prior with standard deviation of 1.0.The phylogenetic analysis was conducted for 50,000,000 generations using BEAST ver.1.6.1 [38] under uncorrelated lognormal relaxed clock parameter and the Yule model of speciation as tree prior.The GTR + G substitution model was selected as the best fit model using jModel test under Akaike's information criterion (AIC).Tree Annotatorver.1.7.5 [39] and FigTreever.1.4[40] were used for the generation and visualization of the maximum clade credibility tree.
The distribution of Lilium was divided into six regions based on the available samples: A (East Asia), B (North Asia), C (West and Central Asia), D (Europe), E (North America), and F (South America).Reconstruct Ancestral State in Phylogenies (RASP) software [41] was used for biogeographic inferences using Statistical Dispersal-Vicariance Analysis (S-DIVA) approaches.To overcome uncertainties associated with S-DIVA analysis out of 10,000 trees from BEAST, the maximum clade credibility tree and distribution file were uploaded into the RASP software for biogeographic reconstructions.In the BBM analysis, only the maximum clade credibility tree was used along with the distribution file to obtain a reconstruction and the Markov chain Monte Carlo (MCMC) was run under the JC + G (Jukes-Cantor + Gamma) model for 5,000,000 generations.

ITS and Plastid Regions (MatK, TrnL-F and TrnH-PsbA)
The nucleotide composition of the Iranian samples of L. ledebourii from four different regions is summarized and compared with available samples of lilies in Table 1.ITS sequences from Iranian samples were longer (644-654 bp) than those retrieved from NCBI (465-485 bp).The matK region showed low variability among Lilium species (length = 853-881 bp, conserved sites = 849).Iranian samples presented 19 variable sites, including nine singletons, whereas species available in NCBI presented 36 variable sites.The trnL-F and trnH-psbA regions showed the high conservation sites (length = 270-255, conserved sites = 246) and nucleotide percent of Temin (Iranian sample = 34.5, sample from NCBI = 35.5).

Phylogenetic Position of L. ledebourii
Bayesian tree and ML trees based on ITS strongly supported a monophyletic genus Lilium and divided the species into six clades (Figure 2A), consistent with the MJ network (Supplementary Material Figure S1).These clades are somewhat consistent with the four main groups defined by Baker et al. [42]

Biogeography of the Genus Lilium
The results of S-DIVA analyses suggest that the current distribution of Lilium species is the outcome of multiple dispersal and vicariance events (Figure 3).S-DIVA suggests 9 dispersal and three vicariance events.They congruently indicate that Lilium ancestors (node I) originated in East Asia, or East-West-Central Asia (Figure 3) at the beginning of Eocene (c.a.50 MYA, 95% HDP: 68.8-36.8,Table 2).Then, the ancestral Lilium from East Asia began to diverge into two clades (Clade II and Clade III) during the early and middle Eocene (Figure 3).During the Late Eocene and Oligocene, the Lilium began to spread to The Bayesian matK tree has a backbone polytomy and four major lineages, with significant differences in composition between these lineages and the six resolved in the ITS phylogeny.The phylogenetic tree based on this region showed a high inconsistency with the Baker et al., 1871 classification, but divided all Hyrcanian specimens along with L. ciliatum from the Caucasus and L. pyrnaicum from Europe into one sub clade.

Biogeography of the Genus Lilium
The results of S-DIVA analyses suggest that the current distribution of Lilium species is the outcome of multiple dispersal and vicariance events (Figure 3).S-DIVA suggests 9 dispersal and three vicariance events.They congruently indicate that Lilium ancestors (node I) originated in East Asia, or East-West-Central Asia (Figure 3) at the beginning of Eocene (c.a.50 MYA, 95% HDP: 68.8-36.8,Table 2).Then, the ancestral Lilium from East Asia began to diverge into two clades (Clade II and Clade III) during the early and middle Eocene (Figure 3).During the Late Eocene and Oligocene, the Lilium began to spread to North America (Clade II) and West and Central Asia (Clade III).The Clade II have been included two obvious lineages (East Asia and North America lineage), and they carried out an independent evolution and diffusion from Miocene.In Clade III, the Lilium species in East Asia have four times the dispersal after Pliocene.Ancestral reconstructions of Clades IV and V further favor an ancestral range in western Asia and Europe (with 93.30% and 78% marginal probability, respectively) in Clade III, including the Iranian L. ledebourri.Lilium ledebourii in West Asia looks like the first lineage separate from Clade IV.The overall scenario favored by S-DIVA (Table 3) supports stepping-stone radiation of the genus Lilium, with a primary center of diversification in eastern Asia (i.e., region A).Our inclusion of the Iranian samples here looks critical, as it supports an expansion and secondary diversification in western Asia (i.e., region C) before further range expansion and diversification towards Europe (i.e., region D).The overall scenario favored by S-DIVA (Table 3) supports stepping-stone radiation of the genus Lilium, with a primary center of diversification in eastern Asia (i.e., region A).Our inclusion of the Iranian samples here looks critical, as it supports an expansion and secondary diversification in western Asia (i.e., region C) before further range expansion and diversification towards Europe (i.e., region D).

Phylogeny and Biogeography of Lilies
Our research sheds light on the natural divergence histories of L. ledebourii in the Caucasian Hyrcanian forest of Iran, which is a hotspot for biodiversity.The results showed that L. ledebourii from this region belonging to the Caucasians are older than the European samples, and thus Iran could be a bridge for transforming L. ledebourii into Europe.Furthermore, we discovered that the phylogenetic position of the Lilium genus did not follow the flower shape.
Phylogenetic inferences based on ITS and matK loci here largely confirmed earliest classifications based on morphology to distinguish either five sections [9] or four main groups [42].Iranian lilies (L.ledebourii), analyzed here for the first time, were closely related to European species (L.pyrenaicum) as well as species from Caucasus [43] to be part of the Martagon section composed of species forming mostly Turk's cap flowers.
The inclusion of the Iranian samples offers a comprehensive sampling with representative species from most regions colonized by lilies to bring fresh insights into their evolution.For instance, the controversial L. lophophorum here appears consistent with the classification of Wilson [44] in a Lophophorum subsection, whereas close relationships with L. humboldtii further supports an early origin out of Chinese lineages.All species comprised here in clade V (sections Archelirion and Martagon of Baker's classification) share an origin in east Asia and support Comber [8] classification.Although the impact of hybridization remains controversial in the genus Lilium [45], present results indicate close relationships between species of Sinomartagon and Archelirion sections, supporting the view of Leslie [46] about hybridization between corresponding species.In particular, L. henrici here is shown to be very close to L. alexandrae and these species may have undergone gene flow.However, inference of reticulate evolution across the genus Lilium is beyond the scope of this work.

Biogeography of Lilium Emphasis on West Asia
Even though Lilium is one of the most famous genera all over the world, the biogeographic study of this genus is limited [5].Gao et al. (2013) tried to analyze the evolutionary events in Lilium.Due to the sampling, this study only focused on the biogeography of the Q-T plateau and the Hengduan Mountains.In this study, we try to cover the distribution area of the genus to explore the biogeography of this genus, adding some species from West Asia.
Biogeographic reconstructions using S-DIVA indicated an ancestral origin of Lilium in Eastern Asia or the extensive area of East-Central-West Asia during the Eocene around 50 MYA (95% HDP: 68.8-36.8).This contrasts with prior dating based on the plastid dataset that concluded on a last common ancestor of Lilium some 13.19 MYA [5].Thus, our results suggest that Lilium may have early origins in the Pan-Asian region.The genus began to spread outside of Asia during the Eocene.Two early radiation events from the ancestral area have particularly supported independent colonization toward Western Asia and towards North America.This event may have benefited from the warm climatic conditions in this Geologic period.This is in keeping with prior hypotheses of migration for diverse plants and animals between the Old and the New Worlds through either the North Atlantic Land Bridge (NALB) or Beringia [47], rather than the suggestion in some works of literature that the Himalayas were the center of origin of the genus Lilium [5,6].This pattern also contrasts with the proposal of Ikinci et al. [48] that European species derive from multiple routes, with L. martagon first colonizing Europe, whereas a second route later gave rise to L. bulbiferum and further diversification of all other European species.Lilium species of Iran (Western Asia) indeed involved radiation of ancestors from East Asia to Europe, with early colonization of Western Asia from East Asia (50-40 MYA) followed by additional events some 20 MYA and 10 MYA.In particular, three radiations were apparent between Western Asia and Europe in our analysis, supporting dispersal from Western Asia to Europe during the Oligocene (28 MYA) and again during the Miocene (18 MYA).
Due to the limited impact of the ice ages during the Pleistocene, it indeed hosts an impressive number of endemic and relict species [49].Our analysis showed that the rare Lilium ledebourii is such a relict species that currently scattered among small populations in cliffs above the tree line.These populations seemingly represent survivors of previously large populations having migrated to the area before the ice age.Accordingly, the northern forests of Iran likely acted as refugia for such species.The presence of tree species from typically high latitudes such as Betula pendula and Sorbus aucaparia in the Hyrcanian forest further corroborates the case of L. ledebourii.

Figure 1 .
Figure 1.Geographical sampling sites of Lilium ledebourii in the Hyrcanian forest, part of Iran.

Figure 1 .
Figure 1.Geographical sampling sites of Lilium ledebourii in the Hyrcanian forest, part of Iran.
based on flower morphology: i.e., L. sect.Eulirion Rchb.(funnel-flowered lilies), L. sect.Archelirion Baker (open-flowered lilies), L. sect.Isolirion Baker (erect-flowered lilies) and L. sect.Martagon Rchb.(Turk's cap lilies).The results of the phylogenetic tree showed that the Caucasian and European species were placed in three separate clade (Figure 2A).All species of the European region were grouped in one clade, while seven species of the Caucasus were grouped in three separate clades.Surprisingly, L. ledebourii was placed in a clade with European species, and the L. ciliatum itself formed a separate clade, while other Caucasian species formed the third clade.Diversity 2022, 14, x FOR PEER REVIEW 6 of 11 the Baker et al., 1871 classification, but divided all Hyrcanian specimens along with L. ciliatum from the Caucasus and L. pyrnaicum from Europe into one sub clade.

Figure 2 .
Figure 2. (A) Phylogenetic tree of Lilium based on ITS sequences (ITS1, 5.8 s, ITS2) with Maximum Likelihood method (model: Tumara3 parameter + Gamma distributed) and bootstrap 1000, (B) Phylogenetic tree of Lilium based on matK sequences with Maximum Likelihood method (model: Tumara3 parameter) and bootstrap 1000.(T and K legends are represented repeated sequences of ITS and matK regions from the Hyrcanian forest, respectively).(C) Four morphological flower types of the genus Lilium based on classification of Baker et al., 1871).

Figure 2 .
Figure 2. (A) Phylogenetic tree of Lilium based on ITS sequences (ITS1, 5.8 s, ITS2) with Maximum Likelihood method (model: Tumara3 parameter + Gamma distributed) and bootstrap 1000, (B) Phylogenetic tree of Lilium based on matK sequences with Maximum Likelihood method (model: Tumara3 parameter) and bootstrap 1000.(T and K legends are represented repeated sequences of ITS and matK regions from the Hyrcanian forest, respectively).(C) Four morphological flower types of the genus Lilium based on classification of Baker et al., 1871).

Figure 3 .
Figure 3. Divergence time estimations and the major divergence events among the Lilium species based on region combined data (ITS + matK regions) and Dispersal scenario of Lilium in the world.Pie charts depict ancestral area reconstruction probability, with the colors of pie slices defined in the legend.(The reader is directed to the web version of this article for an explanation of the color references in this figure legend).

Figure 3 .
Figure 3. Divergence time estimations and the major divergence events among the Lilium species based on region combined data (ITS + matK regions) and Dispersal scenario of Lilium in the world.Pie charts depict ancestral area reconstruction probability, with the colors of pie slices defined in the legend.(The reader is directed to the web version of this article for an explanation of the color references in this figure legend).

Table 1 .
Sequences characterization and nucleotide composition four regions trnH-psbA, trnL-F, matK and ITS without the outgroups.

Table 2 .
Age estimation in million year ago (MYA) of nodes of divergence time clade based on S-DIVA and BBM.

Table 3 .
Dispersal details of different distribution area based on S-DIVA and BBM.