1. Introduction

Species are the fundamental units of biodiversity, playing a crucial role in conservation biology and ecological research (Wheeler and Meier, 2000; Coyne and Orr, 2004; Hong, 2016). However, the concept of species is complex and often controversial. Defining a species has been, and will continue to be, a major challenge for biological studies, affecting the formulation of conservation strategies (Akcakaya et al., 2004). Despite numerous proposed species concepts, no single definition satisfies the entire spectrum of life forms, from unicellular organisms to complex multicellular beings like humans (De Queiroz, 2007; Wilkins, 2009). Therefore, different organisms may require distinct species concepts.

From an evolutionary perspective, species are dynamic entities, continuously evolving through time (Coyne and Orr, 2004; De Queiroz, 2007; Liu, 2016). While complete reproductive isolation has long been considered a key factor in speciation and maintaining species boundaries, particularly in the Biological Species Concept (Mayr, 1963), recent studies have revealed that gene flow is more pervasive than previously thought (Bolnick and Fitzpatrick, 2007; Nosil, 2008; Blanckaert et al., 2020; Wang et al., 2020). This is especially true in plants, where hybridization and introgression occur not only between closely related taxa but also between more distantly related congeners or even across genera (Hu et al., 2021; Wu et al., 2022; Zhou et al., 2022; Mason-Gamer and White, 2024). Surprisingly, despite this gene flow, many species maintain distinct morphological and ecological identities. This paradox raises intriguing questions about the mechanisms of species persistence and boundary maintenance in the face of ongoing genetic exchange.

Contemporary biologists have proposed more integrative approaches to species delimitation. These incorporate multiple lines of evidence, including breeding systems, geography, morphology, DNA sequences, and ecology (Liu, 2016; Donoghue, 2022). Moreover, there is a growing recognition that speciation is a continuous process, with extant species representing various stages of divergence and independence (Liu, 2016; Gao et al., 2019; Zhang et al., 2023). This perspective acknowledges that species boundaries are dynamic and constantly challenged by gene flow, with different lineages adopting unique strategies to maintain their ecological and genetic distinctiveness.

The genus Lilium, spanning the temperate Northern Hemisphere, is renowned for its showy flowers and horticultural importance. Its success in cultivation stems from extensive interspecific hybridization, indicating weak reproductive barriers among species. Lilium species are typically self-incompatible and pollinator-dependent, making them prone to hybridization in both natural and cultivated settings (Marasek-Ciolakowska et al., 2018). Intriguingly, our observations indicate that while hybridization is common in Lilium species, species boundaries remain relatively stable in nature. In this genus, morphological diversity of flowers has led to various sectional classifications (Comber, 1949; Nishikawa et al., 1999; Hayashi and Kawano, 2000). For example, four perigone types correlate with specific pollinators in North America (Givnish et al., 2020), while in East Asia, the diversity is even greater. However, traditional classifications, such as assigning trumpet and recurved forms to sections Leucolirion and Sinomartagon respectively (Comber, 1949), have been challenged by recent phylogenetic studies (Kim et al., 2019; Duan et al., 2022; Li et al., 2022). Our previous research contradicts studies that suggested floral morphology reflects genetic relationships, instead indicating that floral morphology results from parallel evolution in response to similar habitats (Gao et al., 2015). This phenomenon is widely observed throughout the entire genus, making sub-generic classification much more difficult (Liu et al., 2019; Yuan and Gao, 2024; Zhou et al., 2024).

The monophyletic clade Leucolirion proposed by Gao et al. (2013) comprises six taxa exhibiting two distinct floral forms as mentioned above. This clade is endemic to central and southwestern China, with Lilium sulphureum Baker ex Hook. f. extending into neighboring Myanmar. Four members represented by L. regale E.H. Wilson—discovered and introduced by the famous plant hunter E.H. Wilson in the Min River valley—are characterized by the white giant trumpet flower with a strong fragrance. In contrast, the other two taxa (L. henryi Baker and L. rosthornii Diels), from central China, exhibit strongly recurved, inodorous orange to red petals adorned with black speckles. Despite their strikingly different appearances, these taxa share unique purplish bulbs characteristic of Chinese lilies. Our previous research (Liu et al., 2022) revealed that these divergent floral traits correspond to distinct pollinators: the trumpet-flowered species are pollinated by hawk moths at night, whereas butterflies visit the recurved flowers during the day.

After being categorized by floral traits into two groups, these groups were further separated into six taxa, which were accepted by Flora of China (Liang and Tamura, 2000). Although there were pervasive diagnostic criteria to identify each of taxon based on morphology, these taxa still show some overlaps in traits (Liang and Tamura, 2000). Furthermore, their distributions are largely sympatric or partially overlapping, suggesting that gene flow among these taxa is unavoidable. As mentioned earlier, such genetic exchange would impede species formation and blur species boundaries (Mallet, 2008). However, this has not been observed in the field. Most of these taxa have maintained species boundaries, at least at the morphological level. Additional evidence (e.g., genetic and ecological data) may help clarify whether and how these plants maintain species boundaries.

