1. Introduction

Elaeagnaceae, a small family containing ca. 100 species of evergreen and deciduous trees and shrubs, mainly distributed in the Northern Hemisphere with a few species extending their distributions from southern Southeast Asia to northeast Queensland of Australia (Qin et al., 2007). The Elaeagnaceae family contains plants that are economically (Abdalla, 2019), ecologically (Ruan and Li, 2002; Zhang et al., 2018a; Wu et al., 2022a), and pharmacologically (Xu, 1994; Bartish and Swenson, 2004; Bartish and Thakur, 2022) important. However, phylogenetic relationships within the family remain unresolved and the factors that have influenced the diversification and biogeography of Elaeagnaceae have yet to be unexplored.

Elaeagnaceae is a monophyletic family (Bartish and Swenson, 2004; Sun et al., 2016) within Rosales, where it forms a clade with two small and narrowly distributed families, Barbeyaceae and Dirachmaceae (Soltis et al., 2000, 2011; Sytsma et al., 2002; Zhang et al., 2011). Elaeagnaceae includes three widely recognized monophyletic genera: Elaeagnus Tourn. ex L. (2n = 28), Hippophae L. (2n = 24), and Shepherdia Nutt. (2n = 22) (Cooper, 1932; Bartish and Swenson, 2004; Jia, 2013; Zhang, 2016; Jia and Bartish, 2018; Bartish and Thakur, 2022). Phylogenetic analyses based on molecular data (i.e., chloroplast DNA, mitochondrial DNA and nuclear DNA) have indicated that Elaeagnus is sister to the well-supported clade of Hippophae and Shepherdia (Sun et al., 2016; Jia and Bartish, 2018). However, intrageneric relationships within the Elaeagnaceae have yet to be well resolved. For example, intrageneric relationships within Hippophae have varied across several studies (Hyvonen, 1996; Swenson and Bartish, 2002; Bartish et al., 2002; Lian et al., 2003; Jia and Bartish, 2018). The intrageneric relationships within Elaeagnus, the largest genus of this family, also remain poorly resolved. Previous studies divided Elaeagnus into two sections, one consisting of deciduous species (sect. Deciduae Serv.) and one of evergreen species (sect. Sempervirentes Serv.) (Servettaz, 1909; Chang, 1983; Qin et al., 2007). Each of these sections was subsequently divided into two groups based on morphological traits (Sun and Lin, 2010). However, recent molecular analysis has indicated that these clades are non-monophyletic and that multiple reticulation events may have occurred during the rapid radiation of Elaeagnus (Zhang, 2016). Molecular phylogenetic analysis of the intrageneric relationships within Shepherdia, which only contains three species, has yet to be conducted.

Elaeagnaceae have long been of interest to biologists due to the disjunct pattern of distribution in the Northern Hemisphere. Elaeagnaceae fossils dating from the Cretaceous to the Miocene have been discovered in East Asia and North America (Becker, 1960; Grigoreva, 1969; Song and Qian, 1989; Su et al., 2014), indicating that this family might have been more widely distributed, either by gradual dispersal and vicariance or by multiple long-distance dispersals. However, few studies have explored the origin and biogeographical history of Elaeagnaceae or the factors that led to their current distributions. One biogeographical study has indicated that Hippophae originated in the Qinghai-Tibet Plateau (QTP) during the late Eocene and that its extant lineages diverged beginning in the late Oligocene/early Miocene before being dispersed toward Europe during the late Miocene (Jia and Bartish, 2018). To better understand the origin and historical biogeography of this family and its major clades, divergence time estimations must comprehensively sample taxa and use more fossil calibration points.

The aim of this study was to better understand the processes that have shaped the evolution and diversity of Elaeagnaceae. Specifically, we generated phylogenomic data, including nuclear loci and plastome sequences, to clarify phylogenetic relationships among the three genera of Elaeagnaceae and the interspecific relationships within each of the three genera. We also characterized cyto-nuclear discordance between phylogenetic trees to identify incomplete lineage sorting (ILS) and hybridization events within the evolutionary history of the family. Finally, we determined the origin, timing of diversification, and biogeographical history of Elaeagnaceae.

2. Materials and methods 2.1. Taxon sampling

