1. Introduction

Lycium L. (Solanaceae), comprising approximately 80 species, exhibits a remarkable intercontinental disjunct distribution, with most species occurring in the arid or semi-arid regions of South America (ca. 30 spp.), North America (ca. 20 spp.), South Africa (ca. 20 spp.), and Eurasia (from Europe to China and Japan: ca. 10 spp.), while a few others are native to Australia (1), and several islands in the Pacific Ocean (2 spp.; see also Symon, 1991; Fukuda et al., 2001; Levin and Miller, 2005). These plants are typically deciduous spiny shrubs or small trees, including various medicinal and edible species (Yao et al., 2011; Wetters et al., 2018; Tian et al., 2021; Li et al., 2024; Xiong et al., 2025), and some (e.g., Lycium ruthenicum Murray), which are also used for the restoration of degraded or saline soils, particularly in the arid regions of China (Ma et al., 2025). This genus produces small, fleshy berries, typically red or orange, and, more rarely, purple or black, which are considered adaptations for bird dispersal (Fukuda et al., 2001; Yao et al., 2018). Two competing hypotheses have been proposed to explain the intercontinental disjunct distribution of Lycium: (1) ancient vicariance linked to Gondwanan fragmentation (Symon, 1991); and (2) long-distance dispersal (LDD) events (Raven and Axelrod, 1974). Early molecular studies employing nuclear internal transcribed spacer (ITS), chloroplast intergenic spacer matK/trnL–trnF, and granule-bound starch synthase (GBSSI) gene sequences supported the LDD hypothesis. Molecular dating calibrated a crown age of ~29.4 ± 9.7 million years ago (Ma) for Lycium. This supported an American origin (though disputed between North and South America), followed by a single dispersal to Africa and subsequent colonization of Eurasia, with bird dispersal of the fleshy fruits likely facilitating these events (Fukuda et al., 2001; Levin and Miller, 2005; Miller et al., 2011). Cao et al. (2021) suggested a reverse Asian-to-North American dispersal based on single-nucleotide polymorphism (SNP) data. These previous studies were constrained by distinct limitations despite using high-resolution SNP data and had incomplete taxon sampling (Cao et al., 2021), whereas earlier studies (Fukuda et al., 2001; Levin and Miller, 2005; Miller et al., 2011) were limited by the low phylogenetic resolution inherent in using only a few loci. This resulted in weakly supported phylogenies and uncertain biogeographic inferences. Recently, however, the plastome phylogeny of Lycium inferred a North American origin with high confidence (Yisilam et al., 2025).

The intra-generic classification of Lycium has a long history, primarily based on morphological characteristics, such as floral structure and fruit type (Hitchcock, 1932; Chiang-Cabrera, 1981; Bernardello, 1986; Fukuda et al., 2001). Traditional systems, focusing on morphologically diverse American species, recognize several sections (e.g., Eulycium/Lycium, Mesocope, Sclerocarpellum, Schistocalyx) (Hitchcock, 1932; Chiang-Cabrera, 1981; Bernardello, 1986). However, molecular phylogenetic studies have consistently challenged the monophyly of these morphological sections, indicating that key diagnostic traits are either homoplastic or that their evolution is complicated by processes such as incomplete lineage sorting (ILS) and introgression (Fukuda et al., 2001; Levin and Miller, 2005). The recent plastome phylogeny reported by Yisilam et al. (2025) further corroborates these incongruences in several sections. However, because this study only used maternally inherited plastid data, the primary cause of this phylogenetic (taxonomic-molecular) discordance could not be fully clarified.

Interlineage introgression and ILS are the primary contributors to phylogenetic discordance (Delsuc et al., 2005; Schumer et al., 2014; Guo et al., 2023). The genus Lycium presents a compelling case study for such complexity, as it is prone to reticulate evolution. Natural hybridization has been documented (e.g., L. afrum × L. ferocissimum; Minne et al., 1994), and several polyploid species in southern Africa (e.g., L. arenicolum and L. villosum) have been hypothesized to have hybrid origins (Minne et al., 1994; Levin and Miller, 2005). Polyploidy is widespread in this genus (base number x = 12, with known diploid, tetraploid, and hexaploid species), and transitions to gender dimorphism have occurred multiple times, often associated with polyploidization events (Minne et al., 1994; Miller and Venable, 2000; Levin et al., 2007; Gong et al., 2022). These factors can cause significant cytonuclear discordance and conflict among nuclear gene trees; however, their relative contributions within Lycium remain unresolved.

Reticulate processes, such as introgression (Delsuc et al., 2005; Schumer et al., 2014, 2015; Hodel et al., 2021, 2022; Liu et al., 2025), disrupt the shared genealogical history of nuclear and plastid genomes (Dong et al., 2021; Thureborn et al., 2024; Chen et al., 2025). Moreover, ILS may further amplify the topological incongruence between gene trees and species trees (e.g., Cai et al., 2021). These processes create pervasive phylogenetic discordance driven by intrinsic differences in inheritance, effective population size, and evolutionary rate between the two genomic compartments (Hu et al., 2023).

High-throughput sequencing, combined with rigorous ortholog identification, provides the capacity to resolve complex phylogenetic discordance (Brassac and Blattner, 2015; Liu et al., 2022a; Gu et al., 2024; Lin et al., 2025). The integration of biparentally inherited nuclear and uniparentally inherited plastid and mitochondrial genomes data offers an independent genealogical perspective (Qiu et al., 1999; Liu et al., 2022b; Xue et al., 2024). This allows researchers to distinguish the effects of introgression from those of ILS using coalescent models and phylogenetic networks (Degnan and Rosenberg, 2006, 2009; Pease and Hahn, 2015; Solís-Lemus et al., 2017; Hodel et al., 2022), an approach that has clarified similar complexities in other taxa (Ma et al., 2021; Pezzi et al., 2024).

To address these limitations, we established a tripartite phylogenomic framework based on large-scale nuclear gene data (378 single-copy nuclear genes), mitochondrial genomes, and published plastid data (Yisilam et al., 2025). This study aims to achieve several key objectives: (1) to establish a comprehensive nuclear phylogenetic framework; (2) for the first time, to systematically assess the relative impacts of incomplete lineage sorting (ILS) and hybridization/introgression on phylogenetic conflicts, including tripartite nuclear–plastid–mitochondrial discordance and intra-nuclear gene-tree discordance, throughout the entire genus; (3) to identify and model specific reticulate events utilizing a range of complementary gene flow detection and network inference methodologies; and (4) to deduce a biogeographic history informed by nuclear dating techniques. By establishing this multi-genomic phylogenomic framework, we aimed to elucidate the complex evolutionary dynamics of Lycium and provide key mechanistic insights into the formation of its disjunct intercontinental distribution patterns.

2. Methods 2.1. Taxon sampling

In this study, we sampled 43 Lycium species from North America (19 spp.), South America (12), the Pacific islands (1), South Africa (4), Saharan Arabia (1), and Eurasia (6) (see Table S1). Our sampling also included all four traditional sections based on the American morphology (Eulycium/Lycium, Mesocope, Sclerocarpellum, Schistocalyx; see Table S3).

Leaf tissue samples from 36 Lycium species were obtained from herbarium specimens (ARIZ, University of Arizona Herbarium, and WIS, Wisconsin State Herbarium), and fresh leaves of L. qingshuiheense, L. cylindricum, and L. dasystemum were sampled from their natural habitats in Ningxia, Xinjiang, and Qinghai, China (Table S1). Raw genome sequencing data for four Lycium species, namely L. barbarum, L. chinense, L. ruthenicum, and L. ferocissimum, were downloaded from NCBI (Table S2). Raw sequencing data for eight outgroup species representing six genera of Solanaceae–Nolana paradoxa, Atropa belladonna, Datura inoxia, Datura stramonium, Nicotiana tabacum, Physalis grisea, Physalis peruviana, and Solanum dulcamara–were obtained from NCBI [see Table S2 for Sequence Read Archive (SRA) accession numbers].