In the present study, our overall aim is to better understand speciation and the maintenance of species boundaries in plants. For this purpose, we used evidence from genetic analyses, flowering time, and ecological niche modeling to delineate species boundaries and determine the divergence history within Leucolirion, a monophyletic yet dimorphic lily clade. Specifically, we clarified the species status of each taxon described in Flora of China (Liang and Tamura, 2000) and determined a complete set of genetic relationships and differentiation patterns based on population sampling of next-generation sequencing (NGS) data. We then examined how species boundaries are maintained between closely related sympatric species.

2. Materials and methods 2.1. Sampling and sequencing strategy

The extremely large genomes of Lilium species make de novo assembly of a reference genome and population resequencing a challenge. To address this, we generated a full-length transcriptome as a reference for population genomic analysis. Leaf and flower tissues were sampled from a batch of L. regale that had been transplanted to our garden. To maximize transcriptome coverage over the lifespan, tissue from three stages of blooming were collected. Samples were sent to Beijing Novogene Bioinformatics Technology Co., Ltd., where RNA extraction, library preparation, and sequencing using the PacBio Sequel platform were performed. Library preparation, sequencing, data processing, and gene function annotation followed the same procedures as described in Li et al. (2020).

For population sampling of the six Lilium species, fresh leaves and flowers were collected from natural populations and transplanted materials in our lab collected during 2018–2020. Two individuals, L. duchartrei Franch. (voucher: GYD1424) and Cardiocrinum giganteum (Wall.) Makino (voucher: GYD1400), were included as outgroup taxa. Materials for transcriptome sequencing were collected at the blooming stage, and both leaves and flowers were sampled and stored in liquid nitrogen in separated tubes. Samples were sent to Novogene for RNA extraction, library preparation and sequencing. Plastomes of each individual were obtained by genome skimming on the Illumina platform. A total of 142 individuals from 21 populations were sampled (Table S1 and Fig. 1).

2.2. Collection and analysis of flowering time data

