Exploring the systematic positions of two infraspecific taxa of Dioscorea bulbifera based on next-generation sequencing technology, morphology and cytology

Dioscorea bulbifera is a perennial twining herb in the family Dioscoreaceae^[1]. Its tuber, known in traditional Chinese medicine as 'Huangyaozi', has long been used to treat goiter, pharyngitis, pyogenic skin infections, orchitis, and tumors^[2,3]. Recent phytochemical and pharmacological reviews have further reported anticancer, anti-inflammatory, antimicrobial, antidiabetic, antioxidant, and immunomodulatory activities of its extracts and bioactive constituents, particularly steroidal saponins, flavonoids, and diterpenoids^[2−4]. In addition, the species is cultivated as a food crop in parts of Africa and Asia^[2,3] and is grown as an ornamental plant in some regions^[5,6]. Morphologically, D. bulbifera is characterized by ovoid to pyriform tubers, left-twining stems that are green to purplish, and spherical to ovoid bulbils^[1]. Owing to its strong ecological adaptability, D. bulbifera is widely distributed across China, Japan, Korea, India, Myanmar, as well as parts of Oceania and Africa^[1,2]. This broad geographic range has facilitated the accumulation of substantial intraspecific variation, resulting in considerable morphological diversity^[1,2].

The infraspecific taxonomy of D. bulbifera has long been controversial. Since the 19^th century, multiple taxonomic treatments have been proposed based primarily on morphological grounds, yet no consensus has been reached. Prain & Burkill first provided a synopsis of Old World Dioscorea, recognizing D. bulbifera var. kacheo among other infraspecific taxa^[7]; Knuth subsequently recognized eight varieties within D. bulbifera^[8]; and Prain & Burkill later gave a more detailed account of eastern species, further refining infraspecific delimitation, including D. bulbifera var. heterophylla^[9]. Subsequent studies based on molecular and morphological evidence have continued to refine this framework. Using molecular markers, Zheng et al. reported genetic differentiation between samples from mainland China and those from Hainan and Taiwan^[10]. Based on specimen examination and field investigations, Zhou et al. argued for the existence of at least four varieties in China, treating the Hainan entity originally named D. bulbifera var. deltoidea as D. bulbifera var. heterophylla Prain & Burkill, and describing var. brachybotryum Y.Y. Hang & Y.F. Zhou^[11]. These treatments largely document geographically structured morphological variants across the species' wide range^[1,11].

Two historically recognized infraspecific taxa are the focus of the present study. Dioscorea bulbifera var. kacheo Prain & Burkill was first described from the Sikkim Himalaya and adjacent regions to the east^[7]. A lectotype for this name was subsequently designated by Hoque et al.^[12]. It has been retained in regional treatments such as the Flora of Bhutan, where it is diagnosed primarily by stems bearing undulate membranous wings^[13]. In contrast, D. bulbifera f. domestica was originally published as D. sativa f. domestica Makino^[5] and later transferred to D. bulbifera f. domestica^[6]. This taxon has been recognized mainly in Japanese floristic treatments as a cultivated ornamental form characterized by conspicuous variation in bulbil morphology^[6]. These designations follow historical floristic conventions rather than an assessment of evolutionary divergence: in conventional botanical practice, varieties are expected to show geographic, ecological, and/or phylogenetic integrity beyond morphological diagnosability, whereas formae are applied to variants lacking such extramorphological integrity^[14]. Nevertheless, most existing classifications of these taxa still rely heavily on morphological traits, which may yield unstable delimitations when environmental effects are not adequately accounted for^[11]. In particular, robust integrative evidence supporting the recognition of D. bulbifera var. kacheo and D. bulbifera f. domestica remains lacking. Recent low-coverage genome sequencing of Chinese D. bulbifera accessions has revealed geographically structured intraspecific lineages based on plastome data^[15], yet these analyses did not specifically evaluate the taxonomic status of the two taxa.

Cytological evidence, particularly chromosome number and ploidy level, plays a fundamental role in plant systematics and infraspecific delimitation^[16]. Polyploidy is widespread in Dioscorea, and the same species frequently comprises multiple cytotypes; diploid, triploid, tetraploid, and higher-ploidy levels have been documented in cultivated yams such as D. alata and Guinea yams^[17]. For these crops, chromosome counting and flow cytometry provide reliable ploidy estimates^[17]. Cross-ploidy hybridization and introgression are also increasingly recognized as important evolutionary processes in plants^[18]. For D. bulbifera, previous reports likewise indicate marked cytological variation. Asian materials are generally based on x = 10, whereas African and American materials have been associated with x = 9^[16,19]. Somatic chromosome numbers of 2n = 36, 40, 54, 60, 80, and approximately 100 have been recorded in the literature^[16,17,19]. Despite these observations, the relationship between ploidy variation and infraspecific delimitation in D. bulbifera has not been explicitly examined. Whether historically recognized varieties and forms differ in chromosome number or ploidy level therefore remains an open question and provides an important line of evidence for clarifying their taxonomic status.

