1. Introduction

The mulberry family (Moraceae) comprises approximately 1140 species in 48 genera (https://powo.science.kew.org/, last accessed February 2023) and is predominantly distributed in tropical and subtropical regions, with some species in a few genera extending into temperate areas (Zhou and Gilbert, 2003). Moraceae species exhibit remarkable morphological diversity in growth form (trees, shrubs, lianas, to herbs), typically contain milky or watery latex, and have inconspicuous flowers aggregated into distinctive inflorescences, including racemes, cymes, spikes, heads, and syconia (figs) (Zhou and Gilbert, 2003). As a model for breeding system evolution, Moraceae display unparalleled diversity in sex systems, encompassing monoecy, dioecy, gynodioecy, and andromonoecy (e.g., Berg, 2001; Zhou and Gilbert, 2003). Some species are economically important for sericulture (e.g., Maclura spp. and Morus spp.), papermaking (e.g., Broussonetia spp.), and as sources of edible fleshy fruits (e.g., Artocarpus spp., Ficus spp., and Morus spp.) consumed by animals and humans. Moraceae are also ecologically important, with many Ficus species forming dominant components of tropical rainforests (Shanahan et al., 2001; Bain et al., 2014).

Recent phylogenetic analyses using a large number of nuclear gene loci (e.g., Zerega and Gardner, 2019; Gardner et al., 2021a, 2021b, 2023) have consistently supported the division of Moraceae into seven tribes (i.e., Antiarideae, Artocarpeae, Chlorophoreae, Dorstenieae, Ficeae, Moreae, and Parartocarpeae) and the monophyly of most genera (e.g., Datwyler and Weiblen, 2004; Clement and Weiblen, 2009; Zerega et al., 2010; Chung et al., 2017), except for Batocarpus, Clarisia, and Dorstenia (Zhang et al., 2019a; Gardner et al., 2021c). However, the monophyly and circumscription of some infrageneric taxa, such as two of the six subgenera within Ficus (i.e., subg. Ficus and subg. Pharmacosycea), require further investigation (Rasplus et al., 2021; Gardner et al., 2023). While tribe-level relationships are largely resolved (Zerega and Gardner, 2019; Gardner et al., 2021a), the position of Chlorophoreae remains uncertain, recovered as either sister to the Moreae + Artocarpeae clade (e.g., Chung et al., 2017; Zhang et al., 2019b), or sister to a clade comprising Parartocarpeae, Dorstenieae, Antiarideae, and Ficeae (hereafter, the PDAF clade) (Zerega et al., 2010; Zerega and Gardner, 2019; Gardner et al., 2021a). The placement of Bagassa is also contentious (Chung et al., 2017; Williams et al., 2017).

Beyond these nuclear–nuclear phylogenetic incongruences discussed above, significant cytonuclear discordances have been documented in major clades of Moraceae, particularly in Parartocarpeae (Gardner et al., 2021a) and Ficus (Bruun-Lund et al., 2017; Gardner et al., 2023). Parartocarpeae was found to be non-monophyletic in the plastid tree but was strongly supported as monophyletic in nuclear trees (Gardner et al., 2021a, 2023). Similarly, in Ficus, four of the six subgenera were monophyletic in nuclear trees (Rasplus et al., 2021; Gardner et al., 2023), whereas none were monophyletic in plastid trees (Bruun-Lund et al., 2017; Gardner et al., 2023).

Nuclear–nuclear and cytonuclear incongruences may arise from multiple biological processes, including incomplete lineage sorting (ILS) and gene flow (Folk et al., 2017; Cai et al., 2021; Kipkoech et al., 2025). ILS occurs when ancestral polymorphisms are randomly inherited by descendant species, causing some genes to reflect evolutionary histories that differ from the species tree—a pattern particularly common during rapid speciation or recent divergence (Scherz et al., 2022). Introgression, defined as the exchange of genetic material between species or populations after divergence, can blur genetic boundaries and lead to conflicting phylogenetic signals (Liu et al., 2025). Recent phylogenomic studies suggest that hybridization during the early evolution of Ficus was a major contributor to its exceptional cytonuclear discordance (Wang et al., 2021; Gardner et al., 2023). Untangling the sources of discordance is essential not only for clarifying phylogenetic relationships, but also for providing insights into historical evolutionary processes. While introgression patterns have been investigated in specific moraceous clades (e.g., Wang et al., 2021; Gardner et al., 2023), a comprehensive family-wide study of deep nuclear–nuclear and cytonuclear discordance is lacking.

In this study, we conducted a comprehensive phylogenomic study of Moraceae, integrating 93 nuclear loci from 326 samples (319 moraceous species, 45 genera) and 80 plastid genes from 301 individuals (287 moraceous species, 42 genera). The primary objectives of this study were to: (1) Resolve phylogenetic relationships among tribes, genera, and sections; (2) Identify the underlying biological processes contributing to phylogenetic conflicts; and (3) Propose taxonomic revisions based on the phylogenetic results and morphological evidence.

2. Materials and methods 2.1. Taxon sampling

Our nuclear dataset comprised 347 samples (Table S1), including 345 newly sequenced samples and two publicly available complete nuclear genomes from GenBank (Artocarpus heterophyllus: SRR2141198 and Artocarpus nanchuanensis: SRR11659666). Of these, 326 samples represented 319 moraceous species, covering all seven recognized tribes and 45 of the 48 genera (excluding Noyera, Pseudostreblus, and Sloetia; Table 1), with extensive coverage of subgenera, sections, and series within the family. Each species of Moraceae was represented by a single sample, except for seven monotypic taxa (genus or section), each represented by two samples. The remaining 21 samples represented outgroup species from six other families of Rosales (i.e., Cannabaceae, Elaeagnaceae, Rhamnaceae, Rosaceae, Ulmaceae, and Urticaceae). Leaf material was obtained from herbaria (CAS, KUN, MO, NY, OS, TEX; acronyms follow Index Herbariorum) and from silica-gel-dried field collections. Detailed sampling information is provided in Table S1.

2.2. Library preparation and sequencing

Total genomic DNA was extracted from herbarium or silica-dried tissues using a modified CTAB protocol (Doyle and Doyle, 1987). Sequencing libraries were enriched for 100 low-copy or single-copy nuclear genes designed for the nitrogen-fixing clade, which have been successfully applied across several families and genera within the clade, including Moraceae and other related taxa (e.g., Fu et al., 2023; Yang et al., 2023; Tian et al., 2024). DNA quantification, library preparation, hybrid enrichment, and Illumina sequencing (150-bp paired-end reads) were performed by RAPiD Genomics (Gainesville, Florida, USA).

2.3. Nuclear sequence assembly and cleanup

Adapters and low-quality bases were trimmed from raw reads using Trimmomatic v.0.36 (Bolger et al., 2014) using the following parameters: fa: 2:30:10:8:TRUE. Cleaned reads were assembled into the 100 target genes using HybPiper v.1.3.1 (Johnson et al., 2016), with Arabidopsis thaliana amino acid sequences as references. Potential paralogs were rigorously filtered by HybPiper (Johnson et al., 2016) using a two-step ortholog selection procedure. First, paralog detection was conducted with the script exonerate_hits.py, which identifies contigs covering ≥ 85% of the reference sequence as full-length candidates. When multiple contigs were recovered for a given gene, sequencing coverage depth was compared, and contigs with ≥ 10× higher coverage than others were provisionally retained as putative orthologs. Second, when coverage was comparable among contigs, the sequence showing the highest identity to the reference was selected as the ortholog, as it best matched the target locus.

Two datasets were generated to assess the impact of including non-coding sequences on phylogenetic analyses. First, the dataset including exon-only sequences for each gene was extracted via the HybPiper script "retrieve_sequences.py" (https://github.com/mossmatters/HybPiper) (the "FNA" dataset). Second, the dataset including exons with flanking non-coding sequences for each gene was generated using the HybPiper script "intronerate.py" (the "supercontig" dataset).

To alleviate the effect of missing data on phylogenetic analysis, we first filtered out sequences shorter than 20% of the mean length for each gene in each dataset using the shell script "trimbylength.sh" (Gardner et al., 2021b; https://datadryad.org/dataset/doi:10.5061/dryad.v41ns1rrt). The remaining sequences for each gene were aligned using MAFFT v.7.305b (Katoh and Standley, 2013). Columns containing more than 80% gaps in each gene alignment were removed using TrimAl (Capella-Gutierrez et al., 2009). To further reduce the effect of long-branch attraction in phylogenetic analysis, long branches were flagged and removed for each gene in ETE3 (Huerta-Cepas et al., 2016) following the pipeline of KewHyb-SeqWorkshop (Baker et al., 2022) with parameters inlen = 0.30, outlen = 0.70, and leaflen = 0.30. Gene alignments with fewer than 40 samples were excluded from subsequent analyses. The final FNA and supercontig datasets each consisted of 93 genes.

2.4. Plastid sampling, assembly, and cleanup

To maximize taxon sampling in the plastid dataset, we assembled plastid loci from four sources: (1) 152 Hyb-Seq samples from this study (195 samples from which plastid loci could not be reliably assembled were excluded); (2) 70 genome skimming samples (~2 Gb) sequenced following the same approach described in Gu et al. (2024) at the Molecular Biology Experimental Center, Kunming Institute of Botany, Chinese Academy of Sciences; (3) 25 complete plastome sequences obtained from GenBank; and (4) 73 Hyb-Seq samples from Gardner et al. (2021a).

Among the 320 plastid samples, 301 represented 287 species of Moraceae, covering all seven tribes and 42 genera (except Bosqueiopsis, Hijmania, Noyera, Scyphosyce, Sloetia, and Treculia; see Table 1). The remaining 19 samples represented 19 outgroup species from six other families within Rosales, including all outgroups present in the nuclear dataset except Dryas grandis and Ziziphus rhodoxylon. Each species was represented by a single sample, except for monotypic Moraceae genera, each of which was represented by two samples (Table S2).

Plastomes were assembled using GetOrganelle v.1.7.4.0 (Jin et al., 2020) and annotated using PGA (Qu et al., 2019; Zhang et al., 2025), with the plastomes of Ficus altissima (NC_053895.1) and Morus notabilis (NC_027110.1) as references. Subsequently, 90 plastid protein-coding sequences (CDSs) were extracted from the annotated plastomes using the Python script "get_annotated_regions_from_gb.py" (Zhang et al., 2020; https://github.com/Kinggerm/PersonalUtilities/). For samples for which complete plastomes could not be assembled using GetOrganelle, we applied the same methods described for the nuclear dataset to assemble the 90 CDSs from cleaned plastid reads, using sequences from Ficus altissima and Morus notabilis as references. Cleanup and alignment of plastid loci were performed using the same procedures as those used for the nuclear dataset. The final plastid dataset consisted of 90 CDSs representing 80 plastid genes.

2.5. Phylogenetic analysis

Maximum likelihood (ML) phylogenetic analyses, with or without a partitioning scheme, were conducted for the concatenated alignments of two nuclear (i.e., FNA and supercontig) and one plastid dataset using RAxML v.8.2.11 (Stamatakis, 2014) with the GTRGAMMA model and 100 bootstrap (BS) replicates. The best partitioning strategies were determined by PartitionFinder v.2 (Lanfear et al., 2017) with the rcluster algorithm under the AICc criterion. Coalescent-based species trees were inferred for the two nuclear datasets using ASTRAL v.5.6.3 (Mirarab et al., 2014; Zhang et al., 2018), with 93 single-gene ML trees for each dataset being estimated by the same ML analysis. Node support of the ASTRAL tree was measured by local posterior probability (LPP). In total, seven phylogenetic trees were inferred: two ML trees without a partitioning scheme (FNA_ML and supercontig_ML), two ML trees with a partitioning scheme (FNA_ML + P and supercontig_ML + P), two species trees inferred using ASTRAL (FNA_ASTRAL and supercontig_ASTRAL) based on the two nuclear datasets, and a ML tree based on the plastid dataset (plastid tree), respectively. Trees were visualized and manipulated using FigTree v.1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/) and Dendroscope v.3 (Huson and Scornavacca, 2012).

2.6. Phylogenetic hypothesis tests

Six deep nodes exhibiting strong topological conflicts were identified among the seven aforementioned trees. These conflicting topologies involved the positions of the following clades: 1) Chlorophoreae, 2) Bagassa, 3) Sloetiopsis, 4) Artocarpus subg. Prainea, 5) Ficus subg. Sycomorus, and 6) Artocarpus subg. Artocarpus sect. Duricarpus ser. Laevifolii (Figs. 1–Fig. 3). For each conflicting node, we used the approximately unbiased (AU) test (Shimodaira, 2002) to determine whether the focal topology was statistically favored by the FNA dataset. To reduce the impact of other discordant topologies, all branches except the target node were fixed in each AU test. We calculated per-site log-likelihoods for each topology using RAxML under the GTRGAMMA model and then performed the AU test in CONSEL v.1.6 (Shimodaira and Hasegawa, 2001).