2.2. DNA extraction and sequencing

We sequenced the genomes of 39 Lycium species using whole-genome sequencing (WGS). Total genomic DNA was extracted from the dried leaf tissue samples using a Plant Genomic DNA Kit (Tiangen Biotech, Beijing, China) following the manufacturer’s protocol. DNA integrity was assessed via 1% agarose gel electrophoresis, and concentration and purity (A260/A280 and A260/A230 ratios) were measured spectrophotometrically with concentrations of 50 ng/μL. For the 39 Lycium spp., Illumina paired-end libraries (150 bp) were constructed and sequenced on an Illumina HiSeq2500 platform (Illumina, San Diego, CA, USA) to achieve a minimum coverage of approximately 30×. All raw data were subjected to quality control and adapter trimming using Trimmomatic v.0.39 (Bolger et al., 2014). The final sequencing yield ranged between 50 and 90 Gb per species for Lycium (Table S4).

2.3. Probe design and assembly of Orthologous single-copy nuclear genes

Orthologous single-copy nuclear (SCN) genes were identified using a multi-step approach. Transcriptomes of Lycium barbarum, L. chinense, and L. ruthenicum. were obtained from GenBank. Raw read quality was controlled using Fastp v.0.23.2 (Chen et al., 2018), removing adapters, low-quality bases (Phred score < 20), and short reads (< 50 bp). Filtered reads were de novo assembled in Trinity v.2.4.0 (Grabherr et al., 2011) using default parameters, followed by redundancy reduction in Cd-hit v.4.6.5 (Fu et al., 2012) with a 95% identity threshold. TransDecoder v.5.7.1 (Haas et al., 2013) was used to predict the candidate reading frames. Putative orthogroups were inferred using OrthoFinder2 v.2.3.11 (Emms and Kelly, 2019) to generate candidate reference sequences for probe design. Second, candidate orthogroup sequences derived from the transcriptome analysis were used as reference sequences. HybPiper v.1.3.1 (Johnson et al., 2016) was used to assemble potential SCN genes from the cleaned sequencing data of four representative Lycium species, L. cylindricum, L. americanum, L. chilense, and L. californicum, based on the reference sequences. These four species were selected to ensure our probes targeted true orthologs and minimized paralogous capture, based on their broad phylogenetic representativeness across the major geographic distributions (Eurasia, North America, and South America) and evolutionary lineages within Lycium. This pipeline automates read mapping, contig assembly, and reference-guided sequence extraction. Genes flagged in the HybPiper output files, genes_with_paralog_warnings or genes_with_long_paralog_warnings, were manually excluded to establish a high-confidence set of candidate SCN gene probes (orthologous group, OG). HybPiper was used to assemble the SCN genes from the cleaned data of each sample. The filtered OG probe sequences were used in a second HybPiper run to assemble the SCN gene sequences from the resequencing data of all 39 Lycium species and eight outgroup taxa. The assembled sequences were checked using Geneious Prime v.2021.2.2 (Kearse et al., 2012) to exclude the absence of genes in any taxon. After paralog filtering, 378 orthologous SCN genes were assembled for 43 Lycium species (Table S4), which were used for all subsequent phylogenetic and biogeographical analyses.

2.4. Assembly of the Lycium ruthenicum mitogenome and retrieval of mitochondrial sequences across taxa