Advances in molecular systematics have provided powerful tools for resolving complex intraspecific relationships in Dioscorea. Genome-wide SNP data generated by genotyping-by-sequencing (GBS) or DArTseq, often combined with plastome or phenotypic evidence, have been used successfully to reconstruct population structure, domestication history, and cryptic evolutionary relationships that are difficult to resolve from morphology alone. For example, whole-genome resequencing revealed a hybrid origin of white Guinea yam (D. rotundata)^[20], and population genomic analyses have clarified domestication and diversification patterns in African yams^[21]. SNP- and DArTseq-based studies have likewise resolved genetic structure among cultivated yam accessions, including D. alata, D. rotundata, and D. dumetorum^[22]. A high-density genetic map has also been developed for greater yam (D. alata)^[23]. Plastome phylogenomics has long been a cornerstone of plant maternal lineage reconstruction, and in yams it has proven informative for germplasm discrimination and intraspecific phylogeny^[15,24]. For D. bulbifera itself, SSR markers, DArTseq SNPs, and plastome data have revealed substantial genetic diversity and geographically structured lineages among accessions from Brazil, West Africa, and China^[15,25,26]. However, most of these studies were designed for germplasm characterization or breeding rather than taxonomic delimitation. As a result, the evolutionary relationships and taxonomic status of D. bulbifera var. kacheo and D. bulbifera f. domestica remain poorly resolved. Moreover, plastome data alone capture only the maternal genealogy and may be insufficient in groups subject to hybridization or polyploidization^[27], whereas nuclear SNP data provide complementary resolution of nuclear genomic history^[27]. Integrating morphological, cytological, plastome, and nuclear SNP evidence therefore offers a more robust framework for testing infraspecific boundaries in D. bulbifera.

In this study, we employ an integrative approach combining morphological observations, cytological analyses, plastome phylogenomics, and genome-wide nuclear SNP data to address these questions. Specifically, we aim to: (1) assess morphological and cytological differentiation among taxa; (2) reconstruct phylogenetic relationships based on plastome and nuclear genomic data; (3) detect potential gene flow and infer evolutionary processes shaping divergence; and (4) clarify the taxonomic status of D. bulbifera var. kacheo and D. bulbifera f. domestica. Our results provide new insights into the evolutionary history of D. bulbifera and contribute to a more robust framework for its infraspecific classification.

A total of 14 accessions of D. bulbifera var. bulbifera were collected from different geographic regions across China, including Anhui (AHLA), Fujian (FJNJ), Guangdong (GDSG), Guangxi (GXQZ), Henan (HeNSC), Hunan (HuNXX), Jiangsu (JSYX), Sichuan (SCEM), Zhejiang (ZJLS), Chongqing (CQBB), Yunnan (YNLJ), Hainan (HNHK, HNLS), and Taiwan (TWXB). Accession HNLS was used exclusively for chromosome squashing and flow cytometry. With the exception of HNLS, all remaining accessions, together with two D. bulbifera var. kacheo accessions (JCHD1, JCHD2) collected from Yunnan Province and two D. bulbifera f. domestica accessions (YZ1, YZ2) obtained through commercial sources, were subjected to comprehensive analyses (Supplementary Table S1).

Fresh, intact green leaves were collected and immediately preserved by rapid desiccation with silica gel. Voucher specimens corresponding to all accessions were deposited in the Herbarium of the Jiangsu Institute of Botany and the Herbarium of the Chinese Academy of Sciences (NAS) for verification and future reference.

Experimental methods

Ploidy determination by chromosome counting and flow cytometry

Plants of D. bulbifera var. bulbifera, D. bulbifera var. kacheo, and D. bulbifera f. domestica were cultivated in the Dioscorea germplasm nursery of Nanjing Zhongshan Botanical Garden (32°3'32.46'' N, 118°49'51.93'' E). Mature bulbils were collected from October to November and stored in sand until April of the following year. After 3–4 weeks of soil cultivation, fresh lateral root tips were collected between 8:00 and 10:00 a.m. during the period of active mitotic division^[28]. Root tips were pretreated with 0.002 mol·L⁻¹ 8-hydroxyquinoline in the dark for 3 h, fixed in freshly prepared ethanol : glacial acetic acid (3:1) for 24 h at 4 °C (These materials were stored at 4 °C for up to 1–2 weeks before subsequent experiments). After HCl hydrolysis and pectinase treatment, chromosomes were stained with Feulgen's reagent and observed under a microscope at 400×/1,000× magnification^[28]. For each accession, chromosome numbers were determined from 5–10 well-spread metaphase cells in which individual chromosomes could be unambiguously distinguished.