Flowering time data for the six Lilium species were retrieved from the Chinese Virtual Herbarium (CVH, https://www.cvh.ac.cn/) and Plant Photo Bank of China (PPBC, https://ppbc.iplant.cn/). For each accession, the collection date the specimen or the shooting date of the photograph was imported into IBM SPSS Statistics for Windows, v.27.0 (IBM Corp., 2020). Data were then grouped by different species and by species sets that are close in morphology and flowering time. Normality tests were conducted on both the overall and the grouped data using the Kolmogorov–Smirnov and Shapiro–Wilk tests. Results indicated that all tested datasets were non-normally distributed at a significance level of (P < 0.05) (Table S2). Subsequently, the variability of flowering time among subgroups in the monthly data was further analyzed using the Kruskal–Wallis one-way ANOVA (k samples) method in SPSS. Asymptotic significances (for two-sided tests) were utilized, with a significance level of 0.05. The significance values were adjusted using Bonferroni correction to account for multiple tests.

2.3. Ecological niche differentiation

To compare the niches of the six Lilium species, we quantified the niche overlap and width from two different perspectives: environmental space (E-space) and geographic space (G-space). Niche comparisons for G-space are based on the geographic range of species, which in many studies is often inferred from ecological niche modeling (ENM) (Lu et al., 2023), while niche comparisons for E-space typically quantify the niche of species using environmental information obtained directly from occurrence sites (Brown and Carnaval, 2019). Environmental variables used in this study include topographical factors and 19 contemporary bioclimatic factors with a resolution of 30 s (http://www.worldclim.org/). After removing highly correlated environmental variables (|r| > 0.75), eight environmental factors were retained for further analyses: Isothermality (Bio3), Max Temperature of Warmest Month (Bio5), Min Temperature of Coldest Month (Bio6), Precipitation of Wettest Month (Bio13), Precipitation of Driest Month (Bio14), Precipitation Seasonality (Bio16), Slope, and Aspect Index.

To estimate the niche overlap between six Lilium species pairs in E-space, we employed Schoener's D metric (Schoener, 1968). This metric quantifies the degree of niche similarity between different species, with higher D values indicating greater niche overlap. For each species, we generated a minimum convex polygon (MCP) and expanded the MCP's buffer zone based on an average annual diffusion distance of 10 km to define the geographical background. We then randomly chose 1000 points for each species in the background areas and extracted environmental data based on presence records and background points. Subsequently, the values of Schoener's D and niche similarity for the six Lilium species pairs were compared using the PCA-env method provided by the 'ecospat' v.4.10 R package (Di Cola et al., 2017). In the comparison of ecological niche in G-space, we first employed the Maxent v.3.4.4 (Phillips and Dudík, 2008) to predict potential distribution for each species, which was then transformed to binary suitability maps by maximizing the sum of sensitivity and specificity (maxSSS) using 'dismo' R package v.1.3.9 (Hijmans et al., 2022). Finally, we used ENMTools v.1.4.4 (Warren et al., 2010) to assess the Schoener's D value based on the binary suitability maps of Lilium species.

2.4. Species-level phylogeny reconstruction and divergence time estimation

To reconstruct a species-level phylogeny, we merged all individuals from a species and assembled the transcriptomes for each species, treating different samplings within the same species as biological replicates. The Illumina NGS data was first filtered using fastp v.0.21.0 (Chen, 2023) with default parameters, and then assembled using Trinity v.2.9.1 (Haas et al., 2013) with default parameters. Assemblies were filtered using Trinity's built-in script (i.e., 'align_and_estimate_abundance.pl' and 'filter_low_expr_transcripts.pl') for initial quality control. Redundant sequences were removed using CD-HIT v.4.8.1 (Li and Godzik, 2006) with the following parameters: -c 0.95, -n 10, -M 30000, and -T 10. The integrity of the assemblies was evaluated using Benchmarking Universal Single-Copy Orthologs v.5.7 (BUSCO) (Manni et al., 2021), with the ortholog dataset set to Liliopsida (n = 3236). TransDecoder (https://github.com/TransDecoder) was used to translate transcriptome assemblies into protein sequences. OrthoFinder v.2.5.5 (Emms and Kelly, 2019) was used to cluster orthologous genes, with single-copy genes filtered and used for subsequent species tree construction. All coding sequences (CDS) were concatenated by species using custom python scripts. Four-fold sites were obtained with perl script 'obtain_4d.pl' (https://github.com/YanCCscu/comp_genomic).

A Maximum Likelihood (ML) tree was built using IQ-TREE v.2.0 (Minh et al., 2020) with ultrafast bootstrap set to 1000, and SH-like approximate likelihood ratio test set to 1000 replicates. The divergence time of the six species was estimated using BEAST v.2.17 (Bouckaert et al., 2019). BEAST analyses were performed under a strict clock approach with a Yule tree prior, random starting trees, and appropriate nucleotide substitution model (GTR+G+I). Calibration of the crown age of the six species was conducted using a normal distribution with a mean of 2.00 (Gao et al., 2015), standard deviation (s.d.) of 0.5, and 95% quantile ranging from 1.18 to 2.83 million years ago (Mya). Similarly, the crown age of the trumpet species was calibrated using a normal distribution with a mean of 1.90, s.d. of 0.4, and 95% quantile ranging from 1.24 to 2.56 Mya. Markov Chain Monte Carlo (MCMC) chains were run for one million generations, with a sampling frequency of 1000 per generation.

2.5. Plastome phylogeny

To obtain a chloroplast phylogeny, we conducted genome skimming on all individuals. DNA extraction, library preparation, and Illumina sequencing were detailed in our previous study (Yuan and Gao, 2024). Approximately 13 Gb of raw data were filtered and evaluated using fastp with default parameters. Complete plastomes were assembled by GetOrganelle v.1.7.6.1 (Jin et al., 2020). Orientation selection was refined through multiple sequence alignment with MAFFT v.7 (Katoh and Standley, 2013) and was manually verified in Geneious Prime v.2023.1.2 (Biomatters Ltd., Auckland, New Zealand). Haplotypes from each population were selected as representatives for tree construction using DnaSP v.6.12.03 (Rozas et al., 2017).

Sequence alignments were generated using MAFFT v.7. Filtering of ambiguously aligned sites was performed using Gblocks v.0.9b (Talavera and Castresana, 2007) with default parameters. Nucleotide substitution models were estimated using ModelTest-NG v.0.1.7 (Darriba et al., 2020). Phylogenetic trees were inferred using Maximum Likelihood (ML) and Bayesian Inference (BI). ML trees were constructed using RAxML-NG v.1.2.0 (Kozlov et al., 2019) with GTR+I+G4 model and 1000 bootstrap replicates. BI analyses were performed using MrBayes v.3.2 (Ronquist et al., 2012) with the nucleotide substitution model GTR+G+I (lset nst = 6; rates = invgamma). For each analysis, the posterior probability was estimated with two independent MCMC chains (10 million generations) with the preliminary 25% of sampled data discarded as burn-in.

2.6. Transcriptome genotyping, population genetic structure and nuclear phylogeny

Clean transcriptome sequencing data for each individual were mapped to the de novo assembled full-length transcriptome of Lilium regale utilizing BWA-MEM v.0.7.17 (Li, 2013). The alignment process incorporated soft clipping for supplementary alignments. Samtools v.1.2 (Li et al., 2009) was used to filter low-quality reads based on the following criteria: (1) mapping quality score less than 20, (2) multiple best hits, (3) unmapped reads (FLAG tag of 4), and (4) adjusted mapping quality below 50. The Genome Analysis Toolkit (GATK) v.4.0.11.0 (Auwera and O'Connor, 2020) was utilized for read deduplication with MarkDuplicates, while variant calling was performed using HaplotypeCaller. The parameters for HaplotypeCaller included a minimum base quality score of 20, output-mode set to EMIT_ALL_ACTIVE_SITES, and ERC set to GVCF. The resulting GVCF files from all samples were merged via GenomicsDBImport and subjected to genotyping, filtering, and refinement processes using GenotypeGVCFs, SelectVariants, and VariantFiltration. Filtering parameters applied were 'QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0'.

Further SNP filtering was performed with bcftools (Li, 2011) to retain only bi-allelic sites and sites with a minor allele count (MAC) ≥ 2. PLINK v.1.90b6.10 (Purcell et al., 2007) was used to filter out sites with high linkage (r² > 0.1). Population genetic structure was assessed using Bayesian clustering in ADMIXTURE v.1.23 (Alexander et al., 2009) based on unlinked SNPs. The number of genetic clusters (K) was set to vary from 2 to 8, with the optimal K determined by the lowest cross-validation error. To investigate population genetic relationships further, the SNP dataset was subjected to a principal component analysis (PCA) in GCTA v.1.24 (Yang et al., 2011) and a maximum likelihood (ML) tree analysis in IQ-TREE v.2.0 (Minh et al., 2020). For nucleotide alignment, filtered unlinked SNPs were processed using vcf2phylip v.2.0 (Ortiz, 2019). The alignment was trimmed with trimAl v.1.4 (Capella-Gutiérrez et al., 2009) to remove ambiguous sites. Phylogenetic relationships of all individuals were then constructed using IQ-TREE based on the unlinked SNPs.

2.7. Demographic history inference

To estimate the changes in population size over time for each species, we estimated folded site frequency spectrum (SFS) using easySFS.py (https://github.com/isaacovercast/easySFS), incorporating all samples from the six species. The SFS was then input into Stairway Plot 2 (Liu and Fu, 2020) to reconstruct demographic history for each species. A mutation rate (μ) of 2.5e–8 per site per generation and a generation time of one year were used for the calculations. For simulating a comprehensive demographic scenario that includes divergence, population size changes, and gene flow among the six species, we utilized coalescent simulations with the composite-likelihood method implemented in fastsimcoal v.2.6 (Excoffier et al., 2013). Joint folded 2D-SFS were constructed using unlinked SNPs with easySFS.py. To minimize the impact of varying levels of missing data among groups, SNP data for each group were down-projected to an SFS with equal sampling size (26) across groups.

Given the different topologies in the phylogenetic trees reconstructed by SCNGs and unlinked SNPs, we initially tested four models to determine the most likely divergence scenario among the six species (Fig. S1). The existence and magnitude of gene flow between species, as well as changes in population size, were incorporated into the models, with additional manual fine-tuning of the parameters by comparing the differences between simulated and observed SFS (see supplementary dataset). Prior distributions of model parameters were set based on results from BEAST and Stairway Plot 2 analyses. For each scenario, we conducted 100, 000 coalescent simulations per likelihood estimation (-n 100, 000) and performed 40 expectation-conditional maximization (ECM) cycles (-L 40) as command line parameters for each run. The Akaike information criterion (AIC) was used to compare different models, where AIC = 2k – 2ln(MaxEstLhood), with k representing the number of parameters estimated by each model, and MaxEstLhood being the maximum likelihood function value for each model. The best scenario was run 100 times to obtain optimal parameter estimation. Subsequently, 100 independent DNA polymorphism datasets were simulated as joint SFSs conditional on the estimated demographic parameters. Maximum likelihood (ML) analysis was then applied to each joint SFS over 40 ECM cycles to obtain confidence intervals (CIs) for the final estimates.

2.8. Genomic diversity and detection of selection

Nucleotide diversity (π) for each species and population genetic differentiation (F_ST) between each pair of species were calculated using pixy v.1.2.7 (Korunes and Samuk, 2021), focusing on the first 10 kb of each transcript. To detect signatures of positive selection during the speciation of flowers with two distinct types (recurved and trumpet), we grouped the species as follows: Lilium henryi and L. rosthornii (recurved flowers) and L. leucanthum, L. sargentiae, L. sulphureum (trumpet flowers). Given that L. regale is restricted to a small region and is independent both geographically and genetically (see results), this species was excluded from the selection analysis to minimize the confounding effects of genetic drift.

For detecting regions under selection, we identified transcripts showing significantly elevated genetic differentiation (top 5% F_ST) between the flower groups and significant reductions in genetic diversity within either group (top and tail 5% π ratio). These transcripts were considered candidates for divergent selection. Further, functional enrichment analysis, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses, was performed on the orthologous genes in Arabidopsis thaliana using KOBAS v.3.0 (Bu et al., 2021).

3. Results 3.1. Species level phylogeny reconstructed by transcriptome-based orthologous genes

The representative transcriptome assembled for each species ranged from 370 to 516 Mb in length, with scaffold N50 ranging from 823 to 1144 bp (Table S3). The completeness of the assemblies, as assessed by BUSCO, ranged from 89.9% (2909) to 94.9% (3071) (Fig. S2). Through Orthofinder2, a total of 7050 single-copy orthologous genes (SCNGs) shared by the six species were predicted. A maximum likelihood (ML) phylogenetic tree constructed based on these SCNGs revealed two major clades: one clade consisted of the two species with recurved flowers (Lilium henryi and L. rosthornii), while the other included the remaining four species with trumpet flowers (Fig. 2A and B). Within the latter clade, L. regale was identified as the basal taxon, followed by the divergence between L. leucanthum and the clade containing L. sargentiae and L. sulphureum (Fig. 2A). BEAST analysis estimated the crown age of the six species to be approximately 2.36 Mya. The divergence between L. henryi and L. rosthornii was estimated to have occurred around 0.63 Mya, while the diversification of the clade containing the four trumpet-flowered species was dated between 1.08 Mya and 0.67 Mya (Fig. 2A).

3.2. Flowering time and ecological differentiation of six Lilium species

To assess variation in flowering time among six Lilium species, we performed non-parametric tests using monthly data from independent samples. The results revealed significant differences in flowering time between the two recurved species, L. henryi and L. rosthornii (P < 0.05), although these differences were not significant after applying the Bonferroni correction (P > 0.05; Fig. 2C and Table S4). Significant differences in flowering time were also found among the four species with trumpet flowers, with the exception of the pair L. sargentiae and L. leucanthum (Fig. 1C and Table S4). No significant differences were observed in flowering time between the recurved L. henryi and the trumpet-flowered L. sulphureum (Fig. 2C and Table S4). Moreover, the flowering time of L. regale did not resemble that of any other species or species sets (Table S4).

Ecological niche analysis indicated substantial niche overlap between the recurved-flower species Lilium henryi and L. rosthornii in both ecological (E-space) and geographic (G-space) contexts, with Schoener's D values of ≥ 0.446 (P < 0.05; Table 1). Similar patterns were found among the three trumpet-flower species (L. sulphureum, L. sargentiae, and L. leucanthum), where Schoener's D values were ≥ 0.454, most of which were statistically significant. In contrast, L. regale exhibited few niche overlaps with the other trumpet-flowered species (Schoener's D ≤ 0.022) and no overlap (Schoener's D = 0) with the recurved-flower species (Table 1). Members of recurved-flower species also showed relatively minimal overlap with those of the trumpet-flower species, particularly in E-space, with Schoener's D values of ≤ 0.162 (Table 1).

3.3. Full-length transcriptome assembly, population structure and phylogenetic relationships

A full-length transcriptome reference was obtained using the PacBio sequencing, which generated a total of 12, 853, 869 subreads, encompassing 21.36 Gb of data. After removing redundancy, the assembled transcriptome reached a length of 33, 483, 797 bp (N₅₀ = 2743), containing 31, 913 transcripts and 15, 670 genes. Functional annotation of the transcriptome showed that 15, 002 genes were annotated by at least one database (Table S5).

All population samples and outgroups were mapped against this reference transcriptome. Subsequent haplotype calling and filtering procedures resulted in 3, 302, 381 SNPs, with 624, 637 unlinked SNPs retained for further population structure and phylogenetic analyses. Admixture analysis demonstrated a clear population structure from K = 2 to K = 8, with the lowest cross-validation error observed at K = 3 (Fig. S3). When K = 2, individuals of the two flower types formed two distinct clusters, and an additional cluster emerged for Lilium regale at K = 3 (Fig. 3A). At K = 6, the two species with recurved flowers were recognized by different genetic clusters, whereas the three trumpet-flower species (L. sulphureum, L. sargentiae, and L. leucanthum) exhibited signs of admixture and unstable ancestries (Fig. 3A). Principal Component Analysis (PCA) corroborated these findings, revealing three distinct genetic clusters and a close relationship among the three trumpet-flower species (Fig. 3B).

A phylogenetic tree constructed from unlinked SNPs showed a topology (Fig. 3C) similar to the species phylogeny (Fig. 2A). However, the relationships among individuals from the three species with trumpet flowers differed from the species-level phylogeny, with Lilium leucanthum nested within L. sulphureum, and both species forming a clade sister to L. sargentiae (Fig. 3C). The plastome haplotype tree presented a different topology, revealing conflicts between the two recurved-flower species and among the three trumpet-flower species (Fig. 3D). This suggests the possibility of hybridization and/or introgression within these groups.

3.4. Demographic history and gene flow

Stairway Plot 2 analysis revealed relatively stable population sizes for Lilium leucanthum, L. regale and L. henryi. In contrast, the other three species experienced varying degrees of population contraction, notably significant in L. sulphureum and L. sargentiae (Fig. 4A). A composite likelihood test of demographic models (Table S6 and Fig. S1) indicated that L. leucanthum diverged from the most recent common ancestor (MRCA) of L. sulphureum and L. sargentiae. This finding aligns with the species phylogeny constructed using single-copy orthologous genes (SCNGs).

Throughout the diversification of the six species, multiple instances of gene flow were detected, exhibiting various magnitudes and directions (Fig. 4B). Notably, continuous symmetric gene flow (M3 = 6.21e–7) was detected among Lilium leucanthum, L. sulphureum, and L. sargentiae since their divergence. Additionally, a significant pulse of symmetric gene flow occurred between L. henryi and L. rosthornii at approximately 15–20 thousand years ago (Kya) (M1 = 3.87e–5). Significant asymmetric gene flow was observed from the MRCA of species with trumpet flowers to the MRCA of the recurved flower clade at the onset of their differentiation (2.1–1.99 Mya, M9 = 3.28e–5; 1.99–1.23 Mya, M8 = 6.06e–4). Other instances of weak gene flow included that between L. regale and other trumpet-flowered species (M2 = 1.09e–7), as well as from trumpet-flowered species to recurved-flowered species since 1.23 Mya (M7 = 1.3e–8 from 1.23 Mya to 381 Kya; M5 = 7.93e–8 from 381 Kya to 20 Kya; M4 = 5.6e–7 from 20 Kya to 15 Kya). All gene flow events ceased after 15 Kya.

3.5. Genetic diversity and signatures of positive selection

Estimates of nucleotide diversity (π) suggested that the three species with trumpet flowers (Lilium leucanthum, L. sulphureum, and L. sargentiae) have relatively high genetic diversity compared to other species, while L. regale has lower π values. L. rosthornii, with recurved flowers, has the lowest genetic diversity (Fig. 5A). We observed significant genetic differentiation between the two flower shape types, with an average F_ST of 0.55 or higher (Fig. 5B). Genetic differentiation among the three species with trumpet flowers was minimal, likely due to continuous gene flow (F_ST ≤ 0.11). However, L. regale was genetically independent from the other three trumpet-flowered species (F_ST ≥ 0.25; Fig. 5B). The two species with recurved flowers also exhibited genetic divergence from each other (F_ST = 0.23; Fig. 5B).

A genomic scan for elevated genetic divergence between groups and lowered genetic diversity within either group identified 105 transcripts under potential positive selection (Fig. 5C). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses revealed that these genes are involved in responses to environmental stimuli (e.g., circadian rhythm and MAPK signaling pathway), various metabolic and biosynthetic processes (e.g., anthocyanin metabolism), and different cell wall and structural pathways (e.g., cellulose biosynthesis; Fig. 5D).

4. Discussion

The completely different flora forms (two forms in six species) in our study exemplify parallel evolution in the genus Lilium, which has been widely documented in previous studies (Gao et al., 2015; Yuan and Gao, 2024). This phenomenon creates challenges for systematic classification due to discordances between morphology-based taxonomy and molecular phylogeny. These distinct floral morphologies (whitish-trumpet and orange-recurved) were assigned to different sections based on their strong resemblance to these funnel and 'turkey-cap' members assembled as sect. Leucolirion and sect. Sinomatagon, respectively (Comber, 1949). However, molecular studies have confirmed their monophyletic status within Lilium (Gao et al., 2013; Du et al., 2014; Li et al., 2022). This monophyletic yet dimorphic lily group provides an excellent case study for investigating speciation and species boundary maintenance. Our comprehensive analysis of the six Lilium species using transcriptome-based orthologous genes, population genomics, flowering time data, and ecological data has provided valuable insights into their evolutionary history, speciation processes, and adaptive strategies. Our results indicate that diversification of these species has been shaped by a complex interplay between phylogenetic relationships, flower morphology, ecological differentiation, and gene flow.

4.1. History of species divergence within clade Leucolirion

Our transcriptome-based molecular dating estimates the crown age of the six species at early Pleistocene (ca. 2.36 Mya; Fig. 2A), when the recurved and trumpet flower groups diverged from each other. Subsequent diversification events within groups ranged from the early to middle Pleistocene (between 1.08 Mya and 0.63 Mya), a period marked by pronounced climate fluctuations and glacial cycles that promoted rapid angiosperm diversification (Liu et al., 2023). Global climate oscillation during Pleistocene drove numerous speciation and lineage divergence events, causing rapid radiation in many genera, as reported in Rhododendron (Xia et al., 2022), Rodgersia (Ma et al., 2017) and Primula (Ren et al., 2015). However, range expansion-contraction during the glacial–interglacial cycle and/or incomplete range isolation often cause genetic exchange during divergence (Chen et al., 2019). Indeed, our demographic model, inferred by fastsimcoal2 (Fig. 4B), reveals significant intragroup gene flow, indicating historical hybridization or introgression. This genetic exchange likely contributed to reticulated evolution and, in some cases, blurred species boundaries. Gene flow can result in shared maternal genotypes among closely related species, potentially explaining the discrepancies we observed between nuclear and plastid phylogenies, particularly among the trumpet-flowered species (Fig. 3D; Du et al., 2009; Mo et al., 2022). Specifically, significant symmetric gene flow occurred between Lilium henryi and L. rosthornii approximately 15–20 Kya (Fig. 4B), likely causing their plastome haplotypes to nest with each other (Fig. 3D). Additionally, gene flow from trumpet members to recurved members resulted in chloroplast capture, leading some plastome haplotypes of recurved members to cluster with those of trumpet members (e.g., LRO_GYD309-10; Fig. 3D). Similarly, gene flow among the three trumpet-flowered species has shaped complex evolutionary relationships among their plastome haplotypes: some individuals of L. leucanthum were located at the base of the clade, reflecting either earlier divergence or potential gene flow from L. regale (Fig. 4B), while others clustered with L. sulphureum and L. sargentiae, resulting from more frequent gene flow among them.

4.2. Clear boundary between the recurved and trumpet flower groups

The separation of this clade into recurved and trumpet flower groups is strongly supported by multiple molecular datasets. The transcriptome-based phylogeny, constructed from 7050 single-copy orthologous genes, clearly delineates two major clades corresponding to these flower morphologies (Fig. 2). This division is further corroborated by the nuclear SNP data, which shows distinct genetic clusters for recurved and trumpet-flowered species in both admixture and principal component analyses (Fig. 3), with significant genetic differentiation between groups (F_ST ≥ 0.55). Pollinator-mediated reproductive isolation plays a crucial role in maintaining species boundaries (Grant, 1994; Brothers et al., 2014; Van der Niet et al., 2014). Our previous research on pollinator differentiation (Liu et al., 2022) demonstrated distinct pollination syndromes between recurved and trumpet-flowered Lilium species, with butterflies visiting the former during the day and hawk moths pollinating the latter at night. Therefore, the adaptation to different pollinators seems to play a quite important role in isolation of the two major lineages, despite some overlaps in flowering time and geographic distribution (Fig. 2 and Table 1).

Ecological niche analysis further supports this differentiation, showing substantial niche overlap within flower types but limited overlap between them (Fig. 2D and Table 1). The genes under potential positive selection between the two flower types, particularly those involved in responses to environmental stimuli, metabolic processes, and structural pathways, provide a genetic basis for the observed ecological differentiation (Fig. 5C and D) despite the presence of gene flow (Fig. 4B). The most obvious difference between the two groups is flower color (whitish-trumpet and orange-recurved; Fig. 2B), which is reflected in the positive selection of genes related to anthocyanidin metabolic process (Fig. 5D). Notably, these species have evolved different adaptations to attract distinct pollinators—the trumpet-flowered species are pollinated by hawk moths at night, while butterflies visit the recurved flowers during the day—with multiple aspects of genetic evidence supporting this differentiation. For example, the significant difference in nectar volume between groups reported in Liu et al. (2022) is evidenced in the positive selection of genes related to pathways and biological processes included carbon metabolism, carbohydrate, and sugar mediated signaling pathways. Additionally, genes associated with circadian rhythm (Fig. 5D) could contribute to the temporal display of attractive qualities, such as scent emission in a daily rhythmic fashion to attract pollinators (e.g., Fenske et al., 2018). These genetic adaptations highlight how selective pressures on specific pathways can drive ecological differentiation, especially pollinator-related divergence, allowing closely related species to occupy distinct temporal and environmental niches while maintaining reproductive isolation. Such adaptive changes likely facilitate the occupation of diverse niches and contribute to the maintenance of species boundaries despite potential gene flow (Nosil, 2008).

4.3. Asynchronous differentiation within groups

While substantial differentiation was found between the two groups of flower types, more complex patterns emerged within each group. Moderate genetic differentiation was observed between Lilium regale and other trumpet-flower species, as well as between the two recurved-flower species (Fig. 5B). L. regale showed nearly no gene flow with its trumpet-flowered relatives, and gene flow between the two recurved-flower species existed only for a short period (Fig. 4B). This moderate differentiation aligns with our findings on flowering time and ecological niche analysis (Fig. 2C and Table 1). For instance, L. regale blooms nearly a month earlier than the other trumpet members and exhibits little to no ecological and geographic overlap with other species. Despite sharing both distribution ranges and pollinators (hawk moths) with L. leucanthum, L. regale maintains substantial genetic separation from the latter (Fig. 3). In fact, the floral tubes of L. regale were significantly shorter than other trumpet-flowered species, attracting more hawk moths with short proboscises, such as Deilephila elpenor and Ampelophaga rubiginosa (Liu et al., 2022), likely contributing to its more independent evolution. Similarly, the recurved members L. henryi and L. rosthornii demonstrate a tendency for flowering time and niche differentiation. Specifically, L. henryi prefers sunny, drought-prone slopes, whereas L. rosthornii inhabits humid valleys along streams. These lines of evidence indicate a potential contribution to reproductive isolation within groups through temporal isolation (Coyne and Orr, 2004; Gaudinier and Blackman, 2020) and niche divergence (Schluter, 2000; Schluter and McPeek, 2000; Wogan and Richmond, 2015).

In contrast, the other three trumpet-flowered species (Lilium leucanthum, L. sulphureum, and L. sargentiae) show limited genetic differentiation among themselves (Fig. 5B), consistent with the continuous symmetric gene flow (M3 = 6.21e–7) detected since their divergence (Fig. 4B). This limited differentiation is further reflected in their substantial overlap in geographic distribution and ecological niches (Table 1). Additionally, unlike L. regale, with its relatively smaller flowers, these three species exhibit large trumpet-shaped flowers that are not only adapted to long-proboscid pollinators but also attract various hawkmoth species, as those with shorter proboscises can alight and crawl into the flowers in search of nectar, benefiting from an expanded pollination niche breadth (Liu et al., 2022). This enables potential cross-specific visits that may lead to gene flow between species and blur the species boundaries. However, trends toward divergence are evident. For example, the divergence between L. leucanthum and the L. sargentiae + L. sulphureum clade is supported by both molecular (Fig. 2, Fig. 4B) and morphological evidence, with both L. sargentiae and L. sulphureum possessing axil bulblets while L. leucanthum lacks them (Liang and Tamura, 2000). Besides, L. sulphureum blooms significantly later than L. leucanthum and L. sargentiae (Fig. 2C). The demographic history inferences also reveal diverse population size trajectories: while L. sargentiae has experienced recent population contractions, L. leucanthum and L. sulphureum show population stability (Fig. 5B). These varied histories likely reflect species-specific responses to historical climate and geographical changes, illustrating ecological resilience and the influence of historical biogeography on contemporary population dynamics (Avise, 2000).

5. Conclusions

In summary, the evolution of the Leucolirion clade is characterized by an initial split between recurved and trumpet-flowered lineages, driven by morphological adaptations to different pollination strategies. Subsequent diversification within each group involved flowering time shifts, ecological niche differentiation, and localized adaptations. The recurved-flower species (Lilium henryi and L. rosthornii) diverged primarily through temporal isolation, while trumpet-flowered species exhibit more complex diversification patterns. L. regale stands out with its distinct niche and genetic independence, whereas ongoing gene flow among L. leucanthum, L. sulphureum, and L. sargentiae suggests closer relationships and potentially blurred species boundaries. These findings collectively support the idea that species divergence and maintenance in Lilium are driven by a combination of adaptive and non-adaptive processes, including pollinator-mediated selection, ecological niche differentiation, historical demographic events, and natural selection. This research not only advances our understanding of Lilium evolution but also offers insights into the processes that underlie plant speciation and adaptation. Future studies incorporating broader species sampling and more detailed ecological data will further elucidate the evolutionary mechanisms in this diverse plant group. Notably, the recently published L. davidii genome (Xu et al., 2024) could enable more accurate phylogenetic reconstructions and provide deeper insights into the evolutionary history of Lilium. This genomic resource could also reveal structural variations that help elucidate the mechanisms driving adaptation and speciation within the genus, particularly in the context of ecological differentiation and reproductive isolation.

Acknowledgements

We gratefully acknowledge Prof. Xin-Fen Gao and Mr. Jun-Yi Zhang of Chengdu Institute of Biology, Chinese Academy of Sciences, Dr. Pan Li of Zhejiang University, and Mr. Melvyn Herbert for their support with sample collection. This work was supported by the National Natural Science Foundation of China (NSFC Grant Nos. 32171605, 31770406) and Natural Science Foundation of Sichuan Province (Grant No. 2023NSFSC0141) to Y.D. Gao.

CRediT authorship contribution statement

Yu Feng: Writing – review & editing, Writing – original draft, Visualization, Methodology, Formal analysis, Conceptualization. Chaochao Yan: Writing – original draft, Visualization, Methodology, Formal analysis, Data curation. Wen-Qin Tu: Writing – original draft, Methodology, Formal analysis, Data curation. Yu-Mei Yuan: Writing – original draft, Visualization, Formal analysis, Data curation. Jing-Bo Wang: Writing – original draft, Software, Formal analysis. Xiao-Juan Chen: Writing – original draft, Visualization, Formal analysis. Chang-Qiu Liu: Validation, Methodology, Investigation. Yundong Gao: Writing – review & editing, Writing – original draft, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization.

Data availability

The sequencing data generated in this study has been deposited to NCBI SRA database under the BioProject accession PRJNA1141931. Supplementary datasets including configuration files and evaluation of demographic scenarios tested by fastsimcoal2 were deposit to Mendeley Data (10.17632/c7sgdwvrzh.1). Other data will be made available on request.

Declaration of competing interest

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.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2024.12.005.