To conduct mitochondrial genome-based phylogenetic analysis, the complete mitochondrial DNA of Lycium ruthenicum was first assembled de novo from PacBio HiFi reads (https://dataview.ncbi.nlm.nih.gov/object/SRR31802263) using PMAT v.1.5.3, with the “-st HiFi” parameter (Bi et al., 2024). The assembly was annotated using the online website PMGA (Li et al., 2025). All gene boundaries were manually verified and adjusted using Geneious Prime v.2021.2.2 (Kearse et al., 2012), and genome maps were generated using OGDRAW v.1.3.1 (Greiner et al., 2019). Finally, the complete mitochondrial sequence of L. ruthenicum was submitted to the NCBI database under accession number PX778778.

Genomic features were characterized as follows: simple sequence repeats (SSRs) were detected using the online website MISA (https://webblast.ipk-gatersleben.de/misa/), with the minimum repeat thresholds set to 10, 5, 4, 3, 3, and 3 for mono- and hexanucleotides. Tandem and dispersed repeats were identified using Tandem Repeats Finder v.4.09 (https://tandem.bu.edu/trf/trf.html) and REPuter (https://bibiserv.cebitec.uni-bielefeld.de/reputer) (Kurtz et al., 2001), respectively. The relative synonymous codon usage (RSCU) for protein-coding genes was calculated using Mega v.11.0 (Tamura et al., 2021). Potential mitochondrial plastid DNA sequences (MTPTs) were identified by performing a BLASTn (Chen et al., 2015) alignment of the mitochondrial genome against its chloroplast genome (PX778779) counterpart using an e-value cutoff of 1e-5, identity > 80%, and alignment length > 100 bp.

For the phylogenetic analysis, 30 annotated protein-coding genes (PCGs) were extracted from the Lycium ruthenicum mitochondrial sequence. These PCGs sequences served as reference targets in HybPiper v.1.3.1 (Johnson et al., 2016) to retrieve homologous sequences from the whole-genome resequencing data of other sampled species under default parameters. The final concatenated dataset of 30 conserved PCGs shared across all taxa was used for subsequent mitochondrial phylogenetic reconstruction.

2.5. Phylogenetic analyses

Phylogenomic inferences based on 378 single-copy nuclear (SCN) genes were made using concatenation- and coalescent-based methods for Lycium. For the concatenation-based analysis, individual nuclear gene sequences were aligned utilizing MAFFT v.7.407 (Katoh and Standley, 2013) under the automatic setting. Following this, sequences were trimmed using trimAl v.1.2 (Capella-Gutiérrez et al., 2009) with parameters set to -gt 0.2, -st 0.001, and -cons 85 to effectively remove ambiguous sequences. Trimmed alignments were concatenated into a supermatrix using Geneious Prime v.2021.2.2 (Kearse et al., 2012). PartitionFinder v.2.1.1 (Lanfear et al., 2017) was used to estimate the optimal partitioning scheme for maximum likelihood (ML) inference in IQ-TREE v.2.1.4 (Minh et al., 2020). The ML tree ultrafast bootstrap branch support (BS) values were obtained via 1000 ultrafast bootstrap replicates (Hoang et al., 2018). For coalescent-based inference (ASTRAL species tree), nuclear gene trees were first independently reconstructed for each of the 378 SCN genes using IQ-TREE under the settings described above. The species tree was then estimated with ASTRAL-Ⅲ v.5.7.8 (Zhang et al., 2018), and local posterior probabilities (LPPs) were calculated to quantify branch support. Phylogenetic trees were visualized using FigTree v.1.4.4 (https://github.com/rambaut/figtree/releases/tag/v1.4.4).

Mitochondrial phylogenetic inference was conducted using a maximum likelihood framework on a concatenated alignment of mitochondrial PCGs sequences from 50 taxa (43 within Lycium and seven outgroups). Each PCGs was aligned separately using MAFFT v.7.407 (Katoh and Standley, 2013). The aligned regions were trimmed using trimAl v.1.2 (Capella-Gutiérrez et al., 2009) by removing columns with excessive gaps (> 20% of sequences) or low information content (similarity score < 0.001). The filtered alignments for all genes were concatenated into a single dataset using AMAS v.1.0 (Borowiec, 2016). The combined dataset was analyzed in IQ-TREE v.2.1.4 (Minh et al., 2020), employing the MFP + MERGE option to automatically infer the optimal partition scheme and substitution models. Node support was quantified using 1000 ultrafast bootstrap replicates. Phylogenetic trees were visualized using FigTree v.1.4.4.

2.6. Evaluate phylogenetic conflicts and quantify evolutionary dynamics 2.6.1. Quantification and visualization of topological conflicts

To visualize topological conflicts among gene trees, we performed multidimensional scaling (MDS) analysis (Duchêne et al., 2018) based on pairwise Robinson–Foulds distances between the inferred nuclear ASTRAL species tree and all individual nuclear gene trees using the R package phytools v.2.4.4 (Revell, 2024). In the resulting MDS plot, each point represents a gene tree colored according to the average nodal bootstrap support value.

To evaluate the support and consistency of internal branches in our nuclear phylogeny, we applied Quartet Sampling (QS; Pease et al., 2018). This method calculates three branch-specific metrics, namely quartet concordance (QC), quartet differential (QD), and quartet informativeness (QI), based on stochastic subsampling of taxon quartets from the alignment of 378 nuclear genes and the ASTRAL species tree.

Furthermore, topological conflicts within nuclear gene trees and ASTRAL species tree were assessed using PhyParts v.0.0.1 (Smith et al., 2015). Conflicts were calculated using a 50% bootstrap filter to exclude weakly supported relationships and were visualized using phypartspiecharts.py.

2.6.2. Assessment of incomplete lineage sorting

To assess the contribution of incomplete lineage sorting to the observed phylogenetic discordance, we performed an analysis under the ‘T1 model’ using the quartetTreeTest function from the R package MSCquartets v.2.0 (Rhodes et al., 2021), with the inferred nuclear species tree as the reference. This methodology assesses the adequacy of the multispecies coalescent model by computing quartet count concordance factors (qcCFs) for all potential taxon quartets. The findings were represented through two-dimensional probability simplex plots utilizing the quartetTestPlot function, evaluated across four distinct significance levels (α = 0.01, 0.001, 1e−04, and 1e−06).

2.6.3. Quantifying incomplete lineage sorting and historical introgression

To quantify historical introgression in Lycium, we used QuIBL v.1.0 (Quantifying Introgression via Branch Lengths; Edelman et al., 2019) to evaluate introgression signals in the context of both the nuclear gene trees and the nuclear ASTRAL species tree. For each three-taxon triplet, we computed the likelihood using competing models of ILS alone and ILS with introgression. Model selection was performed using the Bayesian Information Criterion (BIC), where BIC > 10 provided strong support for the ILS-only model (Feng et al., 2022). Triplets exhibiting an ambiguous signal (|ΔBIC| ≤ 10) were retained for further quantification of introgression proportions. From these, we extracted the Mixprop₂ parameter (representing introgression contribution, range 0–1) using a mixed model. Genome-wide introgression fractions were derived by weighting Mixprop₂ values against triplet frequencies across gene trees. Branch-specific introgression estimates were calculated as averages of Mixprop₂ for the triplets associated with each node. We used the R package quiblR (Edelman et al., 2019) to visualize the relative contributions of ILS versus introgression at each node as pie charts

2.6.4. Detection and quantification of gene flow

To detect and quantify gene flow (hybridization/introgression) between species, Dsuite v.0.3158 (Malinsky et al., 2021) was used to calculate D-statistics and f4-ratios for all possible species triplets. The analysis was performed using the Dtrios function with input files, including a compressed VCF file containing SNP data derived from the filtered nuclear gene alignment and the ASTRAL species tree as a guide topology. The D-statistics and f4-ratios results were visualized using the plot_d.rb and plot_f4ratio.rb scripts. The output from Dsuite Dtrios (sample_tree.txt) was then used as a direct input for the Dsuite F-branch to calculate the F-branch statistic. Finally, the resulting F-branch statistics were visualized using the dtools.py script.

2.6.5. Reconstruction of species reticulate evolution models

To investigate specific reticulate evolutionary events, we used species networks by applying the quartets (SNaQ) method from the Julia package PhyloNetworks v.0.6 (Bezanson et al., 2017; Solís-Lemus and Ané, 2016; Solís-Lemus et al., 2017) to construct phylogenetic networks. The taxon set for this computationally intensive analysis was strategically selected to represent the major phylogenetic clades and geographic regions, and to include putative hybrid lineages implicated in prior gene flow analyses. Owing to computational limitations, the dataset was reduced to 22 Lycium species and Atropa belladonna as outgroup. We then aligned and trimmed nuclear genes for these 23 taxa using trimAl v.1.2 (Capella-Gutiérrez et al., 2009) with the parameters gt 0.2, st 0.001, cons 85. A subsequent nuclear ML tree was constructed using IQ-TREE v.2.1.4 (Minh et al., 2020), and the coalescent species tree was estimated using ASTRAL-Ⅲ v.5.7.8 (Zhang et al., 2018). Species network inference was conducted based on the genes and species trees. We performed an iterative search for the optimal network, starting from the species tree (hmax = 0), and incrementally tested models with increasing numbers of reticulation events up to hmax = 5. Ten independent runs were conducted for each Hmax value. The optimal reticulate evolution model was selected by comparing the log-likelihood values (-ploglik). The final phylogenetic networks were visualized using Dendroscope v.3.7.6 (Huson and Scornavacca, 2012).

2.7. Divergence time estimation

The divergence times were estimated using BEAST v.1.10.4 (Drummond et al., 2012). The SortaDate package (Smith et al., 2018) was used to select the top 50 best-scoring SCN gene trees from the 378 nuclear gene trees. SortaDate was selected to optimize the dataset for molecular dating by prioritizing clock-like genes, had minimal topological conflicts with the species tree, and contained informative branch lengths (Smith et al., 2018; Stubbs et al., 2020; McLay et al., 2023). Due to the absence of Lycium fossils, three calibration points were applied: (T1) stem age of Solanoideae: 61.28 Ma (lognormal distribution, mean = 0.01, SD = 0.5, offset = 61.28 Ma) (Huang et al., 2023); (T2) stem age of Physalideae: 52.22 Ma (lognormal distribution, mean = 0.01, SD = 0.5, offset = 52.22 Ma) (Wilf et al., 2017); and (T3) crown age of Datura stramonium: 3.6 Ma (lognormal distribution, mean = 0.01, SD = 0.5, offset = 3.6 Ma) (Velichkevich and Zastawniak, 2003). An uncorrelated lognormal relaxed clock (Drummond et al., 2012), the GTR + I + G substitution model, and a birth-death process tree prior (Stadler, 2009) were used for the estimation. Markov chain Monte Carlo (MCMC) analysis was performed for 500 million generations, sampling every 1000 generations. The effective sample size (ESS) (> 200) was assessed using Tracer v.1.7 (Rambaut et al., 2018). TreeAnnotator v.2.4.4 (Drummond et al., 2006) was used to generate a maximum clade credibility (MCC) tree (discarding 0.5% burn-in). Node ages and their 95% highest density distribution (HPD) are shown in FigTree v.1.4.4.

2.8. Diversification rate analysis

Bayesian analysis of macroevolutionary mixtures (BAMM v.2.5.0) was performed to infer the speciation, extinction, and net diversification rates of Lycium based on the BEAST-derived chronogram (MCC tree from the 50-gene dataset) (Rabosky et al., 2014). The MCMC analysis was configured for 50 million generations with sampling every 10,00 generations. After discarding the initial 10% burn-in, results were analyzed using BAMMtools v.2.17 (Rabosky et al., 2014). Convergence was verified by ensuring an ESS > 200 for log-likelihood values.

2.9. Ancestral area and ancestral character reconstructions

We used the Bayesian binary Markov chain MCMC (BBM) methods in RASP v.4.2 (Yu et al., 2020) to reconstruct the ancestral ranges of Lycium species based on the BEAST-derived chronogram (MCC tree set from 50 genes). Due to limited sampling and uncertainty regarding the root area of the outgroup, it was removed from the biogeographic analyses. Distribution data for Lycium species were compiled from the Chinese Virtual Herbarium (CVH; https://www.cvh.ac.cn/), Global Biodiversity Information Facility (GBIF; https://www.gbif.org/), and International Plant Names Index (IPNI: https://www.ipni.org/). Following the floristic regionalization by Liu et al. (2023) and Lycium distribution patterns (Fukuda et al., 2001; Levin and Miller, 2005), we defined six biogeographic distribution areas: (A) North America, (B) Pacific Islands, (C) South America, (D) South Africa, (E) Saharo-Arabia, and (F) Eurasia. The Saharo-Arabian region was analyzed as a distinct unit to explicitly test its historical role as a potential corridor or independent area (following Liu et al., 2023). BBM analysis was run for 500,000,000 generations using 100 MCMC chains with the JC + G model.

Ancestral reconstruction of fruit color was performed using the R package phytools v.2.4.4 (Revell, 2024), based on the BEAST-derived chronogram (MCC tree set from the 50 genes). Information on fruit color per species (Table S1) was retrieved from online sources, including Flora of China (https://www.iplant.cn/foc), Flora of North America (both http://www.efloras.org), and Global Plants (https://plants.jstor.org/), and supplemented by our own field observations.

3. Results 3.1. Characteristics and comparative genomics of the Lycium ruthenicum mitochondrial genome

The complete mitochondrial genome of Lycium ruthenicum de novo assembly yielded a circular DNA molecule with a total length of 408,252 bp and a GC content of 45.4% (Fig. S1). The annotation revealed that the genome contains 65 genes, comprising 38 protein-coding genes, four rRNA genes, and 23 tRNA genes (see Table S5).

Repeat sequence analysis identified 33 simple sequence repeats (SSRs), the vast majority of which were mononucleotide repeats (27) with a distinct A/T bias (Fig. S2). Furthermore, 50 dispersed repeats (total length 24,672 bp) and 15 tandem repeats were detected, all of which were located in noncoding regions (Fig. S2). Codon usage bias analysis indicated a preference for synonymous codons ending in A or U in the mitochondrial protein-coding genes of Lycium ruthenicum (30 codons with RSCU > 1) (Fig. S3). Notably, the diversity of mitochondrial plastid DNAs (MTPTs) analysis identified 31 shared sequence fragments (total length 16,451 bp), accounting for approximately 4.03% of the mitochondrial genome (Fig. S4), which may have contributed to their evolutionary dynamics.

3.2. Nuclear phylogenomic analyses

Nuclear phylogenomic analysis based on concatenation (ML) (Figs. 1B, S5 and S6) and coalescence (ASTRAL species tree) (Figs. 2C and S5) strongly supports the monophyly of Lycium and its sister relationship with Nolana paradoxa (bootstrap percentage (BS) = 100, LPP = 1). Within Lycium, the nuclear concatenated ML tree and nuclear coalescent species tree recovered five well-supported clades (clades Ⅰ–Ⅴ), exhibiting largely congruent backbone topologies, with only slight differences among some nodes at the species level (Fig. S5). Among the six North American species (L. andersonii, L. californicum, L. fremontii, L. pallidum, L. parishii, L. torreyi), clade Ⅰ was strongly supported as monophyletic (BS = 100, LPP = 1) and sister to all remaining Lycium species (BS = 100, LPP = 1). Clades Ⅱ and Ⅲ formed a strongly supported monophyletic group (BS = 100, LPP = 1) containing seven South American species (L. ameghinoi, L. cestroides, L. chilense, L. ciliatum, L. cuneatum, L. morongii, L. nodosum). Clade Ⅳ included 19 species distributed across North America, the Pacific Islands, and South America. However, the internal topology of this clade differed between the ML and coalescent-species trees. The former resolved two primary subclades within clade Ⅳ: one consisting solely of nine North American species (BS = 100), and the other comprising four North American species (L. carolinianum, L. rachidocladum, L. richii, L. tenuispinosum), five South American species (L. berlandieri, L. carinatum, L. elongatum, L. exsertum, L. megacarpum), and the Pacific islands species L. sandwicense. In contrast, the coalescent species tree subdivided Clade Ⅳ into three distinct subclades: the Pacific Islands species L. sandwicense was resolved as sister to the rest of Clade Ⅳ; a second subclade contained the four North American and five South American species; and a third subclade comprised the nine North American species. Clade Ⅴ included species from South Africa, Saharan Arabia, and Eurasia. Species from South Africa (L. acutifolium, L. austrinum, L. ferocissimum) and Saharo-Arabia (L. schweinfurthii) formed a strongly supported monophyletic group (BS = 100, LPP = 1) as sister to the Eurasian species. Within the Eurasian subclade, both trees resolved L. ruthenicum as sister to all other Eurasian species (BS = 100, LPP = 1), whereas the placement of L. dasystemum and L. qingshuiheense differed between the ML and coalescent tree species (Fig. S5).

3.3. Phylogenetic conflicts among nuclear, plastid, and mitochondrial genomes

To investigate the cytonuclear conflicts and their underlying causes, we reconstructed a mitochondrial phylogeny (Fig. 1A) based on the newly assembled mitochondrial genome of Lycium ruthenicum and its orthologs, thereby providing a third independent perspective on the evolutionary history of Lycium. This mitochondrial tree was then compared with the nuclear ML tree (Fig. 1B) from this study and the plastid genome ML tree (Fig. 1C) from our previous study (Yisilam et al., 2025). All three phylogenies strongly supported the monophyly of Lycium (BS = 100) and consistently identified the North American clade of the six species as sister to the rest of Lycium (BS = 100). However, there were at least five significant phylogenetic discordances between these nuclear, plastid, and mitochondrial ML trees (Fig. 1): (1) the nuclear ML tree resolved Lycium into five major, well-supported clades, while both the plastid, and mitochondrial ML tree recovered only three major clades (clades Ⅰ–Ⅲ); (2) within the North American clade Ⅰ, internal relationships among the six species differed substantially between the three datasets; (3) both the plastid and mitochondrial phylogenies grouped species from South America, North America, and the Pacific Islands into a single major clade Ⅱ, which itself comprised several regionally defined subclades whose branching order conflicted between the two organellar genomes; however, the nuclear phylogeny subdivided this assemblage into three distinct lineages, with clades Ⅱ and Ⅲ comprising only South American species and clade Ⅳ containing a mixture of South American, North American, and Pacific Islands taxa; (4) regarding the Eastern Hemisphere species, corresponding to nuclear clade Ⅴ and plastid and mitochondrial clade Ⅲ, all three genomes supported their monophyly but conflicted in their substructure. Both nuclear and mitochondrial phylogenies supported a monophyletic group containing four South African species and the Saharo-Arabian L. schweinfurthii, which was sister to a monophyletic clade of Eurasian species. Conversely, the plastid tree depicted these South African species as paraphyletic to L. schweinfurthii, while also recovering a monophyletic Eurasian clade. (5) Within the Eurasian clade, the phylogenetic positions of L. dasystemum, L. qingshuiheense and L. cylindricum were incongruent among the three genome trees. Notably, all major conflicting nodes in the nuclear, plastid, and mitochondrial phylogenies had high bootstrap values (Fig. 1), highlighting robust but discordant signals from the two genomes. The dynamics of this conflict can be quantitatively elucidated through the subsequent analyses of ILS and gene flow.

3.4. Patterns and quantification of phylogenetic discordance 3.4.1. Visualization of topological conflict patterns

The MDS plot (Fig. 2A) demonstrated that most individual nuclear gene trees displayed average bootstrap values ranging from 50% to 87%. These values were widely distributed across the measured space, indicating a diverse representation of the data. Although some gene trees contained strongly supported nodes, none of the individual gene trees fully recovered the topology of the tree species, and all exhibited a clear distance from it in the MDS space, which intuitively reflects the extensive topological heterogeneity among gene trees (Fig. 2A).

Quartet-sampling (QS) analysis strongly supported the monophyly of Lycium based on the nuclear gene tree (QC = 1, QI = 1; Fig. S7). While most internal branches showed positive concordance (mean QC > 0), significant discordance was detected at several key nodes, as evidenced by negative QC values. The QD values for most of the nodes approached zero, suggesting a single predominant alternative topology among the conflicting quartets. Specifically, within the fully resolved clade Ⅰ (QS: 1/NA/1), sister relationships of L. schaffneri, L. puberulum, L. cooperi, and L. macrodon showed significant conflict (QC = −0.23, QD = 0.17). The placement of clade Ⅱ as sister to the rest of the tree received moderate support (QC = 0.68, QD = 0.12). Clade Ⅲ was moderately supported as sister to clades Ⅱ and Ⅳ (QC = 0.35, QD = 0.97, respectively). A notable and strongly conflicting signal was found for the sister relationship between clades Ⅴ and Ⅳ (QC = −0.23, QD = 0.4), indicating substantial topological discordance. Within clade Ⅴ, the divergence between South African and Eurasian species lineages was moderately supported (QC = 0.40, QD = 0.19) (Fig. S7).

The PhyPart analysis (Fig. 2C) demonstrated that 336 of the 378 informative nuclear gene trees supported the monophyly of Lycium and 307 supported the monophyly of clade Ⅰ. Significant conflict existed within clades Ⅱ–Ⅳ. However, some branches received strong support from the trees. Monophyly of the Eastern Hemisphere clade was supported by 76 gene trees and 177 Eurasian species (Fig. 2C). These pervasive conflicts suggest significant contributions from ILS, introgression, or both.

3.4.2. Assessment of incomplete lineage sorting

The MSCquartets analysis results revealed that under varying significance thresholds (α = 0.01 to 1e−06), most quartet concordance factors clustered near the region representing the dominant tree topology in the simplex plot (Fig. 2B), indicating that ILS is a primary factor driving gene tree discordance within Lycium. Concurrently, a substantial proportion of the quartets rejected the pure multispecies coalescent model (represented by red triangles in Fig. 2B), indicating that additional evolutionary processes, notably gene flow, also played a significant role.

3.4.3. Quantifying the relative contributions of ILS and introgression

We performed QuIBL analysis on 37,023 triplets, of which 35.3% (13,093 triplets) strongly supported the pure ILS model (dBIC > 10), while only 11.3% (4,190 triplets) supported the ILS + introgression model (dBIC < -10). At multiple points throughout the phylogeny, ILS has been shown to account for a significant proportion of the gene tree-species tree discordance. Many nodes displayed ILS contributions exceeding 40%, with several nodes reaching notable peaks above 60%, specifically recorded at 61.00% (Fig. 3). The node-level ILS contribution ranged from 2.53% to 61.00 %, with an average of 32.03%, whereas the average introgression contribution was 3.58% (range: 0.28–8.55%). Although most species pairs showed low introgression levels, significantly elevated gene flow was detected in clade Ⅱ, comprising L. ameghinoi, L. ciliatum, and L. chilense, and all six species derived from North American radiation: L. andersonii, L. californicum, L. fremontii, L. pallidum, L. parishii, and L. torreyi (Fig. S8). Collectively, these results confirm that ILS is the primary driver of gene tree discordance in Lycium.

3.4.4. Detection of historical gene flow

ABBA-BABA (D-statistic) tests revealed widespread and significant evidence of historical gene flow. Of 12,341 tests, 55.14% (6805) showed significant signals (P < 0.05), with 34.7% (4281) being highly significant (|Z| > 3, P < 0.001) (Fig. 4A). Notably, some of the strong signals were detected among several Asian species. For example, the trios L. dasystemum, L. barbarum, and L. ruthenicum showed highly significant results (D = 0.782, Z = 20.162, P < 2.3e-16), indicating exceptional genetic exchange between L. dasystemum and L. ruthenicum. The f4-ratio analysis provided quantitative estimates of introgression proportions, indicating that approximately 29–31% of the L. sandwicense genome was derived from a specific donor lineage (e.g., f4-ratio = 0.315 relative to L. tetrandrum, Z = 7.029, P = 2.08e-12) (Fig. 4B). F-branch analysis further resolved the lineage-level introgression patterns (Fig. 4C), revealing strong signals both between and within continents. These included gene flow between L. austrinum and the African lineages (represented by L. tetrandrum, L. ferocissimum, L. acutifolium), Eurasian L. dasystemum and its lineage members (L. barbarum and L. chinense), South American L. elongatum and African L. ferocissimum, and North American L. schaffneri and its relatives (L. pallidum and L. shockleyi).

3.4.5. Phylogenetic network inference

Model comparisons indicated that a network incorporating three reticulation events (h = 3) best explained the phylogenetic conflicts observed in the data (Fig. 4D). The optimal network model (Fig. 4E) consistently resolved three key historical gene-flow events across multiple independent runs. First, a major introgression event was inferred from an unsampled or extinct ancestral lineage distributed in Africa and the Saharo-Arabian region into the Eurasian L. dasystemum, contributing approximately 87.4% of its genetic composition, whereas the ancestral lineage of the widely distributed Eurasian L. barbarum and L. chinense contributed the remaining 12.6%. This event coherently explains the extreme D-statistic signals detected in L. dasystemum, L. barbarum, and L. ruthenicum (D = 0.782, Z = 20.162). The second reticulation event formed the backbone of the core evolutionary network in Lycium. The descendant lineage of this event gave rise to a large clade containing species such as L. sandwicense and L. schweinfurthii. The genomic composition of this lineage was derived from two sources: 42.6% was contributed by the South American species L. elongatum, and 57.4% originated from a composite lineage that included an unsampled or extinct ancestor together with L. ruthenicum. This event quantifies the substantial South American genetic contribution to L. sandwicense previously suggested by f4-ratio analysis. The final reticulation event occurred within the African-Saharo clade, indicating that L. ferocissimum has experienced substantial historical introgression, with its genome derived approximately 64.7% from L. acutifolium and 35.3% from the lineage represented by L. austrinum. This event, consistently identified in all models with h ≥ 2, underscores the role of gene flow in the diversification of African Lycium. Furthermore, the network model supported the grouping of L. pallidum, L. schaffneri, and L. shockleyi as a recently and rapidly diverging North American clade, which was not directly involved in the major reticulation events described above.

3.5. Divergence time estimation

The BEAST-derived MCC chronogram, based on the 50 top SCN dataset (Fig. 5), estimated the Lycium and N. paradoxa at 35.83 Ma [node 0, 95% highest posterior density (HPD): 26.01–45 Ma]. The most recent common ancestor (MRCA) of Lycium viz., its crown age, was estimated at 21.84 Ma [node 1, 95 % HPD: 15.67–28.64 Ma]. Divergence between North American and South American species occurred at 18.70 Ma [node 2, 95% HPD: 1.83–22.05 Ma]. Subsequently, at 14.83 Ma [node 3, 95% HPD: 10.49–19.58 Ma], the Lycium lineages of the American continent diverged from those of the Eastern Hemisphere (South Africa, Saharo-Arabia, and Eurasia). The divergence of Eurasian species from the rest of the Eastern Hemisphere occurred at 9.52 Ma [node 4, 95% HPD: 6.25–13.25 Ma] (Fig. 5).

3.6. Species diversification rates

The BAMM analysis results indicated that the genus Lycium experienced time-dependent rate changes (Fig. 6A); however, no significant discrete rate-shift events were detected (Fig. 6B). Specifically, both the speciation rate and net diversification rate increased gradually from the crown age (~21.84 Ma) until about 4 Ma, after which they stabilized, while the extinction rate declined over the same period and reached a stable state after approximately 4 Ma (Fig. 6A). This pattern indicates that lineage accumulation in Lycium was relatively rapid during the early evolutionary stage, and then gradually slowed and stabilized.

3.7. Ancestral area reconstruction and ancestral character reconstructions of fruit colors

The BBM ancestral area reconstructions showed that the genus Lycium originated in North America (region A, node 1; probability = 52.17%) (Fig. 5). Subsequently, species from North America (A) dispersed throughout South America (C) (probability = 50.67%). Following this isolation, South America (C) rapidly became the core area for accelerated species diversification within Lycium (nodes 2 and 3; probability = 97.66%). Subsequently, from this South American center, lineages dispersed to the Pacific Islands (B) and back to North America (A) (Nodes 4 and 5). We inferred at least four dispersal events between North America (A) and South America (C). Furthermore, South American species dispersed across the Atlantic to South Africa (D) (node 6, probability: 86.56% to 34.94%). Finally, from South Africa, the lineages reached Eurasia (F, node 7) via the Saharo-Arabian region (E) (Fig. 5).

Ancestral state reconstruction revealed that the fruit color of the most recent common ancestor of Lycium was most likely ‘red/orange’, and this state is preserved in most extant species (Fig. 6C).

4. Discussion

Understanding the evolutionary history of Lycium is fundamental to unraveling its adaptive radiation, biogeography, and speciation. While previous studies have established a broad pattern of intercontinental disjunction using limited genetic markers, the relative roles of key evolutionary processes such as incomplete lineage sorting (ILS) and hybridization in shaping phylogeny remain unclear. This study moves beyond descriptive biogeography to provide a process-oriented understanding of Lycium’s evolution. By integrating, for the first time in this genus, data from three genomic compartments, nuclear genes, mitochondrial, and plastid genomes (Yisilam et al., 2025), and applying a suite of phylogenomic methods, we achieve three key advances: (1) a tripartite assessment of phylogenetic discordance; (2) a quantitative partition of conflict causes, demonstrating that ILS is the predominant driver of intra-nuclear discordance, a signature of rapid Miocene radiation; and (3) the identification and modeling of three specific hybridization events critical to shaping major transcontinental lineages, via robust gene-flow detection and explicit network analysis. Collectively and within the framework of our sampling, our findings indicate that the evolutionary history of Lycium was shaped by a complex interplay of rapid radiation, ILS, hybridization/introgression, and long-distance dispersal (LDD), all set against the backdrop of Miocene geoclimatic upheavals.

4.1. A robust nuclear phylogenomic framework for Lycium

Previous research on the phylogenetic relationships within the genus Lycium relied on nuclear gene fragments (ITS) and single plastid genes (matK, trnL) along with their interspaces (trnT–trnL, trnL–trnF) (Fukuda et al., 2001; Levin and Miller, 2005; Levin et al., 2009). The present results show that both the concatenated ML tree and the ASTRAL species tree of nuclear genome (SCN) data (Fig. S5) strongly support the monophyly of Lycium (BS/LPP = 100/1). This is consistent with the results obtained for plastid genomes (Yisilam et al., 2025). The nuclear ML tree divided Lycium into five major clades (Ⅰ–Ⅴ) with high support and the topological structures were highly consistent with the coalescent tree species. This robust framework provides a foundation for critically reevaluating traditional classifications (Fig. S6) and to understand evolutionary processes.

Nuclear phylogeny indicated significant discordance with traditional classifications based on morphology (Fig. S6; Fukuda et al., 2001; Levin and Miller, 2005). The traditionally defined section, the Sclerocarpellum, which is characterized by hardened fruits (Miller, 2002), showed a polyphyletic distribution in the nuclear gene tree. Lycium ameghinoi clustered with South American species (clades Ⅱ/Ⅲ), and L. californicum was nested among North American taxa (clade Ⅳ) (Fig. S6). The remaining members of this section formed clades with L. pallidum and L. shockleyi (Clade Ⅰ), indicating that this trait is the result of convergent evolution during adaptation to arid environments (Fukuda et al., 2001). Species of the section Mesocape also showed a scattered distribution in the phylogenetic tree, supporting the hypothesis that their shared morphological traits likely represent retained ancestral characteristics or parallel evolution (Fig. S5; Miller, 2002; Levin and Miller, 2005).

The nuclear framework resolves long-standing controversies regarding relationships between species. The concatenated ML and the ASTRAL species tree were robustly supported (BS = 100, LPP = 1) L. chilense and L. ciliatum as sister species, forming a distinct monophyletic clade (Fig. 1, Fig. 2C). This high-resolution result confirms their close phylogenetic affinity and resolves previous taxonomic uncertainties regarding their relationships (Levin and Miller, 2005; Levin et al., 2009). Furthermore, within clade Ⅴ (encompassing species from South Africa, Saharan Arabia, and Eurasia), L. schweinfurthii (Saharan Arabia) clustered robustly with South African species such as L. austrinum, L. acutifolium, L. ferocissimum, L. tetrandrum. This phylogenetic placement suggests a South African origin or a deeply shared ancestry for this Saharo-Arabian/South African group. This is supported by the distribution of L. schweinfurthii in the eastern Mediterranean and Saharan Arabia, as well as the vast range of its close relative, L. shawii. Although L. shawii was not sampled in our study, its broad distribution, spanning northern and southern Africa, Mediterranean Europe, the Middle East, and Asia (Levin et al., 2007), provides a critical biogeographic link between South Africa and Eurasia. In the Eurasian group, L. ruthenicum was sister to the rest of the clade, which is consistent with its representation as an ancient lineage in the region (Zhang et al., 2024; Yisilam et al., 2025). Based on single individuals per species, these findings provide a crucial molecular foundation for future population-level studies to validate species boundaries and integrate morphological and ecological evidence.

4.2. Driving factors of phylogenetic conflicts and reticular evolutionary history of Lycium species

By integrating signals from three independent genomic compartments, this study revealed two types of significant phylogenetic conflicts in the genus Lycium: (1) tripartite topological conflicts among nuclear, mitochondrial and plastid genome trees (Fig. 1), and intra-nuclear discordance (gene trees vs. coalescent species tree) (Fig. 2C). These conflicts are driven primarily by ILS and historical hybridization or introgression. Our analysis disentangles the primary causes of these two distinct types of conflict.

First, the tripartite discordance among the nuclear, mitochondrial, and plastid genomes stems from their distinct genetic and evolutionary modes. The nuclear genome, with biparental inheritance and frequent recombination, rapidly accumulates mutations and records recent divergent events (Dong et al., 2021). By contrast, uniparentally inherited plastid and mitochondrial genomes lack sexual recombination and evolve at different rates (Liu et al., 2018; Wang et al., 2024a; Zou et al., 2025). For instance, the mitochondrial genome of L. ruthenicum exhibits features such as pronounced AT-biased codon usage, highly repetitive sequences, and significant plastid DNA insertions, which may influence its evolutionary trajectory. Moreover, during the early rapid diversification of Lycium, ILS likely caused incomplete sorting of ancestral polymorphisms in organellar genomes, further complicating their phylogenetic signals (Dobrogojski et al., 2020). Critically, the strong observed cytonuclear discordance was most parsimoniously explained by historical hybridization/introgression. This was strongly supported by our genome-wide D-statistics and explicit phylogenetic network models. Although ILS during rapid radiation may have contributed, the clear and localized signals of gene flow between key lineages provide a more direct explanation for the topological conflicts between the biparentally inherited nuclear genome and uniparentally inherited organellar genomes. Such cytonuclear discordance, well documented in numerous angiosperm families and genera, such as Vitaceae (Liu et al., 2021), Eucalyptus (McLay et al., 2023), Catalpa (Dong et al., 2022), and rapidly radiating lineages (e.g., Petunia, Solanaceae) (Pezzi et al., 2024), clearly reflects the heterogeneous histories captured by different genetic systems.

Previous studies have suggested the possibility of a hybrid origin for some South African tetraploid species (e.g., Lycium villosum) based on cytological and morphological evidence (Venter et al., 2003). Our integrated analyses provide robust genomic evidence of widespread gene flow in Lycium. D-statistics indicated significant signals in over half of the taxon combinations, and f4-ratio analysis quantified substantial introgression, such as 29–31% of the L. sandwicense genome derived from a specific donor lineage (Fig. 4A and B). F-branch analysis revealed cross-continental introgression patterns (e.g., between the South American and African lineages) (Fig. 4C). To integrate these signals, we constructed a phylogenetic network (Fig. 4E) that inferred three key hybridization/introgression events, indicating that the Eurasian L. dasystemum largely (87.4%) derives from an African/Saharo-Arabian ancestor, explaining strong D-statistic signals, and quantifies the South American contribution to L. sandwicense, corroborating f4-ratio estimates. This profound African/Saharo-Arabian ancestry in L. dasystemum suggests that its establishment in Eurasia involved not just dispersal, but also a major episode of genetic exchange. This introgression likely occurred during the late Miocene, when intensified aridification could have expanded contiguous dry habitats, facilitating contact between African and incipient Eurasian lineages across the Saharo-Arabian region. These events substantiate that historical hybridization is the primary force behind the major cytonuclear discordance.

Second, we observed extensive discordance within the nuclear genome. Multiple lines of evidence indicate that ILS is the dominant factor in this intranuclear conflict (Fig. 2, Fig. 3). The broad distribution of gene trees in multidimensional scaling (MDS) space visually reflects their topological heterogeneity (Fig. 2A). QuIBL analysis showed that ILS explained a substantial portion of the observed discordance, with contributions reaching peaks of 61.00% at several nodes and exceeding 40% at many nodes across the phylogeny (Fig. 3), whereas the average introgression contribution was much lower (3.60%). This high ILS level is attributable to the rapid early radiation of Lycium, as supported by the BAMM analysis, which detected a significant increase in the net diversification rate near the crown node (21 Ma; Fig. 5A and B), which aligns with the period of global cooling and increased aridity following the Oligocene-Miocene transition. The short time intervals between successive lineage splits during this period hindered the complete sorting of ancestral polymorphisms, thereby promoting ILS. Although introgression showed localized significance between the intercontinental lineages (Fig. S8), the overall contribution to tree species discordance across the nuclear dataset was relatively limited. Thus, the rapid Miocene radiation provides an evolutionary context for the prevalence of ILS in shaping nuclear gene tree incongruence.

The spatial and temporal patterns of these key hybridization events (e.g., involving African-Saharo-Arabian and Eurasian lineages) coincide with the major phases of intercontinental dispersal and aridification discussed in the following biogeographic history. In summary, hybridization/introgression exacerbates phylogenetic discordance through two complementary mechanisms in Lycium: (1) the distinct inheritance patterns of nuclear and organellar genomes directly create topological conflicts between their trees, explaining cytonuclear discordance; and (2) introgressed genomic segments alter local gene genealogies, contributing to the complex patterns of intra-nuclear gene tree variation. Such introgression-driven perturbations are common in angiosperms, as seen in Lachemilla (Morales-Briones et al., 2018), Malpighiales (Cai et al., 2021), Diabelia (Ke et al., 2021), Vitaceae (Ma et al., 2021), Polygonatum (Qin et al., 2024), and Morus (Wang et al., 2024b), highlighting the universality of this mechanism.

4.3. Biogeographic history of Lycium

Our robust nuclear phylogenomic framework and biogeographic analyses revealed that the evolutionary history of Lycium is closely linked to key geoclimatic and paleogeological events that govern its origin, radiation, and intercontinental dispersal. Molecular dating and ancestral area reconstructions (Fig. 5) show that Lycium originated in the Late Eocene/Early Oligocene, with its crown group beginning to diversify in North America during the early Miocene (21.84 Ma; 95% HPD: 17.72–30.34 Ma). A significant increase in the diversification rate began at 21 Ma (Fig. 6B), which corresponds to a period of global cooling and increased aridity following the Oligocene-Miocene transition (Zachos et al., 2001). We infer that this climate shift created new ecological opportunities that drove the initial radiation of Lycium in North America. However, the interpretation of this BAMM result requires caution due to the limited genus-level sampling coverage, as undersampling may affect the estimation of rate heterogeneity (Rabosky et al., 2017).

Our divergence time estimates for Lycium are consistent with the wide confidence interval of an earlier gene fragment-based estimate (29.4 ± 9.7 Ma; Fukuda et al., 2001), while slightly older (6 Ma) than those obtained from nuclear genes (Huang et al., 2023) and substantially older than those based on plastid markers (Yisilam et al., 2025; Särkinen et al., 2013). These discrepancies can be attributed to methodological differences, including our use of early fossil calibration (Wilf et al., 2017), as well as inherent disparities in evolutionary rates and lineage sorting histories between nuclear and plastid genomes. The older ages inferred from our nuclear dataset likely reflect the true species divergence times better, as the larger effective population size and biparental inheritance of nuclear genes make them less susceptible to complete lineage sorting and capture more ancient genealogical signals than the organellar genomes. Notably, such marked discrepancies between nuclear and plastid-derived divergence times have also been reported in other solanaceous groups (e.g., Solanum; Messeder et al., 2024). Our use of multiple nuclear genes and careful calibration provided a robust chronological framework for biogeographic inferences. Notably, while our primary calibration points were necessarily external to Lycium, the estimated crown age (~21 Ma) and its broad 95% HPD interval (17.72–30.34 Ma) consistently support an early Miocene origin, suggesting that this central conclusion is not unduly sensitive to the specific implementation of these secondary calibrations. Future studies that combine genomic-scale data with improved fossil calibrations could further refine these estimates.

Our results further indicate that shortly after its origin in North America (~21 Ma), Lycium dispersed into South America (18.7 Ma). The resulting amphitropical disjunction between the arid regions of both continents conforms to the classic American Amphitropical Disjunction (AAD) pattern (Simpson et al., 2017). This pattern is most parsimoniously explained by direct long-distance dispersal (LDD) (Raven, 1963; Simpson et al., 2017), rather than by the younger Central American corridor, which is strongly supported by the small bird-dispersed fruits of Lycium. Subsequently, South American Lycium underwent diversification between 18 and 10 Ma (Fig. 5), likely driven by heterogeneous canyon plateau habitats created by the main Andean orogeny (10–15 Ma) (Hoorn et al., 2010). The orogeny created extensive rain shadows and a complex mosaic of arid and intermontane basins and valleys along the western side of the mountains (Pérez-Escobar et al., 2022), facilitating allopatric speciation in these newly formed and heterogeneous arid environments. This diverse South American clade subsequently served as a source of trans-Pacific dispersal and of multiple dispersal events back to North America (at least four times). Together, these formed a dispersal pattern connecting the American continent and the Pacific Islands. These dispersal processes are strongly supported by the adaptive traits of Lycium, such as its small, brightly colored fruits (e.g., red, orange, or black-purple), which facilitate bird-mediated dispersal (Fig. 6C).

Transatlantic dispersal is another key event in the biogeographic history of Lycium, which is consistent with earlier phylogenetic hypotheses (Fukuda et al., 2001). Our nuclear ancestral area reconstructions indicate that, during the Middle Miocene (14.83 Ma, 95% HPD: 10.49–19.58 Ma), South American species dispersed over-sea to South Africa (Fig. 5). This event occurred long after the breakup of Gondwana (96–105 Ma; McLoughlin, 2001), which firmly supports the long-distance dispersal hypothesis (Raven and Axelrod, 1974). Similar transoceanic LDD patterns have been documented in other plant groups, such as dispersal from South America to Africa in the Miocene within the tribe Ranunculeae (Ranunculaceae) (Emadzade and Hörandl 2011), further underscoring the importance of trans-Atlantic dispersal in shaping intercontinental distribution patterns.

After arriving in the Eastern Hemisphere around the early Late Miocene (7–10 Ma), intensified aridification in South Africa, coupled with the expansion of the Saharan Desert, promoted the development of a continuous belt of arid-adapted plants. This, in turn, facilitated a range shift that triggered species differentiation and northward range expansion within Africa (see also Calviño et al., 2016). This aridification also promoted the diversification and dispersal of the southern African genus Thamnosma (8.53 Ma; 95% HPD: 5.28–12.11 Ma; Thiv et al., 2011). Thus, based on our ancestral area reconstruction results, we inferred that during the Late Miocene period of increasing aridity, Lycium may have expanded its range northward from Southwest Africa, accompanied by species differentiation within the continent. Our nuclear (6.66 Ma) and plastid 7.88 Ma (Yisilam et al., 2025) genomic data suggest that Late Miocene aridification was a pivotal driver of local diversification in Lycium. The inference that arid corridors facilitated expansion finds support in later periods, wherein the Pleistocene African Dry Corridor (ADC) explains disjunctions in Senecio flavus (Milton et al., 2022). Our earlier Lycium event implies that similar preexisting arid corridors functioned as conduits for range expansion during the Late Miocene.

Shortly thereafter, around 5.69 Ma (Fig. 5), South African Lycium dispersed into the Saharo-Arabian region. These dispersal events likely occurred via the Levantine corridor, a well-documented land connection between Africa and Arabia (El Zaatari, 2018), which has received considerable attention from biologists because of its role as a conduit for shared evolutionary history (e.g., Portik and Papenfuss, 2015; Meheretu et al., 2024). This Saharo-Arabian lineage subsequently served as a biogeographic steppingstone for the dispersal of Lycium into Eurasia, although the possibility of bird-mediated long-distance dispersal cannot be ruled out.

The rapid rise of the Qinghai-Xizang Plateau (QXP) during the Miocene and the attainment of its present height before the Pliocene profoundly shaped the evolution of temperate flora across the Northern Hemisphere, with pronounced impacts in Eurasia (Wang et al., 2008; Zhang and Fritsch, 2010; Xu et al., 2010; Liu et al., 2014; Meng et al., 2015). This significant geological event fundamentally altered the regional topography and climate systems, driving major shifts in plant distribution, speciation, and dispersal routes, as evidenced in lineages such as Sorbus (Li et al., 2017) and Glycyrrhiza (Duan et al., 2020). Specifically, our results indicate that the rapid uplift of the plateau during the Late Miocene (8 Ma) promoted the evolution and diversification of L. ruthenicum at approximately 6.23 Ma, marking the successful establishment and initial specialization of this lineage in Eurasia (Fig. 5). These findings are consistent with those of previous studies (Zhang et al., 2024; Yisilam et al., 2025). Furthermore, during the Pliocene to Quaternary, extreme aridification in inland Eurasia, particularly in northwestern China (Guo et al., 2002), driven by the most recent rapid uplift of the QXP (3.4–1.8 Ma), and the Quaternary glacial–interglacial cycles of the Tianshan Mountains significantly promoted further diversification within the genus Lycium. This evolutionary pattern has also been documented in multiple arid-adapted plant lineages, including Gymnocarpos przewalskii (Ma and Zhang, 2012), Hexinia polydichotoma (Su et al., 2012), Zygophyllum xanthoxylon (Shi and Zhang, 2015), as well as in other taxa, such as Populus euphratica (Zeng et al., 2018), Malus sieversii (Zhang et al., 2021), and L. ruthenicum (Yisilam et al., 2022).

5. Conclusion

This study advances the understanding of Lycium evolution by moving beyond singular genomic perspectives to integrate and compare nuclear, plastid, and mitochondrial phylogenetic signals. We constructed the most robust phylogenetic framework to date, using hundreds of nuclear genes that resolved deep relationships and revealed significant cytonuclear discordance. Crucially, for the first time in Lycium, we quantified the relative contributions of ILS and hybridization to this conflict. Our analyses demonstrated that ILS were the predominant force (> 60%) behind gene tree discordance within the nuclear genome, a consequence of rapid Miocene radiation. Simultaneously, by applying a suite of phylogenetic networks and gene flow detection methods, we identified and modeled three specific hybridization events that shaped major transcontinental lineages, offering a concrete historical narrative for previously observed discordance. Biogeographic reconstructions anchored by nuclear gene dating reaffirm a North American origin and detail the sequence of long-distance dispersal enabled by bird-dispersed fruits. While the broad intercontinental disjunction pattern aligns with earlier studies, our work provides a critical mechanistic layer for quantifying the ILS and modeling hybridization, which explains how this pattern emerged from a complex evolutionary process. While our sampling covered the major phylogenetic lineages and geographic regions, increasing representation from under-sampled areas (e.g., Australia and the Pacific Islands beyond Hawaii) in future studies will be crucial to resolve recent radiation and complex regional patterns fully. Analyses with expanded taxon sampling may provide further insights. The explicit hybridization models proposed herein present testable hypotheses for future population genomic studies within the implicated lineages.

Acknowledgements

This work was supported by the Key Research and Development Program of Guangxi Province of China (GuikeAB25069162, GuikeAB25069252), the National Natural Science Foundation of China (31760102, 32360058, 32570239), the Central Government Guides Local Science and Technology Development Projects (2023ZYZX1224), and Science and Technology Projects of Xizang Autonomous Region, China (XZ202402ZD0005). We would like to thank the University of Arizona Herbarium (ARIZ) and the Wisconsin State Herbarium (WIS) for providing Lycium specimens.

CRediT authorship contribution statement

Gulbar Yisilam: Methodology, Formal analysis, Software, Data curation, Writing – original draft, Writing – review & editing. Yu-Zhu Gao: Data curation. Zuan Wei: Data curation. Hans Peter Comes: Writing – review & editing. Zheng Li: Writing – review & editing. Pan Li: Conceptualization, Investigation, Resources, Supervision, Writing – review & editing. Xin-Min Tian: Conceptualization, Funding acquisition, Methodology, Investigation, Resources, Supervision, Writing – review & editing. All authors have read and agreed to the published version of the manuscript.

Data availability

The raw sequencing data from all 39 Lycium samples in this study have been submitted to the NCBI database under BioProject ID PRJNA1371369.

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.2026.02.005.