Ploidy levels were further estimated by flow cytometry following published protocols for plant DNA flow cytometry^[29]. Approximately 20 mg of fresh leaf tissue was placed in a plastic Petri dish and chopped with a sharp razor blade in 1 mL of ice-cold woody plant buffer (WPB: 200 mM Tris-HCl, 4 mM MgCl₂, 2 mM Na₂EDTA, 86 mM NaCl, 10 mM sodium metabisulfite, 1% (w/v) PVP, 1% (v/v) Triton X-100, pH 7.5)^[30]. The homogenate was mixed gently by pipetting and filtered through a 42 μm nylon mesh into a sample tube. Nuclei were stained with propidium iodide (PI; 50 μg/mL) and treated with RNase (50 μg/mL)^[29,30]. Samples were incubated on ice for several minutes to 1 h with occasional shaking before analysis. Relative nuclear DNA content was measured by flow cytometry based on the fluorescence intensity of PI-stained nuclei, and ploidy levels were inferred by comparison with standards of known ploidy.

DNA extraction and next-generation sequencing

Genomic DNA was extracted from approximately 50 mg of silica-gel-dried leaf tissue from each of the 17 samples using the DNAsecure Plant Genome DNA Kit (TIANGEN Biotech, Beijing, China). The purity and integrity of the extracted genomic DNA were assessed by agarose gel electrophoresis and spectrophotometry. The DNA was then fragmented using a Covaris E210 ultrasonicator (Covaris Inc., USA), and fragments sized 200–400 bp were selected using the Agencourt AMPure XP-Medium Kit (Thermo Fisher Scientific, USA). Subsequent steps included end repair, A-tailing, adapter ligation, purification, and PCR amplification to construct Illumina sequencing libraries with an insert size of 150 bp. Library quantification and quality control were performed using Qubit 2.0 and Agilent 2100. The precise concentration of the libraries was determined using q-PCR. Finally, high-throughput sequencing was performed on the Illumina NovaSeq 6000 platform.

Plastome assembly, annotation, and whole plastome sequence comparison

Complete plastomes of D. bulbifera were assembled from Illumina paired-end genome skimming data. Raw data were quality-controlled using fastp v0.20.1^[31]. Plastomes were assembled de novo using GetOrganelle v1.7.6^[32] with the parameters: -R 20 -k 21, 51, 91, 121, 131 -F embplant_pt. Plastome sequences of all 17 samples were successfully assembled into complete circular structures. All assembled plastomes were inspected and confirmed to form complete circular molecules with the typical quadripartite structure, including a large single-copy region (LSC), a small single-copy region (SSC), and two inverted repeats (IRs).