Sampling of Elaeagnaceae in this study followed the nomenclature of Plants of the World Online (POWO, https://powo.science.kew.org/, 2023-12-12). The current sampling scheme covers all three accepted genera of Elaeagnaceae and represents 67/92 (= number of species sampled/number of species recognized) species of Elaeagnus, 5/7 species of Hippophae, and 3/3 species of Shepherdia. Ten Rhamnaceae species, one Barbeyaceae species and one Dirachmaceae species obtained from Tian et al. (2024) were used as outgroups in phylogenetic analyses based on a nuclear dataset. In addition, another 16 samples of five other families of Rosales and a sample of Fagaceae were also included in dating analysis (Tables S1 and S2). Leaf materials for the samples were obtained from the field or from herbarium specimens from CAS, KUN, NY, OS, and TEX (acronyms following Index Herbariorum). Additionally, the complete plastome sequences of two subspecies of Hippophae rhamnoides L. and two Rhamnus L. species were obtained from GenBank (Table S3).

2.2. Hybrid enrichment and genome skimming

For hybrid enrichment and sequencing (Hyb-Seq) (92 samples for Hyb-Seq are shown in Table S1), total genomic DNA was extracted from fragments of herbarium specimens or silica-dried material using a modified Cetyltrimethyl Ammonium Bromide (CTAB) protocol (Doyle and Doyle, 1987; Folk et al., 2021). We used an exonic bait set for 100 single-copy or low-copy genes (133, 433 bp in total) developed for phylogenetic analyses across the rosid clade (Folk et al., 2021; Fu et al., 2022), which includes Elaeagnaceae. This bait set was developed using 78 rosid transcriptomes and Arabidopsis thaliana (L.) Heynh. as a reference genome. Genomic DNA quantification, library preparation, target enrichment, and Illumina sequencing with 150-bp paired-end reads were conducted by Rapid Genomics (Gainesville, Florida, the United States). For genome skimming (63 samples are as shown in Table S2), total DNA was extracted from fresh or silica gel dried leaves or herbarium specimens using the CTAB method. cDNA libraries were constructed following Illumina's technical guide and sequenced on the Illumina HiSeq 2500 platform, which generated ~2 Gb of raw sequence data per sample. The construction of cDNA libraries and sequencing were performed at BGI (Shenzhen, China).

2.3. Nuclear DNA assembly, filtering, and alignment

Raw data from Hyb-Seq were first filtered to remove low-quality reads and adapters using Trimmomatic v.0.38 (Bolger et al., 2014) (Phred score = 33; ILLUMINACLIP: TruSeq3-PE-adapters.fa: 2:30:10:8:TRUE; SLIDINGWINDOW: 20:20). We then used HybPiper v.1.3.1 (Johnson et al., 2016; https://github.com/mossmatters/HybPiper/) to assemble the cleaned reads using 100 protein sequences from Arabidopsis thaliana as the reference, which correspond to the gene sequences used for probe design for Hyb-Seq in this study. Paralog detection was performed for all exons with the 'paralog investigator, ' and all assembled loci (with and without paralogs detected) were processed following Sun et al. (2023). After paralog processing, all assembled genes were aligned using MAFFT v.7.409 (Katoh and Standley, 2013) with default parameters, and the poorly aligned sites were trimmed using the plugin trimAl (Capella-Gutiérrez et al., 2009) in Phylosuite (Zhang et al., 2020a) using the automated1 option. We removed genes with alignments shorter than 100 bp following Du et al. (2023). Additionally, we detected and deleted samples with long branches for each gene in ETE3 (Huerta-Cepas et al., 2016) using the pipeline of KewHybSeqWorkshop (Baker et al., 2022; Yang et al., 2023). After filtering, all single gene alignments were concatenated into a super-matrix using the python script "concatenate_fasta.py" (Zhang et al., 2020b; https://github.com/Kinggerm/PersonalUtilities).

2.4. Plastid genome assembly, annotation, and alignment

We used GetOrganelle v.1.6.2d (Jin et al., 2020) for high-quality plastome assembly of the genome skimming data with default settings. Annotation of the plastome was performed with the software PGA (Qu et al., 2019). All annotated data were checked and modified with Geneious v.9.1 (Kearse et al., 2012) to confirm that the start and stop codon of each gene were correctly annotated. In all, 78 coding sequences (CDS) of the plastome were extracted using the python script "get_annotated_regions_from_gb.py" (Zhang et al., 2020b; https://github.com/Kinggerm/PersonalUtilities). To increase taxon sampling in the plastid dataset, we also assembled the same 78 plastid CDS from off-target reads of Hyb-Seq for another 14 Elaeagnaceae species without genome skimming data (Table S4) using the same approach described in Fu et al. (2022). Each CDS was independently aligned and filtered using the same method as applied to the nuclear genes. After filtering, the 78 plastid CDS from the genome skimming data and the off-target data were concatenated for each sample and aligned in a super-matrix.

2.5. Phylogenetic analyses

For Hyb-Seq nuclear data, both concatenated and coalescent approaches were applied to reconstruct the phylogenies of Elaeagnaceae. In the concatenated analysis, we applied ML, implemented in RAxML v.8.2.12 (Stamatakis, 2014) with the GTRGAMMA model and 1000 bootstrap replicates; in the coalescent analysis, ASTRAL v.5.7.8 (Mirarab et al., 2014; Zhang et al., 2018b) was applied. Single-gene ML trees were reconstructed in RAxML v.8.2.12 using the GTRGAMMA model and 100 bootstrap replicates. Following previous studies (e.g., Larson et al., 2020; Morales-Briones et al., 2022), we collapsed nodes with poor support [bootstrap support (BS) < 50%] in individual gene trees prior to using these gene trees as input in ASTRAL for species tree inference. For plastome data, we inferred the ML tree for Elaeagnaceae based on the concatenated alignment of CDS using RAxML under the GTRGAMMA model with 1000 bootstrap replicates. All phylogenetic trees were visualized using FigTree v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/), TVBOT (Xie et al., 2023), and Adobe Illustrator (AI) 2020 software.

2.6. Conflict analysis

We made a one-to-one comparison between the plastid tree and the nuclear coalescent tree, rather than the nuclear concatenated tree due to the higher support values across the coalescent tree. Plastid and nuclear gene matrices were pruned to ensure that the two datasets have the same species. Pruned matrices of plastid CDS and nuclear genes were utilized to perform phylogenetic analysis using the methods above. Both the plastid concatenated tree and the nuclear coalescent tree were rooted with Rhamnus globosa Bunge as the outgroup using the "pxrr" command in Phyx (Brown et al., 2017). After collapsing nodes with BS < 50%, we mapped the nuclear coalescent tree to the plastid concatenated tree using a custom R script (Figshare: https://doi.org/10.6084/m9.figshare.26403217). We also investigated topological concordance and conflict among the nuclear gene trees using PhyParts v.0.0.1 (Smith et al., 2015) and considered BS < 50% to be uninformative following previous studies (e.g., Larson et al., 2020; Morales-Briones et al., 2022). Results were visualized with the script "phypartspiecharts.py" (https://github.com/mossmatters/phyloscripts/tree/master/phypartspiecharts).

2.7. Coalescent simulations and network analysis

We used a coalescent simulation strategy to explore potential causes of the cyto-nuclear discordance detected via conflict analysis, with a special focus on Elaeagnus based on those results. We simulated 1000 plastid trees under the coalescent model using the python package DendroPy v.4.1.0 (Sukumaran and Holder, 2010) with the nuclear ASTRAL species tree as the guide tree with branch lengths rescaled by a factor of four to simulate typical chloroplast inheritance in angiosperms (McCauley, 1994; Fu et al., 2022). Simulated trees were then mapped to the empirical plastid tree by PhyParts and visualized by "phypartspiecharts.py" to calculate expected frequencies for the observed clades in simulated plastid trees. At branches with cyto-nuclear discordance, a substantial number of simulated trees supporting the topology of the empirical plastid tree indicates that ILS can adequately account for the observed conflict. Conversely, if only a few or none of the simulated plastid trees support the empirical plastid topology, it indicates that the observed plastid topology is not expected under ILS alone, and an alternative explanation such as gene flow or chloroplast capture should be invoked to explain the observed cyto-nuclear discordance (e.g., García et al., 2017; Wang et al., 2021a; Fu et al., 2022).

To further explore the possibility of reticulation in the evolutionary history of Elaeagnus as evidenced by the nuclear genome given evidence for high levels of nuclear gene-tree discordance in this genus (see Results), we also inferred phylogenetic networks on the nuclear gene-tree set using PhyloNet v.3.8.0 (Than et al., 2008; Yu et al., 2012; Wen et al., 2018). Due to computational limitations, we reduced the taxon sampling to a computationally tractable size (Than et al., 2008; Wen et al., 2018; Fu et al., 2022; Komarova and Lavrenchenko, 2022); 13 representative species (including one outgroup and 12 species from three Elaeagnus clades) were chosen in the network analyses. We inferred the networks using a maximum pseudo-likelihood approach with the command "InferNetwork_MPL" following Yu et al. (2012). The best network was determined by the Akaike Information Criterion (AIC), the AICc, and the Bayesian Information Criterion (BIC) scores, following the calculation method of Fu et al. (2022). All phylogenetic networks were visualized using Dendroscope (Huson et al., 2012) and the Julia package PhyloPlots (https://github.com/cecileane/PhyloPlots.jl).

2.8. Divergence time estimation

Dating analysis was based on the concatenated alignment of 83 nuclear loci using a penalized likelihood approach in treePL (Smith and O'Meara, 2012). Following the empirical guide described by Maurin (2020), the smoothing parameter was determined using the cross-validation option, and priming was used to determine the best optimization scores. To explore the uncertainty of divergence time estimations, we also conducted the same treePL analysis on 1000 bootstrap ML trees from nuclear data that were generated by fixing the topology of the bestTree in the output of RAxML following Maurin (2020). We used TreeAnnotator v.2.6.7 (Bouckaert et al., 2014) to summarize the dated bootstrapped trees and produce the confidence intervals (CI) for a consensus tree. The configuration files used to run these analyses are provided in the Figshare (https://doi.org/10.6084/m9.figshare.26403217).

Reliable fossils of Elaeagnaceae are few (Su et al., 2014). Multiple fossils have long been considered to belong to Elaeagnaceae (Grigoreva, 1969; Ananiashvili and Purtseladze, 1976; Kovar-Eder, 1984; Song and Qian, 1989) but cannot be placed into an extant genus with confidence. Therefore, we only used accepted fossils that can be placed into a certain clade or were used in previous studies. Apart from two fossils and one secondary calibration in Elaeagnaceae, we additionally used three fossils and four secondary calibrations from other families of Rosales. The first is a fossil of Shepherdia weaveri Becker, discovered in the Mormon Creek region of Montana, the United States, and which has been dated to 33.9 Ma (Becker, 1960). We assigned this fossil to the stem node of Shepherdia (min = 33.9 Ma) following Jia and Bartish (2018). Another is a fossil of Elaeagnus tibetensis T. Su et Z.K. Zhou, discovered in eastern Tibet (Su et al., 2014) and dated to the late Miocene. Su et al. (2014) argued that this fossil most closely resembles E. conferta Roxb., E. difficilis Servettaz, E. luxiensis C.Y. Chang, and E. umbellata Thunb., all of which belong to Clade N3 of Elaeagnus. We therefore assigned this fossil to the crown node of Clade N3 (min = 5.33 Ma, max = 11.63 Ma according to two boundaries of late Miocene) in Elaeagnus. The three fossils used outside Elaeagnaceae include a fossil of Coahuilanthus belindae Calvillo-Canadell & Cevallos-Ferriz, discovered in Coahuila, Mexico in late Campanian that considered to have similar floral morphology with extant Rhamneae (e.g., Rhamnus and Sageretia Brongn.) and Zizypheae (e.g., Berchemia Neck. ex DC.) (Calvillo-Canadell and Cevallos-Ferriz, 2007; Onstein et al., 2015). This fossil was placed at the crown node of the Rhamnaceae (min = 72.1 Ma, max = 83.6 Ma according to two boundaries of the Campanian) following Tian et al. (2024). The fossil of Triorites minutipori Muller was discovered in Malaysia and dated to 89.8 Ma (Muller, 1981). This fossil was assigned to the stem node of Cannabaceae (max = 89.8 Ma) following Fu et al. (2022). Our last non-Elaeagnaceae fossil is of Prunus aspensis Brown, found in the Aspen Shale Formation, Wyoming, the United States, and dated to the Albian, early Cretaceous (Peppe et al., 2008). We assigned this fossil to the stem node of Rosaceae (min = 100.5 Ma, max = 113.0 Ma) following Xiang et al. (2017). In addition, a few nodes were calibrated using previously estimated dates. The crown age of Hippophae was calibrated at 18.3–24.0 Ma [95% highest posterior density (HPD)] as estimated by Jia and Bartish (2018); the crown age of Moraceae was calibrated at 73.2–84.7 Ma [95% HPD] as estimated by Zhang et al. (2019); the crown age of Ulmaceae was calibrated at 75.99–95.21 Ma [95% HPD] as estimated by Zhang et al. (2022); the crown age of Rosaceae was calibrated at 94.5–96.4 Ma [95% HPD] as estimated by Zhang et al. (2017); and the stem age of Rosales was calibrated at 104–115 Ma [95% HPD] as estimated by Wang et al. (2009).

2.9. Biogeographic analysis

The geographic distributions of sampled species were obtained from multiple online open sources, including the Flora of China (Qin et al., 2007), Global Biodiversity Information Facility (https://www.gbif.org/, 2023-12-12), and POWO. Given species endemism of some Elaeagnaceae and following the floristic delineations proposed by Wu and Wu (1996), we delimited eight biogeographic regions for analysis: (A) Holarctic region in Eurasia; (B) Qinghai-Tibet Plateau (QTP); (C) Sino-Himalayan subregion; (D) mainland East Asia; (E) Taiwan (China); (F) Japan; (G) Indochina-Malaysia-Australia; and (H) North America.

The ancestral area reconstruction of each node in the phylogeny was inferred using BioGeoBEARS (Matzke, 2018), implemented in the GUI software RASP (Yu et al., 2020). Six biogeographic models (i.e., DEC, DEC + J, DIVALIKE, DIVALIKE + j, BAYAREALIKE, and BAYAREALIKE + j) provided by BioGeoBEARS were tested to select the best-fit model using the log-likelihood (LnL), AIC, and AIC_wt values. The number of maximum areas for ancestral nodes was set at five, which represents the maximum number of areas observed in extant Elaeagnaceae species. We carried out the biogeographical analyses based on the dated ML tree from nuclear data with the outgroup pruned.

3. Results 3.1. Assembly and alignments

The number of assembled nuclear genes for each sample ranged from 32 in Elaeagnus cinnamomifolia W.K. Hu & H.F. Chow to 79 in Trema amboinensis (Willd.) Blume, with an average of 66 nuclear genes recovered (Table S5). The number of exons with paralog warnings ranged from zero in Hippophae rhamnoides and E. montana Makino to 20 in E. bockii Diels (Table S6). After cleaning, the dataset used in phylogenetic analysis retained 83 genes for 76 samples (including 64 samples from three genera of Elaeagnaceae and 12 outgroup samples of Rhamnaceae) with a minimum gene size of 537 bp and a maximum size of 4044 bp. These genes were assembled into a concatenated matrix of aligned length 105, 876 bp with 12, 696 singleton sites, and 74, 042 constant sites. The number of assembled CDS from off-target data for each of the 14 samples for which this was possible ranged from four in E. montana to 66 in E. luoxiangensis C.Y. Chang, with an average of 27 CDS being recovered (Table S4). Plastomes of all 63 samples from genome skimming data were successfully assembled into circular genomes with sequence length ranging from 150, 390 bp (E. angustifolia L.) to 156, 512 bp (H. rhamnoides), and 78 CDS were successfully extracted from each sample. After data filtering and combining, the concatenated plastid matrix, including 71 samples from three genera of Elaeagnaceae and two outgroup samples, had an aligned length of 52, 488 bp with 1364 singleton sites and 46, 763 constant sites. More detailed information on the assembled results is available in Appendix A.

3.2. Phylogenetic relationships resolved by nuclear data

The results from both concatenated and coalescent phylogenetic analyses fully supported (BS = 100%, LPP = 1.00) the monophyly of Elaeagnaceae and each of the three genera (Fig. 1). Both analyses fully supported Elaeagnus as sister to a clade comprising Shepherdia and Hippophae (BS = 100%, LPP = 1.00).

For intrageneric relationships, both analyses resolved Elaeagnus into three moderately to strongly supported (BS > 70%, LPP > 0.7) clades (Clades N1–N3 in Figs. 1 and S1), with Clades N1 and N2 being successive sisters to Clade N3. Within Clade N1 of Elaeagnus, E. mollis Diels was resolved as the sister of a strongly (BS = 100%, LPP = 0.92) supported clade comprising E. commutata Bernh. ex Rydb. and E. angustifolia; Clade N2 of Elaeagnus includes eight species (E. angustata (Rehder) C.Y. Chang, E. wushanensis C.Y. Chang, E. magna (Servett.) Rehder, E. guizhouensis C.Y. Chang, E. stellipila Rehder, E. bambusetorum Hand.-Mazz., E. jingdonensis C.Y. Chang, and E. calcarea Z.R. Xu) in both concatenated and coalescent analyses, but the two analyses show different topologies; Clade N3 of Elaeagnus includes more than 70% of the species of the genus in both the concatenated and coalescent trees, but the support values of many internal nodes are relatively poor (BS < 50% or LPP < 0.50) (Figs. 1 and S1). Hippophae was resolved into two clades by concatenated and coalescent analyses but showed different interspecific relationships (Figs. 1 and S1). The concatenated analysis placed H. gyantensis (Rousi) Y.S. Lian and H. tibetana Schltdl. in a clade with strong support (BS = 99%) and placed H. rhamnoides, H. salicifolia D. Don, and H. neurocarpa S.W. Liu & T.N. He into another clade with low support (BS = 50%). However, the coalescent analysis placed H. gyantensis and H. neurocarpa in a clade with strong support (LPP = 1.00) and the remaining three species into another clade with strong support (LPP = 0.97).

3.3. Phylogenetic relationships resolved by plastid data

Consistent with nuclear phylogenetic analyses, the plastid phylogenomic analysis also fully supported the monophyly of Elaeagnaceae and each of the three genera, with Elaeagnus being sister to a fully supported (BS = 100%) clade comprising Shepherdia and Hippophae. Elaeagnus was resolved into three well-supported (BS > 90%) clades (Clades P1–P3 in Fig. S2), with Clades P1 and P2 being strongly supported (BS = 100%) as successive sisters to Clade P3. Within Clade P1 of Elaeagnus, E. mollis was resolved as the sister of a strongly (BS = 100%) supported clade comprising E. oxycarpa Schltdl. and E. angustifolia; in Clade P2, plastid data well resolved and supported (BS > 95%) interspecific relationships with E. angustata as sister to the remaining species, which further resolved into a clade comprising E. bambusetorum, E. calcarea, and E. jingdonensis and a clade of E. nanchuanensis C.Y. Chang, E. wushanensis, and E. multiflora var. tenuipes C.Y. Chang; Clade P3 also includes more than 70% of the species of the genus, but the plastid tree is better resolved and supported than the nuclear tree, with ca. 50% of the internal nodes of the plastid tree receiving strong support (BS > 90%). Hippophae was resolved into two well-supported (BS = 100%) clades: one comprising H. rhamnoides and H. tibetana and the other H. neurocarpa, H. salicifolia, and H. gyantensis. In Shepherdia, S. rotundifolia Parry was resolved as the sister to a fully supported (BS = 100%) clade of S. argentea (Pursh) Nutt. and S. canadensis (L.) Nutt.

3.4. Conflict analysis

The conflict analysis based on the nuclear gene trees, as assessed by PhyParts, showed that most gene trees (68/80 = supporting nuclear gene trees/informative nuclear gene trees [i.e., BS > 50%]) supported the monophyly of Elaeagnaceae (Fig. 2). The monophyly of Elaeagnus, Hippophae, and Shepherdia was each supported by 39/73, 54/69, and 44/64 nuclear gene trees, respectively. Within Elaeagnus, Clades N1, N2, and N3 were supported by 28/42, 13/43, and 3/50 nuclear gene trees, respectively. Within Clade N3 of the nuclear coalescent tree, consistent nuclear gene trees were observed infrequently; instead, gene tree conflicts were commonly observed at many nodes, leading to low support (LPP < 0.5). Additionally, we also mapped the nuclear coalescent tree to the plastid concatenated tree (Fig. 3), and cyto-nuclear discordances were observed in Hippophae and Clades N2/P2 and N3/P3 of Elaeagnus. In Hippophae, the nuclear coalescent tree supported H. neurocarpa and H. gyantensis forming a clade, and the remaining three species formed another clade, while the plastid tree supported H. neurocarpa, H. gyantensis, and H. salicifolia forming a clade and the remaining two species formed another clade. In Elaeagnus, although Clades N2 have almost the same species composition as P2, Clades N3 have almost the same species composition as P3, significant cytonuclear conflicts were detected within each clade and two species were even found to have conflicting positions between two clades (Fig. 3).

3.5. Coalescent simulations and network analysis

We used coalescent simulations and network analysis to determine whether conflicts between nuclear and plastid trees could be explained by ILS (Fig. S3). Overall, the phylogenetic backbone was found to be highly probable given the species tree and therefore within ILS expectations. We found little gene tree heterogeneity in the Elaeagnaceae plastid backbone. Specifically, all simulations supported Hippophae, 94.6% supported Shepherdia, and 99.2% supported Elaeagnus. Likewise, within Elaeagnus, 92.5% of simulations supported Clade P1, and 99.7% supported the clade comprising Clades P2 and P3 in the plastid tree. However, at more shallow levels within Elaeagnus, simulations showed far less support for the plastid clades, including clades that were well-resolved with nuclear data. For example, no simulations recovered Clades P2 and P3; similarly, no simulations recovered the internal nodes of Elaeagnus except for the subclade comprising E. latifolia L. and E. conferta, which was supported by 73.9% of the simulated trees. These findings indicate that the cyto-nuclear discordance observed at these nodes cannot be explained by ILS alone, and an alternative explanation (e.g., gene flow or chloroplast capture) should be invoked. In Hippophae, all clades similarly had low probability of ILS, although these were not close to zero and therefore less decisive; the sister relationship between H. rhamnoides and H. tibetana in the plastid tree had 14.8% probability, and the clade composed of H. gyantensis, H. salicifolia, and H. neurocarpa had 4.2% probability.

The PhyloNet analysis based on nuclear sequence data inferred a phylogenetic network with five reticulations as the best-fit model with the lowest AIC, AICc, and BIC scores (2170.6, 2063.7, and 2189.2, respectively) (Fig. S4 and Table S7). This best network suggested gene flow between Elaeagnus angustifolia (belonging to Clade N1) and an unsampled or extinct lineage ancestral to E. courtoisii Belval (belonging to Clade N3). An unsampled or extinct lineage was also inferred to have contributed to a reticulation event between E. bambusetorum (belonging to Clade N2) and E. tarokoensis S.Y. Lu & Yuen P. Yang (belonging to Clade N3). Within Clade N3, E. courtoisii was inferred to have parental inheritance from an unsampled or ancestral lineage of the same species; the clade comprising E. umbellata and E. multiflora Thunb. was inferred to have resulted from a reticulation event between two unsampled lineages; and the clade of E. formosana Nakai and E. sarmentosa Rehd. was inferred to have parental inheritance from E. formosensis Hatusima and an unsampled lineage ancestral to Clades N2 and N3. It is worth noting that all five detected reticulations in Elaeagnus involved or occurred within Clade N3, a finding concordant with high nuclear gene tree heterogeneity (Fig. 2) and strong cytonuclear discord (Fig. 3) largely localized to clade N3.

3.6. Divergence time estimation

Based on our nuclear dataset, the optimal smoothing value for the treePL analysis (as indicated by the cross-validation tests) was 1, and the stem and crown ages of Elaeagnaceae were estimated to be in the late Cretaceous (89.41 Ma; 95% CI = 88.41–90.45 Ma) and the late Eocene (40.78 Ma; 95% CI = 40.32–41.30 Ma), respectively (Fig. 4 and Table S8). The stem ages of Shepherdia and Hippophae were estimated as 33.9 Ma (95% CI = 33.90–33.91 Ma), the crown age of Shepherdia as 25.10 Ma (95% CI = 24.08–26.26 Ma), the crown age of Hippophae as 18.3 Ma (95% CI = 18.30–18.31 Ma), the stem age of Elaeagnus as 40.78 Ma (95% CI = 40.32–41.30 Ma), and the crown age of Elaeagnus as 14.72 Ma (95% CI = 14.03–15.33 Ma). Within Elaeagnus, Clade N1 diverged from the remaining species of this genus at 10.51 Ma (95% CI = 9.96–11.05 Ma), and Clades N2 and N3 diverged at 6.63 Ma (95% CI = 6.24–7.06 Ma).

3.7. Ancestral area reconstruction

Based on our nuclear dataset, the comparison among AICc_wt scores of different biogeographical models favored the BAYAREALIKE + j model as the best-fit model for Elaeagnaceae (AIC_wt = 0.72; Table S9). The BAYAREALIKE + j model indicates that the most recent common ancestral area (crown node) of Elaeagnaceae is the QTP (Figs. 5 and S5; Table S10). The common ancestor of Hippophae + Shepherdia was distributed in the QTP; the ancestor of Hippophae was found in the QTP; the most recent common ancestor of extant Shepherdia was distributed in the North America; the ancestor of Elaeagnus was distributed in mainland East Asia (Figs. 5 and S5; Table S10).

A series of dispersals was inferred to have shaped current distributional patterns across Eurasia to North America in the northern temperate zone: the common ancestor of Shepherdia and Hippophae dispersed from the QTP toward North America, and Hippophae dispersed from the QTP toward the Holarctic region in Eurasia as well as from the QTP toward mainland East Asia. Our results also support one dispersal in Clade N1 of Elaeagnus from Eurasia toward North America. Finally, in Clade N3 of Elaeagnus we detected five dispersals from the Sino-Himalayan subregion and/or mainland East Asia toward Indochina-Malaysia-Australia, four dispersals from the Sino-Himalayan subregion and/or mainland East Asia toward Taiwan (China), and five dispersals from the Sino-Himalayan subregion and/or mainland East Asia toward Japan (Fig. 5).

4. Discussion 4.1. Intergeneric relationships in Elaeagnaceae

Both nuclear and plastid datasets fully supported the monophyly of Elaeagnaceae and each of its three genera (Figs. 1, S1 and S2), which agrees with previous molecular results (Soltis et al., 2000, 2011; Sytsma et al., 2002; Zhang et al., 2011). Our results indicate that Elaeagnus is sister to a clade comprising Shepherdia and Hippopahe. The three genera can be easily distinguished according to their morphological traits (Bartish and Swenson, 2004). Elaeagnus has alternate leaves, and the plants are monoecious, whereas Hippophae and Shepherdia have opposite leaves, and the plants are dioecious. The shared features between Hippophae and Shepherdia are consistent with their closer phylogenetic relationship. Hippophae and Shepherdia can be distinguished by the shape and size of their leaves (linear or linear lanceolate vs. larger lanceolate or oval leaves) and the number of lobes on their calyxes (2-lobed vs. 4-lobed) (Zhang, 2016).

4.2. Intrageneric relationships in Hippophae

Intrageneric relationships within Hippophae have been inconsistent across previous phylogenetic studies: Bartish et al. (2002) resolved H. rhamnoides as sister to the clade comprising the remaining four species using plastid RFLPs and morphological characters; Sun et al. (2002) supported H. tibetana, H. rhamnoides as successive sisters to the remaining species using the ITS data; Jia and Bartish (2018) resolved Hippophae into two clades, one comprising H. tibetana and H. rhamnoides, and the other clade comprising the remaining species using both cpDNA and nDNA. Our coalescent analysis based on nuclear data placed H. gyantensis and H. neurocarpa in one clade, and H. rhamnoides, H. salicifolia, and H. tibetana in another clade, whereas our concatenated analysis based on nuclear data placed H. gyantensis and H. tibetana in one clade and the remaining three species into another clade (Figs. 1 and S1). Both topologies differ from upper mentioned three previous phylogenetic results. However, our concatenated tree based on plastid data (Fig. S2) recovered a clade of H. rhamnoides and H. tibetana, while the remaining Hippophae species formed another clade, which is congruent with findings based on five nDNA fragments (Jia and Bartish, 2018). We only sampled five species and limited intraspecific taxa; thus, further phylogenetic analysis using a large number of nuclear genes and samples is needed to clarify the intrageneric relationships of this genus.

4.3. Intrageneric relationships in Elaeagnus

Our phylogenomic analyses based on nuclear and plastid DNA sequence data with broad taxonomic sampling provide new insights into intrageneric relationships of Elaeagnus, the species of which form three strongly supported clades. Clades N1/P1 largely correspond to group IIA of sect. Deciduae proposed by Sun and Lin (2010) based on leaf duration and fruit texture. This clade is also congruent with Clade I resolved by Zhang (2016). Clades N2/P2 in this study largely correspond to Clade II resolved by Zhang (2016), containing both semi-evergreen taxa (e.g., E. stellipila, E. guizhouensis, E. calcarea, and E. jingdonensis) and deciduous taxa (e.g., E. angustata and E. wushanensis). Clades N3/P3 in this study are composed of Clades III, IV, and V of Zhang (2016) and contain evergreen Elaeagnus species (e.g., E. conferta, E. triflora Roxburgh), some semi-evergreen species (e.g., E. argyi H.Léveillé), and also deciduous species (e.g., E. courtoisii, E. umbellata, and E. multiflora). The nuclear sequence data poorly resolved most interspecific relationships, whereas the plastid data provided better resolution and support with ca. 50% of the internal nodes receiving strong support (BS > 90%). Clades N3 and P3 show significant incongruence; however, most relationships in the nuclear trees are poorly supported, as described above.

4.4. Hybridization in Elaeagnaceae

Our coalescent simulation analysis showed zero or nearly-zero simulated trees support the plastid topology (Fig. S3), indicating that ILS alone cannot explain the cyto-nuclear discordance observed in Elaeagnus. Instead, the observed conflicts are likely explained by gene flow (see also Wu et al., 2015; Feng et al., 2019a; Stubbs et al., 2020; Xiao et al., 2022). This explanation is further supported by PhyloNet analysis, which inferred five reticulation events, with all detected reticulations in Elaeagnus occurring within or related to Clade N3. However, few studies have explored hybridizations in Elaeagnus. Servettaz (1909) reported a hybridization between E. pungens Thunb. and E. macrophylla Thunb. (both belong to Clade N3 in this study). Both nuclear and plastid markers have indicated that hybridization gave rise to E. × submacrophylla Servett. (Zhang, 2016). Morphological characters as well as nuclear and plastid markers were used to confirm that a hybridization between E. glabra Thunb. (belongs to Clade N3 in this study) and E. macrophylla gave rise to E. × maritima Koidz. (Ohba, 1999; Jang et al., 2023). Thus, the most plausible cause of cyto-nuclear discordance and conflicts among nuclear gene trees in Clade N3 of Elaeagnus is hybridization.

The Hippophae phylogenetic trees based on nuclear genes using both concatenate and coalescent method conflicted with those based on plastid data. At all internal nodes, coalescent simulation analysis indicated that ILS was an unlikely explanation for these cyto-nuclear conflicts. These conflicts may instead by explained by the occurrence of frequent hybridizations within this genus. The overlapping geographical distributions and largely similar flowering and fruiting phenologies of the five sampled Hippophae species suggest that species within the genus have adequate opportunity for interspecific hybridization (Lian et al., 2003; Bartish and Thakur, 2022). Furthermore, previous studies have identified hybridizations within Hippophae. For example, Bartish et al. (2002) observed a hybridization between H. rhamnoides ssp. sinensis Rousi and H. neurocarpa ssp. neurocarpa, and the hybridization gave rise to H. goniocarpa Y.S. Lian & al. ex Swenson & Bartish based on RAPD markers, which was supported by Sun et al. (2003) and Wang et al. (2008). More recently, Jia et al. (2016) proposed that H. gyantsensis should originated from multiple crosses between H. rhamnoides ssp. yunnanensis Rousi and H. neurocarpa. However, our study could not accurately infer these hybridization events due to insufficient specific and infraspecific sampling.

4.5. Biogeographical history 4.5.1. Origin and early diversification of Elaeagnaceae

Elaeagnaceae are inferred to have originated during the late Cretaceous (the Turonian period, 89.41 Ma; 95% CI = 88.41–90.45 Ma). The three extant genera are estimated to have diverged from each other during the middle Eocene (40.78 Ma; 95% CI = 40.32–41.30 Ma), which is largely congruent with previous studies (Li et al., 2015; Jia and Bartish, 2018; Sun et al., 2020). However, our dating is much older than that proposed by Li et al. (2019), which placed the fossil Coahuilanthus belindae at stem instead of crown of Rhamnaceae. Additionally, we noticed there is a ~50-million-year gap between the origin of the family and the divergence of its current lineages (Fig. 4). Multiple fossils of Elaeagnaceae have been reported in the Northern Hemisphere during that period. For instance, the reported oldest Turonian fossil was excavated in Siberia (Grigoreva, 1969); a Maastrichtian fossil was found in the Taizhou Formation in North Jiangsu Basin, East China (Song and Qian, 1989); multiple late Cretaceous fossils were detected in Georgia and Hungary (Ananiashvili and Purtseladze, 1976; Kovar-Eder, 1984); and Paleogene fossils were found in Germany (Weber, 1852). Therefore, we infer that the long gap between the stem and crown divergence in Elaeagnaceae might be due to extinctions. Other taxa such as Cycas L. (Mankga et al., 2020) have a similarly long gap between origin and diversification of extant lineages, which is also thought to be caused by considerable extinction events based on both diverse Mesozoic fossils and simulation analyses (Crisp and Cook, 2009; Nagalingum et al., 2011). The discovery of the oldest Elaeagnaceae fossil in Siberia (Grigoreva, 1969) indicates that the family could have originated in mid-to-high latitude regions of the Northern Hemisphere.

Divergence between Elaeagnus and the clade comprising Shepherdia and Hippophae was inferred to have occurred during the middle Eocene in the QTP (Fig. 4, Fig. 5). Divergence between Shepherdia and Hippophae have occurred in the QTP at the Eocene-Oligocene boundary when the global climate began to change significantly with a sharp drop in temperature (Cavelier et al., 1981; Ivany et al., 2000; Zachos et al., 2001; Wang et al., 2021b). The climatic deterioration of the Eocene and Miocene forced many lineages that originated in the north temperate zone to adapt to drier and colder climates or withdraw toward the equator for survival (Morley, 2011). Elaeagnaceae could have survived this climatic transition through both processes, as is evidenced by the extremely strong cold and drought resistance of E. angustifolia and adaptation to tropical dense forest habitats and the emergence of evergreen traits in E. conferta, E. latifolia, and E. triflora.

4.5.2. Dispersal and the formation of the current distribution of the three genera of Elaeagnaceae

The common ancestor of Shepherdia and Hippophae was estimated to have originated in the QTP during the early Oligocene (Fig. 4, Fig. 5, Tables S8 and S11). The disjunction between Hippophae and Shepherdia between Eurasia and North America indicates they have experienced allopatric diversification (Bartish and Thakur, 2022). All excavated Shepherdia fossils are in North America without exception, e.g., the Oligocene fossil in Mormon Creek, the United States (Becker, 1960), the early Miocene fossil in Alaska, the United States (Wolfe, 1969), and the middle Miocene fossil in Yukon, Canada (White et al., 1997). As is the case for many temperate angiosperm taxa (Wen, 1999, 2001; Valcárcel et al., 2019), the common ancestor of Shepherdia and Hippophae may have migrated from East Asia into North America via the Bering Land Bridge, which was open to terrestrial organisms from at least the early Paleocene until its closure between 7.4 and 4.8 Ma (Tiffney and Manchester, 2001; Yi et al., 2004). This ancestor would then have evolved into extant Shepherdia species in North America.

Our analyses indicate that extant Hippophae lineages originated in the QTP and diverged during the middle Miocene (Fig. 5, Tables S4 and S5), which is largely congruent with earlier studies (Jia et al., 2012; Jia and Bartish, 2018). Most Hippophae species are distributed in the QTP and its adjacent areas, as is the case for many other angiosperm taxa (Ren et al., 2018; Feng et al., 2019b; Xu et al., 2019) that originated in the QTP and experienced rapid speciation during the middle Miocene. The timing of Hippophae diversification also provides an explanation for the distribution of wide-ranging species. For example, we infer that Hippophae rhamnoides, a species with a broad distribution and multiple subspecies across Eurasia, originated in the QTP. This finding is consistent with a previous study based on more comprehensive sampling of subspecies in H. rhamnoides that suggested this species started to diverge in the QTP during the middle Miocene and dispersed out of the QTP in the late Miocene (Jia and Bartish, 2018). The dispersal into Europe and mainland East Asia of multiple plants groups that originated in the QTP and previously thrived in alpine environments (e.g., Rhodiola L., see Zhang et al., 2014; Physochlaina G. Don, see Lei et al., 2021; Actaea L., see Ling et al., 2023; Pleione D. Don, see Wu et al., 2023; and Lilium Tourn. ex L., see Zhou et al., 2024) coincides with a substantial decrease in global temperatures that began in the middle Miocene (Zachos et al., 2001).

Ancestral area reconstruction and dating analysis suggest that extant Elaeagnus lineages started to diverge in mainland East Asia (area D) during the middle Miocene (14.72 Ma, 95% CI = 14.03–15.33 Ma; Fig. 4, Fig. 5, Tables S8 and S10). During this time, the Himalayas reached an elevation of 2.3 km before rapidly uplifting to a higher elevation (ca. 5.5 km) in the middle Miocene (ca. 15 Ma), which may have intensified the East Asia monsoon (Mao et al., 2021) and in turn led to the diversification of some plant genera, e.g., Lepisorus (J. Sm.) Ching (Wang et al., 2012) and Fagus L. (Wan et al., 2023). These significant geological and climatic changes might also have facilitated the diversification of Elaeagnus.

Elaeagnus was inferred to have experienced multiple dispersals toward other regions, including one dispersal toward North America after the middle Miocene, five toward Indochina and Australia after the Pliocene, four toward Taiwan (China) and five toward Japan after the Pliocene (Fig. 5). Other plant taxa that originated on the QTP or Sino-Himalayan subregion also experienced multiple dispersals to regions in the Northern Hemisphere during the same period, e.g., Saxifraga Tourn. ex L. (Ebersbach et al., 2017) and Polygonatum Mill. (Xia et al., 2022).

Our analyses indicate that Elaeagnus species from Clade N1 dispersed from East Asia into North America (Fig. 5) and that this dispersal could not have occurred earlier than 10.51 Ma (95% CI = 9.96–11.05 Ma). Fossil evidence indicates that at least one Elaeagnus species belonging to Clade N1 (i.e., E. angustifolia, excavated in Alaska) migrated to North America via the Bering Land Bridge during this period (Leopold and Liu, 1994). As the Bering Land Bridge is thought to have acted as an important corridor for multiple temperate plant taxa during the Miocene (Wen, 1999; Tiffney and Manchester, 2001; Yi et al., 2004), it is highly likely that it facilitated the dispersal of species in Clade N1 of Elaeagnus. Alternatively, long-distance dispersal in Elaeagnus may have been facilitated by birds (Choi and Chae, 2007; Edwards et al., 2014).

Additional dispersal events were identified in a subclade of the Elaeagnus Clade N3. This subclade, which consists of Elaeagnus conferta, E. triflora, and E. latifolia, is estimated to have begun its crown divergence 5.50 Ma (95% CI = 5.15–5.95 Ma). The subclade, which was inferred to have originated in East Asia during the late Miocene (Fig. 5, Tables S8 and S10), shows a broad distribution across mainland East Asia, the Indochina Peninsula and Indonesian islands, and Australia (areas referred to in our biogeographic analysis as Indochina-Malaysia-Australia). The current distribution of Elaeagnus species in this subclade was likely facilitated by long-distance dispersals. Studies have suggested that the intermittent formation of land bridges between islands in Southeast Asia from approximately 2.7 Ma (Woodruff, 2010) would have facilitated the southward dispersal of multiple plant taxa, e.g., Ilex Tourn. ex L. (Yao et al., 2021), Oreocnide Miq. (Wu et al., 2022b) and Cymbidium Sw. (Chen et al., 2024). Our findings suggest that the southeastward dispersals of three Elaeagnus species may have occurred via these land bridges. The remaining four dispersals in Clade N3 may also have been facilitated by land bridges, although birds might also have played an important role in their dispersals.

Elaeagnus was inferred to have migrated into Japan five times after the Pliocene (Fig. 5, Tables S8 and S10), most likely by long-distance dispersal. Japan has been separated from mainland East Asia by the Japanese Sea since the Miocene (25–15 Ma; Maruyama et al., 1997), and this timeframe is much earlier than the inferred dispersal events. Therefore, we infer that the five dispersals from mainland East Asia toward Japan might have been facilitated by frugivorous birds rather than land bridges.

Four dispersals were detected from mainland East Asia to Taiwan after the Pliocene. A land bridge formerly known as "Dongshan" (Lin, 1982; Zeng, 1993) emerged many times due to sea level regression during the Pleistocene. The Dongshan land bridge may have facilitated floristic exchanges between mainland China and Taiwan (Zeng, 1993), including the dispersals toward Taiwan during that period in Elaeagnus. However, we cannot disregard the possibility of these dispersals occurring through transport by birds.

5. Conclusion

This study presents the most comprehensive phylogenetic framework for Elaeagnaceae to date and supports Elaeagnus as sister to the clade of Hippophae and Shepherdia. This study successfully resolves intrageneric relationships, especially for Elaeagnus, which is resolved into three strongly supported clades. High levels of cyto-nuclear discordance as well as nuclear gene-tree conflicts in Clade N2/P2 and Clade N3/P3 of Elaeagnus were observed. Our coalescent simulations suggested that ILS alone cannot explain the high level of cyto-nuclear discordance observed in Elaeagnus. PhyloNet results indicated that these observed conflicts could be explained by hybridization events. Fossil evidence indicates that Elaeagnaceae may have originated in mid-to-high latitude regions of the Northern Hemisphere during the late Cretaceous. The observed 50-million-year gap between the origin of Elaeagnaceae and the divergence of its current lineages may have been caused by a series of extinctions. Our findings also suggest that rapid uplift of the Himalayas and the intensified formation of the East Asian monsoon likely facilitated diversification of Elaeagnaceae. In addition, estimated divergence times of Elaeagnaceae suggest that long-distance dispersal has played an important role in the formation of current intercontinental distributions, and both bird-mediated dispersal and land bridges likely facilitated these patterns.

Acknowledgements

This research was supported by the National Natural Science Foundation of China, Key International (regional) Cooperative Research Project (no. 31720103903), the Science and Technology Basic Resources Investigation Program of China (no. 2019FY100900), the Strategic Priority Research Program of Chinese Academy of Sciences (no. XDB31000000), the National Natural Science Foundation of China (no. 31270274), the Yunling International High-end Experts Program of Yunnan Province, China (no. YNQR-GDWG-2017-002 and no. YNQR-GDWG-2018-012), the China Scholarship Council (202004910775), the Chinese Academy of Sciences (CAS) President's International Fellowship Initiative (no. 2020PB0009), the China Postdoctoral Science Foundation (CPSF), and the United States Department of Energy (grant no. DE-SC0018247 to PSS, RPG, and DES). We are grateful to the following institutions for providing specimens or silica-dried materials: the herbarium of the California Academy of Sciences (CAS); the herbarium of Kunming Institute of Botany, Chinese Academy of Sciences (KUN); the Germplasm Bank of Wild Species and Molecular Biology Experiment Center, Kunming Institute of Botany, Chinese Academy of Sciences; the Missouri Botanical Garden Herbarium (MO); the New York Botanical Garden Herbarium (NY); the Ohio State University Herbarium (OS) and the University of Texas Herbarium (TEX). We are also grateful to Jiajin Wu for help with sampling and DNA extraction; to Germplasm Bank of Wild Species (Kunming institute of botany, Chinese Academy of Science) for providing genome skimming data of Elaeagnaceae; to Tiantian Xue (Institute of Botany, Chinese Academy of Sciences) for his suggestions on the project; and to the iFlora High Performance Computing Center of Germplasm Bank of Wild Species (iFlora HPC Center of GBOWS, KIB, CAS) for computing.

Data availability

The target enrichment bait set, nuclear and plastid DNA alignments, phylogenetic trees, and other biogeographic data and results from this study are provided in the Figshare data repository https://doi.org/10.6084/m9.figshare.26403217.

CRediT authorship contribution statement

Wei Gu: Visualization, Software, Methodology, Investigation, Formal analysis. Ting Zhang: Resources, Methodology, Investigation, Data curation. Shui-Yin Liu: Writing – review & editing, Software. Qin Tian: Writing – review & editing, Resources, Investigation, Data curation. Chen-Xuan Yang: Writing – review & editing, Software, Methodology, Formal analysis. Qing Lu: Writing – review & editing, Software. Xiao-Gang Fu: Writing – review & editing, Software, Formal analysis. Heather R. Kates: Resources, Investigation. Gregory W. Stull: Software, Resources, Data curation. Pamela S. Soltis: Writing – review & editing. Douglas E. Soltis: Writing – review & editing. Ryan A. Folk: Writing – review & editing. Robert P. Guralnick: Writing – review & editing. De-Zhu Li: Resources, Methodology, Investigation, Funding acquisition, Data curation, Conceptualization. Ting-Shuang Yi: Writing – review & editing, Supervision, Resources, Methodology, Investigation, Funding acquisition, Data curation, Conceptualization.

Declaration of competing interest

No conflict of interest exits with the submission of this manuscript, and it has been approved by all authors for publication.

Appendix A. Supplementary data

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