We used PhyParts (Smith et al., 2015) to quantify concordance and conflict among gene trees. Specifically, we mapped the 93 gene trees from the FNA dataset onto the FNA_ASTRAL tree and considered branches with BS < 50% in the gene trees to be uninformative. Results were visualized using the Python script "Phypartspiecharts.py" (https://github.com/mossmatters/MJPythonNotebooks/blob/master/phypartspiecharts.py).

To distinguish strong conflict from weakly supported branches, we evaluated the phylogenetic information across the species tree using Quartet Sampling (QS; Pease et al., 2018) with 1000 replications and the concatenated alignment of the FNA dataset as input. The QS method is a quartet-based evaluation system that estimates the confidence, consistency, and informativeness of internal relationships by randomly subsampling quartets that are subsampled from the input species tree and concatenated supermatrix (Pease et al., 2018). The QS results provide four values for each internal node: Quartet Concordance (QC), Quartet Differential (QD), Quartet Informativeness (QI), and Quartet Fidelity (QF). The values of QC, QD, and QI were plotted on the FNA_ASTRAL tree using the R script "plot_QC_ggtree.R" (https://github.com/ShuiyinLIU/QS_visualization).

2.7. Phylogenetic signal

To further investigate the gene-wise and site-wise phylogenetic information for each uncertain node, we quantified phylogenetic signals for the six sets of conflicting topologies using the methods of Shen et al. (2017). Phylogenetic signal was calculated as the difference in log-likelihood scores between alternative nodal resolutions (two or three topologies) for a given node in a phylogenetic tree. We first calculated site-wise log-likelihood scores (SLS) and gene-wise log-likelihood scores (GLS) for each alternative topology of a given node based on the FNA dataset. We then computed the difference in site-wise log-likelihood scores (ΔSLS) among topologies for each site, the difference in gene-wise log-likelihood scores (ΔGLS) among topologies for each gene, and the average ΔSLS value for each gene.

Two methods were used to detect "outlier" genes that may have disproportionately influenced phylogenetic inference due to high phylogenetic signal. In the first method, we defined the gene with the highest phylogenetic signal (i.e., highest ΔGLS) as an outlier gene using the Perl script "Strong_genes_selector.pl" (Shen et al., 2017; https://figshare.com/articles/dataset/Contentious_relationships_in_phylogenomic_studies_can_be_driven_by_a_handful_of_genes/3792189). In the second method, average ΔSLS values of genes were assumed to follow a Gaussian-like distribution. Genes with average ΔSLS values outside the two-tailed bounds (lower or upper) of this distribution were defined as outliers (Zhang et al., 2020). Using the standard deviation, the bounds were defined as:

\text{upper bound} =\min (\max (\mathrm{x}), {\rm{\mathsf{μ}}}+3 {\rm{\mathsf{σ}}})\\ \text{lower bound} =\max (\min (\mathrm{x}), {\rm{\mathsf{μ}}}-3 {\rm{\mathsf{σ}}})

where x is the set of average ΔSLS values, and μ and σ denote the mean and standard deviation of this set, respectively.

For each focal node, two reduced datasets were generated by excluding outliers identified by each method. The ML and ASTRAL trees were then reconstructed based on these reduced datasets using the same approach as described above.

2.8. Coalescent simulations and phylogenetic networks

We conducted coalescent simulations to test whether ILS alone was sufficient to explain the nuclear–nuclear and cytonuclear discordances following previous studies (Garcia et al., 2017; Stull et al., 2020). Two simulations were performed: 1. We first examined potential nuclear–nuclear discordance by simulating 1000 nuclear gene trees under the coalescent with the FNA_ASTRAL tree as the guide tree using the Python package DendroPy v.4.1.0 (https://github.com/jeetsukumaran/DendroPy; Sukumaran and Holder, 2010). We then mapped the simulated trees onto the nuclear alternative topologies of these six uncertain clades to calculate the clade probabilities using RAxML ("-f b" option). 2. We tested for potential sources of cytonuclear discordance by pruning both the FNA_ASTRAL tree and the plastid tree to retain 209 shared species, including Urtica fissa as the sole outgroup, to minimize phylogenetic interference. The branch lengths of the pruned FNA_ASTRAL tree were separately scaled by factors of two and four to approximate the expected branch lengths of an organellar tree (McCauley, 1994). We then simulated 1000 plastid gene trees using the same approach, with the scaled FNA_ASTRAL tree as the guide tree. The simulated trees were then mapped onto the plastid tree to calculate clade probabilities. Under the scenario of ILS, any relationship observed in the actual nuclear or plastid tree should appear with high or non-negligible probability among the simulated trees. Conversely, when introgression is present, relationships exclusive to the phylogenetic tree would be observed at a very low probability or may be absent (Folk et al., 2017; Garcia et al., 2017; Stull et al., 2020).

To further explore the reticulate evolutionary history of Moraceae, we inferred phylogenetic networks from individual gene trees of the FNA dataset using the "InferNetwork_MPL" command (Yu and Nakhleh, 2015) in PhyloNet v.3.8.2 (Than et al., 2008; Wen et al., 2018). To generate taxon-reduced datasets for each focal node, we selected one or two species representing major clades relevant to the branch of interest. This approach addressed computational restrictions while focusing on potential reticulation associated with six phylogenetically uncertain nodes (Chlorophoreae, Bagassa, Sloetiopsis, Artocarpus subg. Prainea, Artocarpus subg. Artocarpus sect. Duricarpus ser. Laevifolii, and Ficus subg. Sycomorus). In total, six taxon-reduced datasets were employed, excluding gene trees lacking outgroups.

Network analysis was conducted with five network searches, allowing for the maximum number of reticulations ranging from one to five. We performed 10 independent runs for each search, collapsing gene tree nodes with a BS value below 50%. We also calculated the likelihood scores for bifurcating FNA_ML and FNA_ASTRAL trees using the "CalGTProb" command (Yu et al., 2012). AIC values were calculated to determine the optimal number of possible reticulations, with the number of parameters being equal to the number of branch lengths plus the number of estimated inheritance probabilities (γ).

3. Results 3.1. Nuclear and plastid datasets assembly

After filtering, the number of nuclear genes per sample in the FNA dataset ranged from seven (Ampalis dimepate) to 90 (Maquira guianensis), and 226 of 347 (65.13%) samples had more than 80 genes. The concatenated alignment of the FNA dataset had 97, 917 bp, of which 64.2% were parsimony-informative. The number of genes per sample in the supercontig dataset ranged from six (A. dimepate) to 90 (Ficus prostrata and M. guianensis) nuclear genes, with 240 of 347 (71.21%) samples having more than 80 genes. Its concatenated alignment was 144, 610 bp, of which 57.9% were parsimony-informative. The final plastid dataset contained 90 plastid CDSs; its concatenated alignment was 69, 545 bp, of which 56.1% were parsimony-informative.

3.2. Nuclear phylogeny and gene trees discordances

Our nuclear trees fully supported the monophyly of Moraceae and its seven tribes (BS = 100, LPP = 1; Figs. 1, 2 and S1–S5). Of the 45 sampled genera, 15 were represented by one taxon, whereas 30 were represented by two or more taxa. Among the 30 genera with sampling percentages ranging from 18.2% to 100% (Table 1), 27 were monophyletic, while the remaining three (Batocarpus, Clarisia, and Dorstenia) were non-monophyletic. In addition, 14 sampled subgenera and 37 sampled sections were represented by two or more taxa, and one sampled subgenus and eight sampled sections were represented by one taxon. Our nuclear trees consistently supported the monophyly of 11 of 14 subgenera and 19 of 37 sections with multiple taxa (Figs. 1, 2 and S1–S5), and the remaining subgenera and sections were non-monophyletic (e.g., Brosimum subg. Trymatococcus, Ficus subg. Ficus, Dorstenia sect. Dorstenia, and Artocarpus sect. Duricarpus).

Similar tribal, generic, subgeneric, and sectional relationships were resolved for Moraceae based on different nuclear datasets (FNA vs. supercontig), different phylogenetic methods (concatenation vs. coalescent), and with or without partitioning schemes for the concatenated alignment in ML analyses (Fig. 1 and S1–S5). However, several conflicting topologies were observed at or above the sectional level in the six nuclear trees (Figs. 1 and S1–S5). For example, the placement of Artocarpus subg. Prainea exhibited strong conflict, with two ASTRAL trees supported subg. Prainea as sister to the clade comprising subgenera Artocarpus and Pseudojaca (LPP = 0.80, 0.92, Figs. 1 and S3), whereas all four ML trees supported subg. Pseudojaca as sister to the clade comprising subgenera Artocarpus and Prainea (BS = 95–100, Figs. 2, S1, S2, S4 and S5). Poor gene sampling likely contributed to other conflicts, as seen in Ampalis, with only six to seven genes recovered from Hyb-Seq data, it was either supported as sister to Paratrophis in two ASTRAL trees and two ML trees based on the FNA dataset (BS = 77, 81, LPP = 0.42, 0.49; Fig. 1 and S3–S5), or as sister to the remaining Moreae excluding Taxotrophis in the two ML trees based on the supercontig dataset (BS = 54, 65; Figs. 2, S1 and S2).

All ML trees strongly supported the Bosqueiopsis + Trilepisium clade, the Scyphosyce group (i.e., ((Scyphosyce, (Dorstenia africana, D. kameruniana)), (D. djettii, D. oligogyna))), and the Hijmania + Utsetela clade as successive sisters to Dorstenia s.s. (i.e., this genus excluding the aforementioned four species), with BS values of 100, 93–99, and 100 for three successive nodes, respectively (Figs. 2, S1, S2, S4 and S5). An alternative topology with low support was revealed by two ASTRAL trees, in which a clade comprising Dorstenia s.s. and the Hijmania + Utsetela clade was sister to the clade comprising the Scyphosyce group and the Bosqueiopsis + Trilepisium clade (LPP = 0.46, 0.59; Figs. 1 and S3). Within Ficus, a sister relationship between Ficus subg. Oreosycea and subg. Spherosuke (= subg. Urostigma) was supported by three of four ML trees (i.e., FNA_ML, supercontig_ML, and supercontig_ML + P trees; BS = 45–72; Figs. 2, S1, S4 and S5). In contrast, subg. Oreosycea was resolved as sister to the gynodioecious clade (clade name according to Rasplus et al., 2021) comprising subg. Ficus, subg. Terega (= subg. Sycidium), subg. Sycomorus, and subg. Synoecia; this relationship was supported in the FNA_ML + P tree and two ASTRAL trees (BS = 53, LPP = 0.82, 0.83; Figs. 1, S2 and S3).

PhyParts analyses revealed extensive gene tree discordance regarding the phylogenetic positions of Chlorophoreae, Bagassa, Sloetiopsis, Artocarpus subg. Prainea, Artocarpus ser. Laevifolii, and Ficus subg. Sycomorus (Figs. 1 and S6), and their phylogenetic relationships in the FNA_ASTRAL tree received weak to strong QS support (Fig. S7). Specifically, the sister relationship between Bagassa and the Morus + Sorocea + Trophis clade, and that between Ficus subg. Sycomorus and the remaining gynodioecious clade were strongly supported in the QS analysis (Fig. S7). However, the sister relationships between Chlorophoreae and the MA clade (i.e., Moreae + Artocarpeae), between Sloetiopsis and a clade comprising Bosqueiopsis, Brosimum, Dorstenia, Hijmania, Scyphosyce, Treculia, Trilepisium, and Utsetela, between Artocarpus subg. Prainea and subg. Artocarpus, as well as between Artocarpus ser. Laevifolii and sect. Artocarpus, received weak support in the QS analysis (Fig. S7).

3.3. Plastid phylogeny of Moraceae

In the plastid tree, six of seven tribes, and 27 of 30 genera with more than one sample (the other 12 sampled genera contained only one sample) were recovered as monophyletic (Figs. 2 and S8). The non-monophyletic tribe was Parartocarpeae. The non-monophyletic genera were Dorstenia, Naucleopsis, and Pseudolmedia. Most subgenera and sections were monophyletic, except for one subgenus of Brosimum (i.e., subg. Helianthostylis), one section of Artocarpus (sect. Duricarpus), and most subgenera and sections of Dorstenia and Ficus.

The plastid tree generally supported a similar backbone for tribal and generic relationships within Moraceae as those in the nuclear trees (Fig. 2). However, several cytonuclear discordances were observed at the tribal and generic levels, and high levels of cytonuclear discordance were observed at the subgeneric level, particularly in Ficus (Figs. 2 and S8). Chlorophoreae was sister to a clade comprising two of three genera in Parartocarpeae (Hullettia and Parartocarpus; BS = 96; Figs. 2 and S8), but was sister to a clade comprising Artocarpeae and Moreae in the nuclear tree. The third genus (Pseudostreblus) of Parartocarpeae was sister to a clade comprising Dorstenieae, Antiarideae and Ficeae (hereafter, the DAF clade; BS = 99; Figs. 2 and S8), resulting in a non-monophyletic Parartocarpeae. Other instances of genus-level cytonuclear discordance were found in the phylogenetic positions of one genus (Bagassa) in Moreae and four genera (Helicostylis, Maquira, Naucleopsis, and Pseudolmedia) in Antiarideae. Subgenus-level cytonuclear discordance was observed in the phylogenetic positions of nearly all subgenera of Ficus. Cytonuclear discordance was also observed for the positions of several Ficus sections such as F. subg. Pharmacosycea sect. Oreosycea.

3.4. Phylogenetic signal of six sets of phylogenetic hypotheses

The AU tests indicated that topologies derived from our nuclear trees generally received stronger support than those from the plastid dataset or previous studies for the placements of Chlorophoreae, Bagassa, Sloetiopsis, Artocarpus subg. Prainea, Artocarpus ser. Laevifolii, and Ficus subg. Sycomorus (Fig. 3). Because the AU test evaluates plastid-derived topologies only under the nuclear dataset, it cannot directly refute the validity of the plastid tree itself. However, the distribution of gene-wise and site-wise phylogenetic signal showed that all alternative topologies of these six uncertain clades were supported by subsets of genes and sites (Fig. 3). For the position of Chlorophoreae, nuclear topology A1 (i.e., sister to the MA clade) was favored, supported by the highest proportions of genes (51.6%) and sites (67.6%), followed by nuclear topology A2 (i.e., sister to the PDAF clade) and plastid topology A3 (sister to the Parartocarpus+Hullettia clade), which were supported by 33.3% and 15.1% of genes, and 20.6% and 11.8% of sites, respectively (Fig. 3A). Regarding the position of Bagassa, 48.4% of genes and 22.4% of sites supported plastid topology B3 (sister to the Morus + Trophis clade), while 33.3% of genes and 46.1% of sites supported nuclear topology B2 (sister to Sorocea), and the remaining genes and sites supported nuclear topology B1 (i.e., sister to (Sorocea, (Morus, Trophis))) (Fig. 3B). For the position of Sloetiopsis, similar proportions of genes (49.5% and 50.5%) supported nuclear topologies C1 (i.e., sister to ((Brosimum, Treculia), ((Bosqueiopsis, Trilepisium), (Scyphosyce group, (Dorstenia s.s., (Hijmania, Utsetela)))))) and C2 (i.e., sister to Brosimum + Treculia) (Fig. 3C), with the latter topology had higher site-wise support (63.9%). For the position of Artocarpus subg. Prainea, nuclear topology D1 (sister to subg. Artocarpus) was supported by 32.3% of genes and 66.4% of sites, while nuclear topology D2 (sister to subg. Artocarpus + subg. Pseudojaca) was supported by 67.7% of genes and 33.6% of sites (Fig. 3D). For the position of Artocarpus ser. Laevifolii, nuclear topology E1 (sister to sect. Artocarpus) was supported by 37.6% of genes and 82.5% of sites, with the remaining genes and sites supporting nuclear topology E2 (sister to ser. Asperifoli) (Fig. 3E). For the position of Ficus subg. Sycomorus, nuclear topology F1 (sister to the remaining gynodioecious clade) had lower proportions of supporting genes (48.4%) but higher proportions of supporting sites (84.4%) compared to the nuclear topology F2 (sister to (subsect. Frutescentiae, (sect. Eriosycea, subg. Synoecia)) (Fig. 3F).

Additionally, the results of phylogenetic analyses based on reduced nuclear datasets showed that removing outlier genes did not alter the topologies for the six uncertain nodes, regardless of the method used to define the outlier gene (Figs. S11–S34).

3.5. Coalescent simulations and network analyses

Coalescent simulations indicated that clade probabilities for all proposed topologies of the six clades with uncertain positions in the nuclear tree were low to moderate (Figs. 3 and S35–S40). Specifically, nodes related to Chlorophoreae, Bagassa, and Sloetiopsis showed moderate clade probabilities (16–35%), suggesting the presence of ILS. Clade probabilities for nodes related to Artocarpus subg. Prainea, Artocarpus ser. Laevifolii, and Ficus subg. Terega were near zero (1–4%), indicating that ILS alone was insufficient to explain these nuclear–nuclear discordances.

When simulated trees were mapped onto the empirical plastid tree, similar patterns of clade probabilities across the tree were recovered when the guide tree multiplied by 2 (Fig. S42) or 4 (Fig. S43). Regarding deep cytonuclear nodes, the clade probability for the Chlorophoreae + Parartocarpeae clade was near zero (1% or 4%), suggesting the conflicting plastid topology was not within the ILS prediction, whereas the clade of (Bagassa, (Morus, Trophis)) showed moderate clade probabilities (25% or 27%), which were within the ILS prediction (Fig. 3A and B). Different patterns of clade probabilities emerged across different tribes of Moraceae (Figs. S36 and S37). The plastid topologies that conflicted with the nuclear tree had clade probabilities of zero or near zero in the whole Ficeae, Antiarideae, Chlorophoreae, and most of Dorstenieae (particularly in Dorstenia), whereas the conflicting topologies had a moderate or non-negligible clade probability in most of Moreae and Artocarpeae.

PhyloNet results suggested that any network involving one to five reticulation events fitted six taxon-reduced datasets better than a binary tree, based on the AIC, AICc, and BIC scores (Table S3). Optimal phylogenetic networks (Fig. 4) revealed extensive reticulations among major lineages of Moraceae, including among genera of Moreae (Fig. 4B), among subgenera of Ficus (Fig. 4F), and among sections (or series) of Artocarpus (Fig. 4E).

For the six uncertain clades that were our focus in this study, the optimal network suggested no direct reticulation event between Chlorophoreae and the Parartocarpus + Hullettia clade. Instead, it indicated a reticulation event between the Chlorophoreae + MA clade and the Parartocarpus + Hullettia clade, with a high γ value (0.446; Fig. 4A). Two hybridization events with nearly equal γ values were detected: (1) between A. ser. Asperifoli and the ancestor of A. ser. Laevifolii (0.435 vs. 0.565; Fig. 4E), and (2) between A. ser. Laevifolii and the ancestor of A. ser. Asperifoli (0.423 vs. 0.577; Fig. 4E).

Reticulation events with highly unbalanced γ values from two parental lineages were inferred in relation to Bagassa (Fig. 4B), Sloetiopsis (Fig. 4C), Artocarpus subg. Prainea (Fig. 4D), and Ficus subg. Sycomorus (Fig. 4F). For instance, potential hybridization was observed between the stem lineages of Bagassa and Trophis (γ = 0.243; Fig. 4B), and between the ancestor of Sloetiopsis and a clade comprising Bosqueiopsis, Dorstenia, Scyphosyce, Trilepisium, and Utsetela (γ = 0.170; Fig. 4C). Hybridization was detected between the clade of F. subg. Sycomorus and the clade comprising F. subg. Ficus sect. Eriosycea and F. subg. Synoecia, with an inheritance probability of 0.334 (Fig. 4F).

4. Discussion 4.1. Insights into phylogenetic relationships in Moraceae based on nuclear data

Our phylogenetic analyses indicated that neither analytical approach nor the inclusion of non-coding regions in the nuclear dataset significantly impacted the inference of tribal, generic, subgeneric, and sectional relationships in Moraceae (Figs. 1 and S1–S5). This finding was consistent with previous studies examining the influence of methodologies on phylogenetic inference, as exemplified by studies in Artocarpus (Moraceae; Gardner et al., 2021c), Fagaceae (Liu et al., 2025), and angiosperm-wide analyses (Hu et al., 2023). Nevertheless, methodological differences can lead to topological conflicts in specific cases, as observed in Artocarpus subg. Prainea where concatenation and coalescent approaches strongly supported conflicting topologies, with one receiving higher site-wise support and the other greater gene-wise support (Fig. 3D). This conflict underscores a fundamental challenge in phylogenomics: concatenation methods, which optimize the phylogenetic signal across all nucleotide sites, can be misled by strong but misleading signal concentrated in a subset of sites, while coalescent methods, which aggregate individual gene tree histories, are more sensitive to widespread but stochastic gene tree discordance caused by processes like ILS (Folk et al., 2017; Shen et al., 2021; Hughes et al., 2023). The case of Artocarpus subg. Prainea is a prime example where the dominant site signal and the dominant gene tree signal tell different stories, a pattern characteristic of ancient hybrid lineages where a majority of the genome has been homogenized by introgression from one parent, while a smaller subset of genes retains the signature of the other.

While our study recovered a well-resolved backbone of Moraceae that was largely congruent with previous nuclear phylogenies (e.g., Zerega and Gardner, 2019; Gardner et al., 2021a, 2023), strongly supported topological conflicts persist at several deep nodes across nuclear trees from earlier studies (e.g., Zerega and Gardner, 2019; Gardner et al., 2021a, 2021b, 2023) and our six nuclear trees (Fig. 3). In addition to previously documented incongruences regarding the phylogenetic positions of Chlorophoreae (Zerega and Gardner, 2019 vs. Zhang et al., 2019b), Bagassa (Zerega et al., 2005; Clement and Weiblen, 2009; Williams et al., 2017 vs. Gardner et al., 2021a, 2023), and Artocarpus subg. Prainea (Zerega et al., 2010; Williams et al., 2017 vs. Gardner et al., 2021c), our expanded taxonomic and genomic sampling newly revealed nuclear–nuclear conflicts concerning the placements of Sloetiopsis, Ficus subg. Sycomorus, and Artocarpus ser. Laevifolii.

Notably, these persistent deep topological conflicts within Moraceae are prime candidates for investigation using quantitative approaches capable of elucidating their evolutionary basis, such as ILS and hybridization. Integrative strategies combining morphological evidence with genomic datasets, supplemented by AU tests and phylogenetic signal quantification across alternative topologies, provide a powerful framework to resolve these contentious nodes, as successfully demonstrated in metazoans (Shen et al., 2017), Fabaceae (Zhang et al., 2020), and Fagales (Yang et al., 2021).

Quantitative analyses indicate that the phylogenetic resolutions placing Chlorophoreae as sister to the MA clade and Bagassa as sister to Sorocea—supported by Zerega and Gardner (2019), Gardner et al. (2021a), and our study (Figs. 1, 2 and S1–S5)—receive relatively strong support from both AU tests and gene-/site-wise phylogenetic signals (Fig. 3A and B). These results suggest that alternative topologies reported in earlier studies likely represent historical signals that have been overwritten or outcompeted in most of the genome but remain detectable in smaller datasets. Further support comes from a shared morphological synapomorphy—straight stamens in both Bagassa and Sorocea—which corroborates the sister relationship between these genera. In contrast, support for the alternative topologies of Chlorophoreae and Bagassa, reported in earlier studies based on limited loci (Clement and Weiblen, 2009; Chung et al., 2017; Gardner et al., 2017; Williams et al., 2017; Zhang et al., 2019b), likely reflects ILS. These coalescent simulations suggest that ILS may contribute substantially to the discordance (Fig. 3A and B), but alternative processes cannot be ruled out. The fact that removing potential "outlier" genes did not alter these topological inferences (Figs. S11–S34) further indicates that the discordance is not driven by a few aberrant loci but is a genome-wide phenomenon characteristic of ILS. We acknowledge that coalescent simulations depend on the accuracy of the species tree, and thus our ILS assessments based on the FNA_ASTRAL tree could be sensitive to gene flow or uncertainty in that tree. To address this, we collapsed nodes with bootstrap support below 10%, 30%, and 50% in the individual FNA gene trees to construct alternative ASTRAL trees (Figs. S44–S46). Forthe six conflicting nodes of interest, the topologies of these collapsed trees are consistent with the FNA_ASTRAL tree, indicating that our nuclear ASTRAL tree is highly reliable.

Hybridization, alongside ILS, serves as a key evolutionary process explaining persistent phylogenetic discordance. Among the major gynodioecious clades of Ficus, this is evident in the ambiguous placements of subg. Sycomorus and the subsect. Ficus + subg. Terega clade, which are best explained by introgression between subg. Sycomorus and other gynodioecious lineages (Fig. 3F, 4F). Similarly, hybridization between subg. Ficus subsect. Frutescentiae and subg. Sycidium likely underlies the non-monophyly of sect. Ficus. These findings corroborate a previous report of extensive inter- and intra-subgeneric introgression in Ficus (Wang et al., 2021), highlighting reticulate evolution as a critical process shaping its phylogeny. This pervasive introgression may be facilitated by the shared fig-wasp pollination system, which can act as a conduit for gene flow between otherwise morphologically distinct lineages. We acknowledge that other evolutionary processes, such as selection or gene duplication/loss, cannot be completely excluded as contributors to the observed phylogenetic discordance. Our PhyloNet analyses were based on 93 low-copy nuclear genes, a number sufficient for detecting major reticulation events in plants. While larger nuclear datasets could provide even finer resolution, the inferred networks are robust and provide meaningful insights into hybridization and introgression within Moraceae. Furthermore, the interplay of ILS and hybridization can generate even more complex signatures, as exemplified in Artocarpus sections Artocarpus and Duricarpus. There, strongly supported but conflicting species trees (Figs. 1 and 2), extensive gene tree discordance (Fig. 1), and support from coalescent simulations (Fig. 3E) and network analyses (Fig. 4E) collectively reveal a history shaped by both processes. The detection of bidirectional hybridization with nearly equal inheritance probabilities (γ) between A. ser. Asperifoli and A. ser. Laevifolii (Fig. 4E) suggests a history of extensive gene flow and potential hybrid swarm formation during their divergence, making the reconstruction of a strictly bifurcating tree unrealistic.

Some nodes, however, remain deeply unresolved despite quantitative analysis. The placement of Sloetiopsis represents a notable example, where strongly or moderately supported but conflicting topologies from different analyses receive nearly equal gene support and both pass AU tests (Fig. 3C). This impasse suggests that current data and methods are insufficient to resolve the relationships among Sloetiopsis and its relatives, likely due to ILS as suggested by coalescent simulations (Fig. 3C). This "hard polytomy" signal indicates a potential historical rapid diversification event where ancestral lineages diverged over such a short timescale that genealogical histories failed to coalesce, leaving a persistent, unresolvable signature in the genome (Scherz et al., 2022). Such cases underscore the necessity of incorporating more comprehensive genomic data, including complete genomes, to dissect these phylogenetic challenges.

4.2. Ancient hybridization and incomplete lineage sorting together explain deep cytonuclear discordances

Different genomic compartments may retain distinct evolutionary histories due to variations in inheritance mode, effective population size, and exposure to different evolutionary forces. Although the plastid and nuclear phylogenies were not constructed from exactly the same datasets, this does not affect the discussion of the potential causes of cytonuclear discordance. Similar approaches have been adopted in previous studies investigating cytonuclear discordance (e.g., Stull et al., 2020).

Cytonuclear discordance is especially common in plant phylogenetics, often serving as a clue that complex evolutionary processes—ILS or hybridization—have shaped the genome differently across its compartments. Due to stochastic sorting of ancestral alleles, ILS often produces random discordance across loci and shows symmetrical signature in the genome (i.e., similar frequencies for conflicting topologies) (Edelman et al., 2019; Shang et al., 2025). Our coalescent simulations suggested that some deep cytonuclear discordances in Moraceae fall within the ILS prediction, such as those in Artocarpeae, Moreae and part of Dorstenieae (Figs. S18 and S19). This pattern is further supported by evidence of short internodes and extensive gene-tree discordance at these phylogenetic positions (Figs. 1 and S1–S6). However, similar frequencies of supporting genes and sites were not observed for one of our focal cytonuclear discordances—specifically, the position of Bagassa in Moreae (Fig. 3B)—possibly due to the limited number of loci in our dataset. The presence of ILS in these lineages suggests that they likely experienced episodes of rapid speciation.

Unlike ILS, hybridization and introgression reflect the post-divergence movement of genetic material between lineages, leading to asymmetric patterns of allele sharing (Pease and Hahn, 2015). Cytoplasmic (e.g., mitochondrial or chloroplast) capture is one of the prominent examples of introgression. Because organellar genomes are usually maternally inherited and non-recombining, they can introgress into a species more easily than nuclear genomes through hybridization followed by repeated backcrossing (Larson et al., 2026). This uniparental inheritance and lack of recombination mean that a single successful hybrid crossing can lead to complete chloroplast replacement, whereas nuclear introgression is typically piecemeal and limited by selection against incompatible alleles. This results in a mismatch of evolutionary history between the nuclear DNA (reflecting the host species) and the organellar DNA (reflecting the donor species), leading to cytonuclear discordance. Recent phylogenomic studies increasingly identify ancient hybridization as the major source of deep cytonuclear discordance, which is often structured across geographic regions, reflecting past hybrid zones or secondary contact (Acosta and Premoli, 2010; Liu et al., 2025). Our study and previous studies (Wang et al., 2021; Gardner et al., 2021a, 2023) revealed similar cases across Moraceae lineages including Antiarideae, Chlorophoreae, and particularly Dorstenia and Ficus. Both coalescent simulations (Figs. 3A, S42 and S43) and network results (Fig. 4) suggest that ancient hybridization may have played a more prominent role than ILS as the primary driver of plastid-nuclear phylogenetic discordances in these lineages. For example, historical hybridization between the Hullettia+Parartocarpus clade and the Chlorophoreae + MA clade, as supported by network analysis (Fig. 4A) and their overlapping geographic distribution in Asia (Gardner et al., 2023), could have potentially contributed to the non-monophyly of Parartocarpeae observed in the plastid tree.

These findings collectively challenge the reliability of plastid data for resolving intertribal, intergeneric, and intrageneric relationships (e.g., Ficus and Dorstenia) in Moraceae. The plastid genome, while valuable for its high copy number and ease of sequencing, should be interpreted with caution in groups with a history of hybridization, as it may represent the history of pollen recipients or a specific maternal lineage rather than the species' primary evolutionary trajectory. Nevertheless, a comprehensive understanding of the evolutionary history of the focal group requires integrating both nuclear and plastid data. The discordance itself becomes a rich source of historical information, revealing past contact zones and hybridization events that shaped biodiversity. Furthermore, the signature of ancient hybridization is possibly more easily detectable in the relatively long and non-recombining organellar genomes, which maintain a consistent phylogenetic signal and, due to their highly conserved functions, may introgress more freely across species boundaries with fewer functional incompatibilities than nuclear genes (Larson et al., 2026). Consequently, signals of chloroplast genome introgression are more likely to persist over time even in cases of deep divergence, as seen between the Hullettia + Parartocarpus clade and the Chlorophoreae + MA clade in Moraceae, between Heuchera and Mitella in Saxifragaceae (Folk et al., 2018), and between Lithocarpus and Quercus in Fagaceae (Liu et al., 2025). The extensive introgression in Ficus, as suggested by cytonuclear discordance, implies that hybridization may have played an essential role in its rapid radiation in tropical rainforests during the mid-Miocene (Wang et al., 2021; Gardner et al., 2023). Hybridization could have facilitated adaptation to new ecological niches or pollinator regimes, providing the genetic novelty necessary for explosive diversification. However, although ancient hybridization events have been inferred, the extent and impact of nuclear introgression remain largely unknown in Moraceae. Future research integrating population-level sampling with whole-genome sequencing will be crucial for addressing key questions—for instance, what proportion of the nuclear genome has been introgressed, which specific genes are involved, and what functional roles these genes may play in the evolutionary success of lineages in Moraceae (e.g., Dorstenia and Ficus). Additionally, studying the co-evolutionary implications of chloroplast capture on plant–pollinator interactions, especially in highly specialized systems like that of Ficus, presents a fascinating avenue for future research.

4.3. Suggested taxonomic changes

Phylogenetic results in the current study (Figs. 1 and S1–S5) and previous studies (e.g., Gardner et al., 2021a, 2021b, 2023) necessitate taxonomic revisions for three non-monophyletic genera (i.e., Batocarpus, Clarisia, and Dorstenia), and some sections and subsections. We advocate a re-evaluation of morphological characters used in previous classifications to assess their congruence with molecular phylogenetic findings, especially in the case of the three non-monophyletic genera and the largest and taxonomically challenging genus, Ficus.

Clarisia and its sister group, Batocarpus, have been consistently recovered as non-monophyletic across multiple phylogenetic studies (Zhang et al., 2019b; Gardner et al., 2021a, 2021c, 2023), including the current study. Considering their shared morphological synapomorphies, including spike-like male inflorescences and globous-capitate or uniflorous female inflorescences, along with their sympatric distribution in Neotropical regions, we therefore propose subsuming Batocarpus into Clarisia based on integrated phylogenetic and morphological evidence.

The second most species-rich moraceous genus, Dorstenia (ca. 120 species), has been recovered as non-monophyletic in the current study, consistent with previous findings (Misiewicz and Zerega, 2012; Zhang et al., 2019a), with Hijmania, Scyphosyce, and Utsetela nesting within this genus. Zhang Qian et al. (personal communication) proposed merging Hijmania, Scyphosyce, and Utsetela into Dorstenia to restore its monophyly. However, such broad circumscription would render Dorstenia morphologically heterogeneous and obscure the distinctiveness of these three smaller genera, each possessing clear diagnostic features. Instead, we propose to transfer several shrubby Dorstenia species (D. africana, D. djettii, D. kameruniana, and D. oligogyna) to Scyphosyce. This refinement reflects both phylogenetic relationships and growth form (shrubs vs. herbs), and it results in a monophyletic and morphologically coherent Dorstenia composed exclusively of herbaceous species.

The delimitation of sections within Dorstenia requires substantial revisions. Both our study and Zhang et al. (2019a) reject the monophyly of eight out of nine sections recognized by Berg and Hijman (1999)—namely, sects. Aculoma, Kosaria, Dorstenia, Emygdioa, Lecanium, Lomathophora, Nothodorstenia, and Xylodorstenia—with the exception of the monotypic sect. Bazzemia (D. zambesiaca). Vianna-Filho et al. (2021) proposed including all Neotropical Dorstenia species within an expanded sect. Dorstenia, while the non-monophyletic African herbaceous/succulent sections remain unresolved. Based on both current and previous phylogenetic results (Zhang et al., 2019b; Vianna-Filho et al., 2021), we propose recognizing three African clades as separate sections, all of which are successively sister to the Neotropical sect. Dorstenia. The first African clade encompasses most species previously assigned to sect. Aculoma and sect. Kosaria, characterized by succulent or semi-succulent stems and receptacles with two rows of appendages. The second African clade includes most species of sect. Lomathophora, which are rhizomatous herbs with receptacles bearing a broad fringe of appendages. The third African clade consists of a single species, D. elliptica, a rhizomatous subshrub with woody stem bases and receptacles featuring two rows of bracts. This clade was identified in previous studies (Zhang et al., 2019b; Vianna-Filho et al., 2021) but was not sampled in our current dataset. This revised sectional classification provides a phylogenetically informed framework that also corresponds to major morphological and biogeographic discontinuities within the genus.

Berg and Corner (2005) divided Ficus into six subgenera based on morphology, geographic distribution, and earlier studies by Corner (1960a, 1960b, 1960c) and Weiblen (2000). Our findings (Figs. 1 and S1–S5) and recent molecular phylogenetic studies (Rasplus et al., 2021; Gardner et al., 2023) have revealed that four of these subgenera—subg. Spherosuke, subg. Sycomorus, subg. Synoecia, and subg. Terega—are monophyletic, whereas subg. Pharmacosycea and subg. Ficus are not. Despite the non-monophyly of subg. Pharmacosycea, its two constituent sections are monophyletic: sect. Pharmacosycea is sister to the remaining Ficus, and sect. Oreosycea is either sister to the gynodioecious clade or to subg. Spherosuke. Similarly, the non-monophyletic subg. Ficus includes two sections: sect. Eriosycea, which is monophyletic and sister to subg. Synoecia; and sect. Ficus, which is non-monophyletic but contains two monophyletic subsections—subsect. Ficus, sister to subg. Terega, and subsect. Frutescentiae, sister to a clade comprising sect. Eriosycea and subg. Synoecia.

Given the non-monophyly of two Ficus subgenera, a redefinition of subgeneric boundaries is warranted. Two alternative revision strategies are proposed, taking into account monophyly as well as distinct morphological traits and geographic distributions. The first, advocated by Corner (1960a, 1960b, 1960c, 1965) and Gardner et al. (2023), suggests dividing Ficus into a single gynodioecious subgenus—subg. Ficus (formed by merging subg. Sycomorus, subg. Synoecia and subg. Terega)—and three monoecious subgenera: subg. Oreosycea (elevated from sect. Oreosycea of the former subg. Pharmacosycea), subg. Pharmacosycea, and subg. Spherosuke. In this framework, all species within the redefined gynodioecious subgenus are gynodioecious, except for a few monoecious species originally assigned to subg. Sycomorus. This scheme has the advantage of recognizing the major evolutionary split between monoecious and (predominantly) gynodioecious lineages. However, given the large number of species encompassed by this broadly defined subg. Ficus, as well as their considerable morphological and biogeographical diversity, a second, more finely resolved classification may be more practical and informative. In this alternative scheme, we recommend elevating Ficus subg. Ficus sect. Eriosycea, subg. Ficus sect. Ficus subsect. Frutescentiae, and subg. Pharmacosycea sect. Oreosycea to subgenus rank. This classification would recognize five gynodioecious subgenera (Eriosycea, Ficus, Frutescentiae, Synoecia, and Terega), three monoecious subgenera (Oreosycea, Pharmacosycea, Spherosuke), and one mixed subgenus (Sycomorus) that includes both breeding systems.

5. Conclusions

Through integrative analyses of nuclear and plastid data, this study establishes a solid phylogenetic framework for Moraceae and systematically investigates the sources of both nuclear–nuclear and cytonuclear incongruence across genomic compartments and analytical methodologies. Our analyses revealed significant topological conflicts between nuclear and plastid phylogenies, particularly at critical nodes involving rapidly diversified lineages. A combination of ILS and ancient hybridization likely contributes to both deep nuclear–nuclear and cytonuclear conflicts. These findings underscore the complex evolutionary dynamics within Moraceae and highlight the necessity of using coalescent-based approaches for robust phylogenetic reconstruction in this family. We propose several taxonomic revisions to resolve non-monophyletic genera, subgenera, and sections, along with recommendations for further refinement of the family's classification. Overall, our comprehensive phylogeny offers a solid foundation for future studies on the classification, biogeography, and diversification of this ecologically and economically important family.

Acknowledgments

This research was supported by the Yunnan Revitalization Talent Support Program: Yunling Scholar Project (XDYC-YLXZ- 2024-0021), the Basic Research Project of Yunnan Province (No. 202401BC070001), the National Natural Science Foundation of China, key international (regional) cooperative research project (No. 31720103903), the Yunling International High-end Experts Program of Yunnan Province, China (No. YNQR-GDWG-2017-002 and No. YNQR-GDWG-2018-012). We are grateful to the following institutes for providing specimens or silica-dried materials: 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 herbarium of the California Academy of Sciences; the Missouri Botanical Garden; the New York Botanical Garden; the Ohio State University Herbarium; and the University of Texas Herbarium. We are also grateful to Jia-Jin Wu for help with sampling; to Hua-Feng Wang, Diego F. Morales-Briones, Nelson Zamora Villalobos, Rong Zhang, Hui Liu, Si-Yun Chen, Xiao-Gang Fu, Ying–Ying Yang, and Henrique Borges Zamengo for their generous technical support and valuable assistance; to Shuai Liao for assistance in identifying specimens; and to the iFlora High Performance Computing Center of the Germplasm Bank of Wild Species(iFlora HPC Center of GBOWS, KIB, CAS) for computing. We appreciate the photos provided by Qin Tian of the Honghe Tropical Agriculture Institute of Yunnan, China.

Data availability

Raw sequence data have been deposited in the NCBI Sequence Read Archive (BioProject: PRJNA980892). The aligned DNA sequences and trees are available at: https://doi.org/10.6084/m9.figshare.31948566.

CRediT authorship contribution statement

Chen-Xuan Yang: Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis. Shui-Yin Liu: Writing – review & editing, Methodology, Visualization, Software, Formal analysis, Data curation, Investigation. Qin Tian: Writing – review & editing, Resources, Investigation, Data curation. Wei Gu: Writing – review & editing, Methodology, Formal analysis. Qing Lu: Writing – review & editing, Methodology, Formal analysis. Robert P. Guralnick: Writing – review & editing. Gregory W. Stull: Resources, Investigation, Data curation. Heather R. Kates: Resources, Investigation, Data curation. Ryan A. Folk: Resources, Investigation, Data curation. Douglas E. Soltis: Writing – review & editing. Pamela S. Soltis: Resources, Investigation. Elliot M. Gardner: Writing – review & editing, Methodology. Ting-Shuang Yi: Writing – review & editing, Supervision, Resources, Project administration, Investigation, Funding acquisition, Data curation, Conceptualization.

Declaration of competing interest

The authors declare that there are no conflicts of interest regarding the submission of this manuscript, and all authors have approved its publication.

Appendix A. Supplementary data

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