The assembled plastomes were annotated using the online software GeSeq^[33] (https://chlorobox.mpimp-golm.mpg.de/geseq) by uploading the sequences and selecting Dioscorea bulbifera as the reference, generating annotation files and high-resolution circular genome maps. To detect potential expansions or contractions in the inverted repeat (IR) regions within D. bulbifera plastomes, a comparison and visualization of the four junctions between the IR and LSC/SSC regions were conducted using IRscope^[34] (https://irscope.shinyapps.io/irapp).

Phylogenetic analyses within Dioscorea bulbifera

Phylogenetic analyses were performed using three datasets: (1) complete plastome sequences; (2) concatenated coding sequences (CDS) of 79 protein-coding genes; (3) concatenated amino acid sequences of the same 79 protein-coding genes (Supplementary Table S2), employing D. nipponica (OQ525997), D. alata (OP787123), and D. zingiberensis (NC_027090) as outgroups. All datasets comprised 20 samples, including three Dioscorea species as outgroups, and all resulting trees were rooted with D. zingiberensis.

Complete plastome sequences were aligned using MAFFT v7^[35], and ambiguously aligned regions were removed using trimAl v1.4 with the automated1 option^[36]. CDS and amino acid sequences were individually aligned and subsequently concatenated into supermatrices based on gene order and sample IDs. Missing data were treated as gap characters ('−') following standard phylogenomic practice.

Maximum likelihood (ML) phylogenetic analyses were performed using IQ-TREE^[37] v3.1.1. The best-fitting substitution models for each dataset were selected using ModelFinder based on the Bayesian Information Criterion (BIC)^[38]. Branch support was evaluated using 1,000 ultrafast bootstrap replicates and 1,000 SH-like approximate likelihood ratio test (SH-aLRT) replicates^[37]. To reduce potential overestimation of bootstrap support values, the -bnni option was applied, performing additional nearest-neighbor interchange (NNI) optimization during bootstrap analysis^[37]. The best-fitting substitution models selected were TVM + F + I for both the plastome and CDS datasets, and Q.BIRD + F + I for the amino acid dataset, as determined by ModelFinder^[38].

Bayesian inference (BI) of the complete plastome dataset was conducted in MrBayes v3.2.7^[39] under the GTR + I + Γ model, which was used because the BIC-selected TVM + F + I model is unavailable in MrBayes. Two independent MCMC runs with four chains each were run for 1,000,000 generations, sampling every 1,000 generations. Convergence was assessed by an average standard deviation of split frequencies < 0.01 (reported by MrBayes) and by effective sample sizes (ESS) > 200 inspected in Tracer v1.7.2^[40]. After discarding the first 25% of sampled trees as burn-in, a majority-rule consensus tree was summarized, and Bayesian posterior probabilities (PP) were calculated. Because ML and BI analyses recovered identical ingroup topologies, and because the CDS and amino acid trees were congruent with the plastome tree, the plastome ML tree is presented with SH-aLRT, UFBoot, and PP support values.

Phylogenetic analysis based on nuclear genome SNPs

Cytological observation revealed ploidy variation among sampled accessions. Because no subgenome-resolved reference genome is available for D. bulbifera, sample-specific ploidy was not incorporated into variant calling; instead, genotype likelihoods and SNP matrices were generated under a unified diploid (diploidized) model for all accessions to enable joint cross-accession analyses. Diploidized genotype calling has been applied in autopolyploid GBS studies when sequencing depth was insufficient for reliable tetraploid dosage estimation^[41]. Under this diploidized coding, heterozygous states are collapsed into a single heterozygous class, and allelic dosage is not recovered. Comparative analyses indicate that diploidized genotypes can retain analytical utility relative to dosage-aware coding, although some information is lost^[42]. Downstream phylogenetic and population-genetic inferences were therefore interpreted primarily in terms of relative placement and differentiation patterns rather than dosage-dependent parameter estimates.

Genome-wide single nucleotide polymorphisms (SNPs) were identified using a standardized Sentieon-based pipeline^[43]. The haplotype-resolved genome of D. bulbifera var. bulbifera (Hainan accession; unpublished) was used as the reference genome. Raw paired-end reads were quality-filtered using fastp v0.23.4^[31] with default parameters to remove adapter sequences and low-quality bases. Clean reads were aligned to the reference genome using Sentieon BWA-MEM (Sentieon Genomics, v202308), and alignments were sorted using Sentieon utilities. PCR duplicates were removed using Sentieon Dedup, and alignment statistics, including coverage depth, GC bias, and insert size distribution were calculated using built-in Sentieon modules. Variant calling was performed using Sentieon Haplotyper in GVCF mode (--emit_mode gvcf, --emit_conf 30, --call_conf 30), followed by joint genotyping across all samples to generate a multisample VCF. Outgroup species (D. nipponica, D. alata, and D. zingiberensis) were processed using the same pipeline.

Biallelic SNPs were filtered using GATK-style hard filters (QD < 2.0, MQ < 40.0, FS > 60.0, SOR > 3.0)^[44] and PLINK v1.9^[45] (--geno 0.15, --maf 0.05). Linkage disequilibrium pruning was performed using PLINK v1.9^[45] (--indep-pairwise 50 10 0.2), retaining approximately independent SNPs.

SNP genotypes were converted into a pseudo-sequence alignment in PHYLIP format using TASSEL v5.2.90^[46], with heterozygous sites encoded using IUPAC ambiguity codes. Maximum likelihood phylogenetic analysis was conducted using IQ-TREE v3.1.1^[37], under the best-fit model (TVM + F + R2) selected by BIC, with branch support assessed using 1,000 ultrafast bootstrap replicates (-bnni).

Population genetic analysis based on nuclear genome SNPs

The 17 in-group samples were assigned to three populations: HD (D. bulbifera var. bulbifera), JC (D. bulbifera var. kacheo), and YZ (D. bulbifera f. domestica). Nucleotide diversity (π) and genetic differentiation (F_ST) were calculated using VCFtools v0.1.16^[47].

Gene flow was assessed using the D-statistic (ABBA–BABA test) implemented in Dsuite v0.4^[48], with the topology (HD, JC, YZ, O), where O represents the outgroup. Statistical significance was evaluated using block jackknife resampling.

For demographic inference, LD pruning was performed using PLINK (--indep-pairwise 50 10 0.1). TreeMix v1.13^[49] was used to infer population splits and migration events, with the outgroup specified as the root. Migration edges (m = 0–5) were tested with 10 independent runs per model. The optimal number of migration events was determined using OptM^[50].

Population structure was analyzed using LD-pruned SNPs. Principal component analysis (PCA) was performed using PLINK v1.9^[45] (--pca). ADMIXTURE v1.3.0^[51] was run for K = 1–8 with cross-validation (--cv), and the optimal K was selected based on the lowest cross-validation error.

The female inflorescences of the three taxa were essentially uniform in morphology (Fig. 1a–c). They are spicate, solitary or several together in the leaf axils, pendulous, with sessile florets closely adpressed to the rachis; the florets bear six white tepals and six staminodes (Fig. 1d–f). This uniformity indicates that the structure of reproductive organs is relatively conservative at the intraspecific level, whereas the principal differences among the three taxa lie in the vegetative organs, as detailed below.

Figure 1. 

(a), (d), (g), (j) Morphological comparison of D. bulbifera var. bulbifera; (b), (e), (h), (k) D. bulbifera var. kacheo; and (c), (f), (i), (l) D. bulbifera f. domestica. (a)–(c) Female inflorescences; (d)–(f) female flowers; (g)–(i) leaves; (j)–(l) stems.

Leaf morphology (Fig. 1g–i) shows moderate variation. Both D. bulbifera var. bulbifera and D. bulbifera var. kacheo exhibited broadly ovate-cordate leaves, whereas D. bulbifera f. domestica showed broadly cordate leaves with a length nearly equal to the width.

Stem morphology (Fig. 1j–l) exhibited more pronounced differentiation. Dioscorea bulbifera var. kacheo was characterized by ridged stems with membranous wings along the ridges, a feature absent in the other taxa. In contrast, both D. bulbifera var. bulbifera and D. bulbifera f. domestica possessed terete stems, sometimes with longitudinal striations. The stems of D. bulbifera var. bulbifera are often dark red or green with red stripes, while the stems of D. bulbifera f. domestica are usually emerald green and subglabrous.

Bulbil morphology (Fig. 2a–c) varied substantially among taxa. Dioscorea bulbifera var. bulbifera produced subglobose bulbils with irregularly arranged depressions, whereas D. bulbifera var. kacheo bore subglobose to broadly ellipsoid bulbils with regularly arranged protuberances. In contrast, D. bulbifera f. domestica exhibited highly variable bulbil morphology, often ridged and bearing irregularly arranged protuberances of different sizes.

Figure 2. 

Morphological comparison of (a), (d) D. bulbifera var. bulbifera; (b), (e) D. bulbifera var. kacheo; and (c), (f) D. bulbifera f. domestica. (a)–(c) Bulbils. (d)–(f) Tubers. Scale bar = 5 cm.

Tuber morphology (Fig. 2d–f) differed conspicuously among the three taxa. Dioscorea bulbifera var. bulbifera formed well-developed, swollen underground tubers that occasionally reach several kg in weight. These tubers are usually subglobose to ovoid-globose, with a blackish-brown and easily detached epidermis and yellow flesh, and they give rise to several adventitious roots. The tubers of D. bulbifera var. kacheo are essentially similar in morphology and structure to those of D. bulbifera var. bulbifera, differing mainly in their somewhat smaller size (4.0–5.0 vs 8.0–13.0 cm in diameter). In marked contrast, D. bulbifera f. domestica scarcely forms tubers (1.0–1.5 cm in diameter where present); instead, the junction between the stem and rhizome bears well-developed adventitious roots. Quantitative measurements of tuber diameter for the three taxa are provided in Table 1.

Chromosome number and genome size estimates

Chromosomes were extremely small and dot-like across all samples, requiring well-spread metaphase plates for accurate counting (Fig. 3). In D. bulbifera var. bulbifera, chromosome numbers varied among the examined accessions: the Hainan Lingshui accession (HNLS) had 2n = 40 (Fig. 3a), the Hainan Haikou accession (HNHK) had 2n = 60 (Fig. 3b), and the Zhejiang Lishui accession (ZJLS) had 2n = 80 (Fig. 3c). In contrast, the examined accession of D. bulbifera var. kacheo (JCHD1) consistently showed 2n = 80 across all counted cells (Fig. 3d), and D. bulbifera f. domestica (YZ1) uniformly displayed 2n = 60 (Fig. 3e).

Figure 3. 

Chromosome spreads of D. bulbifera var. bulbifera, D. bulbifera var. kacheo and D. bulbifera f. domestica. (a) D. bulbifera var. bulbifera (HNLS, 2n = 40); (b) D. bulbifera var. bulbifera (HNHK, 2n = 60); (c) D. bulbifera var. bulbifera (ZJLS, 2n = 80); (d) D. bulbifera var. kacheo (JCHD1, 2n = 80); (e) D. bulbifera f. domestica (YZ1, 2n = 60). Scale bar = 2 μm.

To further corroborate these ploidy differences, genome sizes of representative accessions of the three taxa were estimated by flow cytometry (Table 2). The two 2n = 80 accessions (ZJLS, 1.31 Gb; JCHD1, 1.38 Gb) had the largest genome sizes, approximately twice that of the 2n = 40 accession (HNLS, 0.66 Gb), whereas the two 2n = 60 accessions (HNHK, 1.02 Gb; YZ1, 1.15 Gb) showed intermediate values, with minor variation among accessions sharing the same chromosome number. Thus, the proportional differences in genome size among accessions are consistent with the ratios of their chromosome numbers.

Plastome features and comparison

A total of 17 plastomes were assembled (Supplementary Table S2; Fig. 4), with genome sizes ranging from 152,970 to 153,256 bp. The smallest genome was observed in the D. bulbifera var. bulbifera accession HeNSC (152,970 bp), whereas the largest genome was found in D. bulbifera f. domestica accession YZ2 (153,256 bp).

Figure 4. 

Plastome map of D. bulbifera. Genes on the outer circle are transcribed clockwise, and those on the inner circle are transcribed counterclockwise. Functional categories are indicated by color-coded genes. The dark grey and the light grey plots in the inner circle correspond to the GC content and AT content, respectively.

All plastomes exhibited the typical quadripartite structure, comprising a pair of inverted repeats (IRs: 25,499–25,505 bp), a large single-copy (LSC: 83,120–83,394 bp), and a small single-copy (SSC: 18,851–18,855 bp) region (Supplementary Fig. S1). The overall GC content was identical across all samples (37.0%), regardless of geographic origin or taxonomic grouping (Supplementary Table S2). Each plastome contained 113 genes, including 79 protein-coding genes, 30 tRNA genes, and 4 rRNA genes, with 19 genes duplicated in the IR regions. Gene content and order were highly conserved across all accessions.

Comparative analysis of IR/SC junctions revealed highly consistent boundary configurations among all samples (Supplementary Fig. S1). Specifically, the LSC/IRb boundary (JLB) was located within the rps19 gene, while the LSC/IRa boundary (JLA) was positioned in the psbA-trnH intergenic region. The SSC/IRb boundary (JSB) occurred within the ycf1 gene, and the SSC/IRa boundary (JSA) was located adjacent to ndhF. Only minor length variation (≤ 6 bp) in IR regions was detected in two D. bulbifera f. domestica accessions (YZ1 and YZ2), without altering the overall boundary structure.

Phylogenetic relationships, population structure, and gene flow

We first inferred maximum-likelihood (ML) trees from the concatenated coding sequences (CDS) and amino acid sequences of 79 plastid genes (Supplementary Figs. S2, S3). In both trees, only the monophyly of D. bulbifera f. domestica received strong support, whereas relationships among the remaining ingroup accessions were poorly resolved. We therefore analysed the complete plastome alignment using both ML and Bayesian inference (BI). These whole-plastome analyses provided stronger resolution of infraspecific relationships within D. bulbifera. Because ML and BI recovered identical ingroup topologies, only the ML tree is presented (Fig. 5a), with SH-aLRT, UFBoot, and PP values annotated on the nodes.

Figure 5. 

Phylogenetic relationships and cytonuclear discordance in D. bulbifera. (a) Maximum likelihood (ML) tree based on complete plastome sequences; node labels show SH-aLRT/UFBoot/PP (SH-aLRT support [%], ultrafast bootstrap support [%], and Bayesian posterior probabilities). (b) ML tree based on genome-wide nuclear SNPs; node labels show SH-aLRT/UFBoot support (%), with * indicating SH-aLRT/UFBoot = 100/100. All trees were rooted with D. zingiberensis.

In the plastome phylogeny (Fig. 5a), the two D. bulbifera f. domestica accessions formed a strongly supported monophyletic clade that was recovered as sister to all remaining D. bulbifera samples (SH-aLRT/UFBoot/PP = 100/100/1.00). Dioscorea bulbifera f. domestica therefore constitutes an early-diverging lineage within the species. By contrast, the two D. bulbifera var. kacheo accessions failed to cluster into a single clade. Instead, they were recovered as a monophyletic clade together with D. bulbifera var. bulbifera populations from southeastern, southern, and southwestern China (Fujian, Taiwan, Guangdong, Guangxi, Chongqing) (SH-aLRT/UFBoot/PP = 100/100/1.00).

In the nuclear SNP phylogeny (Fig. 5b), D. bulbifera f. domestica also formed a distinct, well-supported clade (SH-aLRT/UFBoot = 100/100) that was sister to the remaining ingroup. Dioscorea bulbifera var. kacheo, however, was no longer nested within D. bulbifera var. bulbifera: JCHD1 and JCHD2 formed a strongly supported clade (SH-aLRT/UFBoot = 100/100), which was sister to a clade comprising all D. bulbifera var. bulbifera accessions (SH-aLRT/UFBoot = 100/100). Dioscorea bulbifera var. bulbifera itself was likewise recovered as a single well-supported clade (SH-aLRT/UFBoot = 100/100).

Thus, the two trees are concordant in placing D. bulbifera f. domestica as a distinct early-diverging lineage, but they are discordant for D. bulbifera var. kacheo: this taxon is nested within D. bulbifera var. bulbifera in the plastome tree, yet recovered as an independent sister clade to D. bulbifera var. bulbifera in the nuclear SNP tree (Fig. 5).

PCA and ADMIXTURE analyses were used to further resolve population structure. PCA separated the three groups, with PC1 (16.91%) and PC2 (13.39%) explaining 30.30% of the total genetic variation (Fig. 6a). Dioscorea bulbifera f. domestica formed a distinct cluster completely separated from D. bulbifera var. bulbifera and D. bulbifera var. kacheo, whereas D. bulbifera var. kacheo was positioned at the periphery of the HD cluster.

Figure 6. 

Population genomic analyses of three D. bulbifera infraspecific taxa: HD (D. bulbifera var. bulbifera), JC (D. bulbifera var. kacheo), and YZ (D. bulbifera f. domestica). (a) PCA based on genome-wide SNPs; PC1 and PC2 explain 16.91% and 13.39% of genetic variation. Orange, purple, and blue denote HD (n = 13), JC (n = 2), and YZ (n = 2), respectively. (b) ADMIXTURE analysis (K = 1–8) showing ancestry proportions of each sample. (c) Network of (π) nucleotide diversity and (F_ST) population differentiation; values inside circles indicate π, and lines show pairwise F_ST (white to red indicates increasing differentiation). (d) TreeMix analysis (m = 1) showing population divergence and gene flow; the yellow arrow indicates migration from YZ to JC (weight = 0.183).

ADMIXTURE analysis showed that K = 1 had the lowest cross-validation error (Supplementary Fig. S4). At K = 2, D. bulbifera f. domestica was differentiated from the D. bulbifera var. bulbifera/D. bulbifera var. kacheo groups with minimal admixture (Fig. 6b); this model was retained for its stronger biological interpretability.

Nucleotide diversity and genetic differentiation were consistent with these patterns (Fig. 6c). JC exhibited the highest nucleotide diversity (π = 6.380 × 10⁻³), followed by D. bulbifera var. bulbifera and D. bulbifera f. domestica. Pairwise F_ST values indicated strong differentiation between JC and D. bulbifera f. domestica (0.5051), intermediate differentiation between D. bulbifera var. bulbifera and D. bulbifera f. domestica (0.2801), and the lowest between D. bulbifera var. bulbifera and D. bulbifera var. kacheo (0.0752).

D-statistic analyses (Supplementary Fig. S5) yielded significantly positive values (D[HD, JC, YZ, O] > 0, |Z| > 3, p < 0.001), indicating excess allele sharing between D. bulbifera var. kacheo and D. bulbifera f. domestica. TreeMix analysis identified m = 1 as the optimal model by the Δm method and revealed a single migration edge from YZ to JC (weight = 0.183) (Fig. 6d; Supplementary Figs. S6, S7).

Collectively, these results show that D. bulbifera var. bulbifera and D. bulbifera var. kacheo are genetically closely related and only weakly differentiated, whereas D. bulbifera f. domestica is clearly separated from both groups. Significant D-statistic and TreeMix signals of excess allele sharing between D. bulbifera var. kacheo and D. bulbifera f. domestica were also detected.

Our integrative dataset revealed clear morphological and cytological differentiation among the three taxa, together with pronounced asymmetry in their levels of genomic divergence. Morphologically, both focal taxa are readily diagnosable against D. bulbifera var. bulbifera: D. bulbifera var. kacheo is characterized by ridged stems bearing membranous wings along the ridges, whereas D. bulbifera f. domestica is distinguished by broadly cordate leaves with length nearly equal to width, highly variable and often ridged bulbils, and the near absence of swollen underground tubers. Cytologically, chromosome counts and concordant flow-cytometric genome-size estimates revealed ploidy variation within D. bulbifera: the examined D. bulbifera var. bulbifera accessions spanned 2n = 40, 60, and 80, whereas the sampled accessions of D. bulbifera var. kacheo (2n = 80) and D. bulbifera f. domestica (2n = 60) were each uniform and differed from one another, indicating that cytotype variation is partitioned among as well as within infraspecific taxa.

Phylogenomic and population-genetic analyses further resolved the two taxa, but to markedly different degrees. Dioscorea bulbifera f. domestica was consistently differentiated by every class of molecular evidence: it formed a strongly supported, early-diverging lineage in the plastome tree, a maximally supported clade in the nuclear SNP tree, and a discrete cluster in PCA, carried a nearly undiluted ancestry component in ADMIXTURE, and showed substantial differentiation from D. bulbifera var. bulbifera (F_ST = 0.2801). In contrast, D. bulbifera var. kacheo was not separable from D. bulbifera var. bulbifera in several analyses: it was nested among D. bulbifera var. bulbifera accessions in the plastome tree, occupied the periphery of the var. bulbifera cluster in PCA, and showed only weak genome-wide differentiation (F_ST = 0.0752), although it was recovered as a distinct clade in the nuclear SNP tree. Concordance among independent character systems is widely used to evaluate lineage separation^[52], and integrative morphological–genomic frameworks have proven effective for refining taxon boundaries in angiosperms, including Dioscorea^[53]. On this basis, D. bulbifera f. domestica more closely approximates a genetically differentiated infraspecific lineage than does D. bulbifera var. kacheo.

The D. bulbifera var. kacheo pattern—clear morphological diagnosis without commensurate genome-wide differentiation from D. bulbifera var. bulbifera—indicates incongruence between phenotype and nuclear population structure. Pronounced phenotypic differentiation accompanied by weak genetic isolation has been reported in Dianthus virgineus^[54], and incongruence between phenotypic clusters and genomic lineages has likewise been documented in the Aquilegia viridiflora complex^[55]. Morphological diagnosability alone is therefore insufficient to equate D. bulbifera var. kacheo with D. bulbifera f. domestica in degree of evolutionary independence.

This asymmetry is at odds with historical floristic treatments that assigned varietal rank to D. bulbifera var. kacheo and forma rank to D. bulbifera f. domestica; those assignments reflect regional morphological conventions rather than relative levels of genome-wide differentiation. In botanical practice, varieties are generally circumscribed as morphologically diagnosable entities that additionally possess geographic, ecological, and/or phylogenetic integrity, whereas formae are applied to morphological variants lacking such extramorphological integrity^[14]. Following this conventional distinction, and integrative approaches emphasizing cross-validation among evidence classes^[52,53], we applied operational criteria under which varietal status requires simultaneous support from diagnosable morphology, substantial genome-wide differentiation, and concordance across phylogenomic and population-genetic analyses; morphological divergence without accompanying genome-wide differentiation is regarded as more consistent with forma^[54,55]. Under these criteria, D. bulbifera f. domestica more closely meets the evidentiary expectations of varietal rank, whereas D. bulbifera var. kacheo satisfies the morphological criterion primarily. Given that each focal taxon was represented by only two accessions from narrow ranges, and that the observed cytological uniformity may not hold across their full distributions, formal nomenclatural rearrangement would be premature. We therefore retain the historically applied ranks and treat the ranking implied by our data as a testable working hypothesis; expanded sampling of D. bulbifera var. kacheo in Bhutan and adjacent regions and of wild D. bulbifera f. domestica populations in Japan will be required before formal revision can be justified.

Cytonuclear discordance and reticulate evolution

The plastome and nuclear SNP trees were concordant in recovering D. bulbifera f. domestica as a distinct, early-diverging lineage, but they were strongly discordant in the placement of D. bulbifera var. kacheo. In the plastome phylogeny, the two D. bulbifera var. kacheo accessions were nested among D. bulbifera var. bulbifera samples rather than forming a sister lineage to that taxon. In the nuclear SNP tree, by contrast, D. bulbifera var. kacheo formed an independent, well-supported clade that was sister to all D. bulbifera var. bulbifera accessions. Plastomes reflect maternal inheritance, whereas nuclear genomes capture biparental evolutionary processes^[27]. Cytonuclear discordance is widespread in plants and is typically attributed to incomplete lineage sorting (ILS), hybridization, and introgression^[27,56]. These processes may act concurrently, particularly in recently diverged or polyploid lineages.

In this study, population genomic analyses provide complementary evidence for historical gene exchange. Significantly positive D-statistics indicate excess allele sharing between D. bulbifera var. kacheo and D. bulbifera f. domestica relative to the tested topology^[57], and TreeMix recovered a migration edge from D. bulbifera f. domestica to D. bulbifera var. kacheo^[49]. These results are consistent with introgression contributing to the observed discordance, but they do not exclude a substantial contribution from ILS. Formal discrimination between introgression and retained ancestral polymorphism remains challenging under the multispecies coalescent, especially with limited per-taxon sampling^[56,58]. We therefore interpret the phylogenetic incongruence as more consistent with a reticulate evolutionary model than with a strictly bifurcating history, while acknowledging that the relative contributions of introgression and ILS cannot be fully partitioned with the present dataset^[56,59]. These patterns are compatible with a reticulate history involving episodic gene flow; given the cytological heterogeneity documented in D. bulbifera, interploidy gene flow may also be relevant^[18,59,60]. The congruent placement of D. bulbifera f. domestica in both genomic compartments further indicates that its distinctiveness is not confined to a single inheritance system.

Several limitations should be acknowledged. Only two accessions each of D. bulbifera var. kacheo and D. bulbifera f. domestica were available. Consequently, population genetic metrics based on these groups (π, F_ST, D-statistic, and TreeMix migration weights) should be regarded as exploratory: absolute values may be sensitive to sampling, even though the direction of the patterns was consistent across phylogenies, PCA, ADMIXTURE, differentiation statistics, and gene-flow analyses. For this reason, and consistent with the discussion above, taxonomic conclusions regarding the relative infraspecific rank of D. bulbifera var. kacheo remain tentative. The clear separation of D. bulbifera f. domestica across all analyses indicates that its differentiation is unlikely to be solely an artifact of sampling.

Nevertheless, broader sampling across geographic ranges will be essential for further resolving fine-scale population structure and evaluating the full extent of genetic diversity within each lineage. Future studies integrating larger population datasets and demographic modeling will help to refine the evolutionary history and taxonomic boundaries of D. bulbifera.

This study provides an integrative framework for evaluating infraspecific taxonomic relationships within D. bulbifera by combining genome-wide nuclear SNP data, plastome variation, cytological evidence, and morphological traits. Our results indicate that D. bulbifera f. domestica represents a genetically differentiated infraspecific lineage, supported by concordant patterns across phylogenetic and population-genomic analyses. In contrast, D. bulbifera var. kacheo is morphologically diagnosable but shows only limited genome-wide differentiation from D. bulbifera var. bulbifera, indicating that phenotypic divergence is not necessarily accompanied by substantial genomic differentiation.

The observed cytonuclear discordance further suggests that the evolutionary history of D. bulbifera is unlikely to conform to a strictly bifurcating model. Population-genomic analyses are consistent with historical gene flow contributing to this incongruence; however, incomplete lineage sorting cannot be excluded, and the relative contributions of these processes—and of potential interploidy gene flow—cannot yet be fully partitioned.

Overall, these findings clarify the asymmetric differentiation between historically recognized infraspecific taxa of D. bulbifera and underscore the value of integrating morphological, cytological, and genomic evidence when evaluating infraspecific boundaries in taxonomically complex plant groups. Formal nomenclatural revision, however, should await broader geographic sampling.