Integrative phylogenomics resurrects Unio molleuri (Bivalvia: Unionidae) as a valid species from the transboundary Lancang–Mekong basin

Transboundary biodiversity hotspots present unique challenges for evolutionary research, as political and logistical barriers have long fragmented biological exploration. The Indo-Burma Biodiversity Hotspot, anchored by the Lancang–Mekong River system, exemplifies this dilemma. Southern Yunnan Province, China, shares a 4,060-km international border with Vietnam, Laos, and Myanmar, encompassing rugged mountains and deeply incised valleys that have historically impeded comprehensive field surveys^[1]. More critically, national boundaries have constrained collaborative sampling and integrative taxonomic comparisons, leaving the evolutionary affinities of many cross-border taxa poorly understood. The Lancang–Mekong basin supports exceptional freshwater diversity, ranking among the top three global centers for freshwater biodiversity and harboring 899 fish species alongside rich but understudied invertebrate faunas^[2,3]. The upper reaches in southern Yunnan represent a critical biogeographic transition zone, yet the lack of cross-border phylogenetic integration continues to obscure lineage diversification and faunal connectivity across the basin.

Freshwater mussels (Bivalvia: Unionida) in this region epitomize the intersection of taxonomic uncertainty and limited genomic sampling. Widespread shell plasticity and morphological convergence, compounded by historically fragmented molecular data, have resulted in pervasive nomenclatural instability and unresolved higher-level relationships^[4,5]. The Unionetta fabagina (Deshayes, 1876) species complex illustrates these difficulties. Originally described as three distinct nominal taxa, namely Unio broti Deshayes, 1876, Unio molleuri Morlet, 1891, and Unio fabagina Deshayes in Deshayes & Jullien, 1876, they have undergone a complex nomenclatural history marked by shifting generic placements and morphological synonymy. Early systematic treatments by Simpson^[6,7] classified all three within Parreysia. However, Haas^[8] subsequently synonymized Unio molleuri under 'Unionella' fabagina, using a generic name that is a junior homonym of Unionella Etheridge, 1888, and that he later replaced with Unionetta F. Haas, 1955, thereby establishing Unionetta fabagina (Deshayes, 1876)^[9]. This taxonomic framework was subsequently adopted by Brandt^[10] and formally consolidated in the integrative molecular revision by Pfeiffer et al.^[11] and Bolotov et al.^[12]. However, phylogenetic analysis in Pfeiffer et al.^[11] recovered a highly divergent ?Unionetta fabagina specimen from the Tonle Sekong watershed that rendered the genus Unionetta non-monophyletic. This finding exposed the limitations of morphology-based synonymy and indicated that the Unionetta fabagina species complex harbors unrecognized higher-level diversity.

Subsequently, specimens from Xishuangbanna, China were identified as Unionetta fabagina^[13]. More recently, Chen et al.^[14] acknowledged the non-monophyly of Unionetta and elevated the Xishuangbanna population to a new species, Unionetta xishuangbannaensis. Shortly thereafter, Xiang et al.^[15] established the new genus Chenunionetta to accommodate U. xishuangbannaensis, thereby aiming to preserve the monophyly of Unionetta. However, the taxonomic treatments in both studies relied exclusively on a single mitochondrial marker (COI) and shell morphology, unchecked against historical type specimens of regional synonyms like Unio molleuri. By relying on potentially variable shell characters and bypassing type-based verification, they left both the species-level validity and generic assignment of the Xishuangbanna population unresolved. Consequently, a definitive taxonomic resolution requires reconciling historical nomenclatural priority with robust phylogenetic evidence.

Phylogenomics has fundamentally transformed our ability to resolve evolutionary relationships, particularly in taxonomically complex groups where traditional markers fail to disentangle rapid radiations, morphological convergence, or historical gene flow^[16,17]. However, applying phylogenomic-scale data to non-model invertebrates remains constrained by practical limitations, including DNA quality, sequencing costs, and the scarcity of standardized nuclear markers across disparate studies^[18]. Genome skimming has emerged as a cost-effective alternative, routinely yielding complete mitochondrial genomes and high-copy nuclear regions from low-coverage sequencing^[19−21]. Recent methodological advances further demonstrate that genome skimming data can be leveraged to recover anchored hybrid enrichment (AHE) loci, enabling seamless integration with published phylogenomic datasets^[22,23]. This approach not only maximizes data utility but also provides a scalable framework for quantifying gene tree discordance and distinguishing between incomplete lineage sorting (ILS) and historical introgression, which are key evolutionary processes that frequently obscure phylogenetic inference in rapidly diversifying lineages yet remain underexplored in many freshwater invertebrate groups^[24,25].

Here, we apply an integrative phylogenomic framework to resolve the systematic status of the Xishuangbanna lineage. By combining comparative shell morphology with direct type-specimen comparisons, a newly assembled F-type complete mitogenome, and AHE nuclear loci captured in silico from genome skimming data, we aim to: (1) validate the taxonomic identity of this lineage against historical type specimens; (2) determine its species validity and resolve its higher-level phylogenetic placement within Parreysiinae; and (3) quantify gene tree discordance to infer the relative contributions of ILS and introgression to its diversification. Beyond resolving a long-standing taxonomic puzzle, this study demonstrates the practical utility of extracting standardized nuclear markers from genome skimming data, offering a transferable pipeline for unionid systematics. Ultimately, by bridging cross-border knowledge gaps and applying phylogenomic rigor to a representative Lancang–Mekong taxon, this work advances our understanding of freshwater mussel evolution in the Indo-Burma hotspot and provides a robust foundation for evidence-based conservation in this threatened transboundary river system.

Freshwater mussels were collected from the Luosuo River, a first-order tributary of the Lancang–Mekong River system in Menglun Town, Xishuangbanna Dai Autonomous Prefecture, Yunnan Province, China, in close proximity to the type locality (Muang Sing, Laos) of Unio molleuri (Fig. 1). A total of 10 live individuals were hand-collected and used for molecular sequencing. No endangered or protected species were sampled. Soft tissues for molecular analysis were preserved in anhydrous ethanol at –20 °C. Shells for morphological examination were air-dried, labeled, and deposited in the Museum of Biology, Nanchang University (NCUMB), China (voucher codes: 21_NCU_XPWU_UF01–21_NCU_XPWU_UF10).

Figure 1. 

Sampling locality and species delimitation of Chenunionetta molleuri comb. nov. (a), (b) Map and general view of the sampling locality. The map was generated using ArcGIS (v10.2, Esri, Redlands, CA, USA). (c) Neighbor-joining tree of 224 COI sequences (Supplementary Fig. S1), based on Kimura 2-parameter model with 1,000 bootstrap replicates. Bootstrap support values ≥ 50% are shown at the nodes. Tip labels of Chenunionetta lineage are highlighted in red. Results from three delimitation analyses are shown; colored bars indicate inferred groupings and allow visual comparison of congruence among methods.

Shell morphology was examined using a caliper with an accuracy of ±0.1 mm to measure maximum length (L), height (H), and width (W). Comparative morphological assessments utilized type specimens of Unio fabagina Deshayes, 1876 (syntype: MNHN-IM-2000-1694), Unio broti Deshayes 1876 (syntype: MNHN-IM-2000-1652), and Unio molleuri Morlet, 1891 (syntype MNHN-IM-2000-1751; paratype SMF 3599)^[26] (Fig. 2). Shell shape, umbo position and sculpture, hinge dentition, and nacre coloration were documented using standardized photography.

Figure 2. 

Shells of Unionetta fabagina species complex. (a) Unionetta fabagina (Unio fabagina Deshayes, 1876, syntype MNHN-IM-2000-1694). (b), (c) Sequenced specimens of Unionetta fabagina in Pfeiffer et al. (2018)^[11]. (d) Unionetta fabagina (Unio broti Deshayes, 1876, syntype MNHN-IM-2000-1652). (e) ?Unionetta fabagina in Pfeiffer et al. (2018)^[11]. (f)–(h) Unio molleuri Morlet, 1891 (Syntype MNHN-IM-2000-1751, paratype SMF 3599, and NCUMB 21_NCU_XPWU_UF10, respectively). Scale bar in (h) = 1 cm and applies only to panel (h); remaining panels are scaled proportionally to reflect approximate size differences among specimens.

Institutional abbreviations: MNHN, Muséum National d'Histoire Naturelle, Paris; SMF, Senckenberg Forschungsinstitut und Naturmuseum; UF, Florida Museum of Natural History; NCUMB, Museum of Biology, Nanchang University.

DNA extraction, PCR, and Sanger sequencing

Genomic DNA was extracted from foot muscle tissue using the TIANamp Marine Animals DNA Kit (Tiangen Biotech, China). Three conventional gene fragments were amplified for initial phylogenetic assessment: F-type mitochondrial cytochrome c oxidase subunit I (COI) with primers LCO22me2/HCO700dy, F-type mitochondrial 16S ribosomal RNA (16S rRNA) gene with primers 16Sar-L-myt/16Sbr-H-myt, and nuclear 28S ribosomal RNA (28S rRNA) gene with primers D23F/D4RB (Huang et al.^[20] and references therein). PCR amplification was implemented in 25 μL volume containing 12.5 μL 2× Taq Plus Master MixII (Vazyme, China), 1 μL of each primer (10 μM), 1 μL template DNA, and 9.5 μL ddH2O, using the following thermal cycling conditions: initial denaturation at 94 °C for 3 min; 35 cycles of 94 °C for 30 s, 50–52 °C for 30 s, and 72 °C for 1 min; final extension at 72 °C for 5 min. Purified PCR products were bidirectionally sequenced on an ABI 3730xl DNA Analyzer (Applied Biosystems). New sequences of the unique haplotypes were deposited in GenBank under accession numbers PZ448460, PZ448462, and PZ448463.

Genome skimming and mitogenome assembly

For one representative specimen (21_NCU_XPWU_UF10), total genomic DNA was sheared to 350-bp fragments and used to construct a paired-end library. The library was sequenced on the DNBSEQ-T7 (BGI, Shenzhen, China) platform (2 × 150 bp paired-end reads) at Novogene (Tianjin, China), yielding 19.11 Gb of raw data. Raw reads were first processed with fastp v0.19.7^[27] to generate basic quality statistics and to remove unreliable data. The filtering criteria were: (1) read pairs containing adapter contamination in either read were discarded; (2) read pairs in which more than 10% of the bases in either read were uncertain (N bases) were discarded; (3) read pairs in which more than 50% of the bases in either read were of low quality (Phred quality < 5) were discarded. The resulting clean reads were then subjected to quality trimming (quality limit = 0.001; minimum number of nucleotides in reads = 50) and de novo assembled using CLC Genomic Workbench 12.0 (Qiagen, Germany). Contigs identified as mitogenome sequences were concatenated and manually inspected for overlap at both the beginning and end, resulting in the circularization of the mitogenome. Mitogenome annotation was performed using the GeSeq web server^[28], applying the invertebrate mitochondrial genetic code (translation table 5). This annotation was further validated manually for protein-coding genes (PCGs) using Geneious v.11^[29]. Additionally, PCGs and ribosomal RNA (rRNA) genes were confirmed through BLAST searches of the NCBI database and alignment with homologous genes from other published Unionidae mitogenomes. Transfer RNA (tRNA) genes were identified using ARWEN v.1.2^[30]. The complete mitogenome was deposited in GenBank under accession number PZ272566. To analyze strand asymmetry, GC skew and AT skew were calculated using the formulas: GC skew = (G − C)/(G + C) and AT skew = (A − T)/(A + T).

Recovery of AHE loci from genome skimming data

To maximize the phylogenomic utility of genome skimming reads, we recovered Unioverse (the freshwater mussel-specific anchored hybrid enrichment probe set^[31]) using GeneMiner2^[23]. Of the 811 targeted loci originally designed, 569 loci of Unioverse had a gene occupancy of at least 70% across Unionidae^[31] and were used as the reference target set in this study. Quality-trimmed genome skimming reads were processed through the GeneMiner2 pipeline with adaptive k-mer selection enabled and fine-grained read selection with strand-orientation verification. Contigs were extracted and aligned per locus using MAFFT v7.52^[32]. Alignments were trimmed with trimAl v1.2^[33].

Loci recovered in the target specimen were merged with published AHE data for 10 Parreysiinae and two Nodularia (subfamily Unioninae) outgroup taxa^[31] to construct a concatenated supermatrix. Published AHE datasets were downloaded from NCBI Sequence Read Archive (SRA) under accession numbers including SRR8473063 (Harmandia somboriensis), SRR8473072 (Scabies phaselus), SRR8473032 (Indochinella pugio), SRR8473076 (Indonaia caerulea), SRR8473042 (Coelatura choziensis), SRR8473053 (Prisodontopsis aviculaeformis), SRR8473030 (Leoparreysia tavoyensis), SRR8473038 (Leoparreysia olivacea), SRR8473060 (Lamellidens corrianus), SRR8473059 (Lamellidens generosus), SRR8473033 (Nodularia douglasiae), and SRR8473031 (Nodularia jourdyi).

Phylogenetic analyses

Conventional marker datasets

Two datasets were constructed: (i) the COI dataset (224 sequences; 648 bp) (Supplementary Table S1); (ii) the three-gene dataset (COI + 16S + 28S; 66 sequences; 1,866 bp) (Supplementary Table S2).

COI was codon-aligned by MUSCLE^[34], implemented in MEGA X^[35]. MAFFT v7.490 with the Q-INS-i algorithm was used for the multiple sequence alignments of ribosomal genes (16S and 28S rRNA genes). Divergent and ambiguously aligned blocks were removed from each gene alignment using Gblocks v0.91b^[36] integrated in PhyloSuite v1.2.3^[37]. Otherwise, the default settings were applied. Subsequently, trimmed alignments were concatenated in PhyloSuite for downstream analyses. The absent sequences (16S and/or 28S) of some taxa were treated as missing data.

MEGA X was used to perform neighbor-joining (NJ) analysis and calculate inter- and intraspecific distances based on the COI dataset using K2P (Kimura 2-parameter) distances. We used three distinct approaches to determine the number of molecular operational taxonomic units (MOTUs). These approaches were the Poisson Tree Process model (bPTP)^[38], haplotype network reconstructions in TCS v1.21^[39], and the Automatic Barcode Gap Discovery (ABGD) method^[40] implemented via the SPART web portal (https://spartexplorer.mnhn.fr/delimitation) using the K80 (Kimura 2-parameter) distance model with default settings (P_min = 0.001, P_max = 0.1, Steps = 10, NbBins = 20).

The best-fit partitioning schemes and substitution models were selected by PartitionFinder2 v2.3.4^[41] under greedy search with unlinked branch lengths. Predefined data blocks for partitioning scheme searches were designated by gene region (rRNA gene) or codon position (PCG). The optimal model selection and partitioning schemes were determined using the AICc method (Supplementary Table S3). RAxML v8.2.12^[42] was used for ML analyses by using the GTRGAMMA model with 1,000 rapid bootstrap replicates. Maximum parsimony (MP) was performed in PAUP* v4.0a169^[43] using heuristic searches with 100 random sequence addition replicates, conducted with tree-bisection-reconnection branch-swapping and 1,000 bootstrap replicates. MrBayes v3.2.7a^[44] was used for Bayesian inference (BI) analyses by four independent Markov chain Monte Carlo (MCMC) runs with 10 million generations. Trees were sampled every 1,000 generations with a 25% burn-in. Sufficient mixing of the chains was considered to be reached when the average standard deviation of split frequencies was below 0.01, and the potential scale reduction factor (PSRF) was 1.

Mitogenome phylogeny

The complete mitogenome of Chenunionetta molleuri was aligned with 92 published mitogenomes retrieved from GenBank (Supplementary Table S4), representing all major subfamilies of Unionidae (Unioninae, Gonideinae, Ambleminae, and Parreysiinae) and Margaritiferidae (Margaritiferinae and Gibbosulinae). The nucleotide dataset consisted of 12 concatenated protein-coding genes (PCGs) and two ribosomal RNA (rRNA) genes; we excluded atp8 due to high sequence variability. The mitochondrial phylogenomic analyses employed an analytical framework consistent with the three-gene approach described above, including codon-aware alignment by MUSCLE, ambiguous-site removal with Gblocks, optimal partitioning and model selection via PartitionFinder2, and tree inference using RAxML (ML), MrBayes (BI), and PAUP (MP).

AHE nuclear phylogeny

The concatenated AHE supermatrix was analyzed under both concatenation and coalescent-based frameworks. Phylogenetic trees were inferred using both IQ-TREE v2.2.0^[45] and FastTree 2.1^[46] based on the concatenated trimmed alignments, employing 1,000 bootstrap replicates for branch support estimation. Additionally, a coalescent-based species tree was constructed using ASTRAL-III^[47]. Single-gene trees were inferred from the individual gene alignment datasets in FastTree under the GTR + G model, and branch support was assessed using 1,000 bootstrap replicates. Quartet support values were calculated using ASTRAL-III to assess the proportion of quartets from gene trees that support each internal branch.

Gene tree discordance analysis

To quantify the relative contributions of incomplete lineage sorting (ILS) and introgression to gene tree heterogeneity, we analyzed quartet topologies of the coalescent-based species tree using Phytop v1.0^[24]. Phytop defines ILS and introgression/hybridization (IH) indices based on the proportions of quartet trees that support alternative topologies among triplets of taxa. For each relevant triplet involving Chenunionetta, Harmandia, and their closest relatives, we computed ILS-i and IH-i (the calculated ILS and IH indices) as well as ILS-e and IH-e (the proportions of gene tree topological incongruence explained by ILS and IH, respectively). The significance of quartet asymmetry was assessed using a χ² test of whether the two alternative quartet topologies (q₂ and q₃) occur at equal frequencies.

At the species level, the COI barcode dataset recovered Chenunionetta as a distinct clade, separated from all other Indochinellini by deep genetic divergences (Fig. 1c). Species delimitation analyses of this dataset delineated 15 and 16 putative species using the TCS and bPTP methods, respectively. The ABGD analysis, based on the K80 distance model, consistently recovered 11 putative species across all recursive partitions. However, we integrated these results with pairwise genetic distances and obtained an estimate of 12 distinct species. The minor discrepancy between the ABGD result (11 species) and our final estimate (12 species) arises from the treatment of three Scabies taxa. In the ABGD analysis, Scabies scobinatus, S. phaselus, and S. anceps were merged into a single MOTU. Nevertheless, the K2P genetic distances between S. scobinatus and S. phaselus (2.72% ± 0.56%) and between S. scobinatus and S. anceps (2.69% ± 0.56%) approached the commonly adopted 3% threshold for congeneric species delimitation in Unionida (Table 1), supporting the retention of S. scobinatus as a distinct species. In contrast, the K2P genetic distance between S. phaselus and S. anceps was only 1.94% ± 0.52%. We therefore provisionally treat these two taxa as a single unit, as validating their specific status requires further morphological evidence that is beyond the scope of the present study. In contrast, genetic distances demonstrated substantial divergence between Chenunionetta and its closest relatives: 7.30% ± 1.10% against Unionetta fabagina, 6.41% ± 1.05% against Harmandia somboriensis, and 7.09% ± 1.14% against ?Unionetta fabagina from the Tonle Sekong watershed (Table 1). These values markedly exceed the species-level threshold typically diagnostic of congeneric divergence in Unionidae, strongly supporting the recognition of Chenunionetta molleuri as a valid species.

Phylogenetic reconstruction based on the concatenated three-gene dataset (COI + 16S + 28S) recovered Chenunionetta as a lineage within the Parreysiinae (Fig. 3). Rather than clustering with Unionetta fabagina or Harmandia somboriensis, it formed a sister relationship with ?Unionetta fabagina, albeit with low statistical support (BS < 50%) (Fig. 3; Supplementary Figs S2 and S3). Furthermore, the maximum parsimony (MP) tree also failed to resolve the relationships among Scabiellus songkramensis, ?Unionetta fabagina, Unionetta fabagina, Harmandia somboriensis, and Chenunionetta molleuri (Supplementary Fig. S4).

Figure 3. 

Maximum likelihood phylogenetic tree based on the concatenated three-gene dataset (COI + 16S + 28S; 1,866 bp). Tip label of Chenunionetta lineage is highlighted in red. Nodal support values are shown as ML bootstrap/BI posterior probability, with dashes indicating nodes that conflict between the two inference methods. Bootstrap support values ≥ 50% or posterior probability ≥ 0.5 are shown at the nodes.

Mitogenome phylogeny

The complete F-type mitochondrial genome of Chenunionetta molleuri (specimen 21_NCU_XPWU_UF10) was assembled as a circular molecule of 15,721 bp (GenBank: PZ272566). As a member of the tribe Indochinellini, C. molleuri represents the second completely sequenced mitogenome for the subfamily Parreysiinae and is currently the only known representative of this subfamily in China. It contains the typical set of 37 mitochondrial genes: 13 protein-coding genes (PCGs), 22 transfer RNA (tRNA) genes, and two ribosomal RNA genes (rrnS and rrnL) (Fig. 4a). The overall base composition was A = 25.5%, T = 39.4%, G = 23.6%, and C = 11.5%, resulting in an A+T content of 64.9%. The gene order was identical to that of Lamellidens marginalis (Parreysiinae: Lamellidentini), confirming the absence of gene rearrangements within this subfamily. Strand asymmetry was evaluated using AT and GC skews, which were calculated as −0.21 and 0.34, respectively.

Figure 4. 

Analysis based on mitochondrial data. (a) Gene map of the F-type mitochondrial genome of Chenunionetta molleuri comb. nov. (b) Maximum likelihood phylogeny based on the F-type complete mitochondrial genome dataset (12 PCGs and 2 rRNAs). Tip label of Chenunionetta lineage is highlighted in red. Nodal support values (bootstrap from RAxML/posterior probability from MrBayes) are annotated on branches, with dashes indicating nodes that conflict between the two inference methods.

Maximum likelihood (ML) and Bayesian inference (BI) analyses of the mitogenome dataset yielded congruent topologies, with most nodes receiving maximum statistical support (BS = 100%, PP = 1.00) (Fig. 4b; Supplementary Figs S5 and S6). Maximum parsimony (MP) analysis recovered a largely congruent topology, with most nodes receiving high bootstrap support (Supplementary Fig. S7). All trees strongly supported the monophyly of all recognized subfamilies, including Parreysiinae, to which C. molleuri belongs. Within the family Unionidae, the ML and BI trees recovered Parreysiinae and Unioninae as a sister group; this combined clade was in turn sister to Gonideinae, while Ambleminae formed the sister group to the remaining subfamilies.

AHE nuclear phylogeny and gene tree discordance

GeneMiner2 successfully recovered 569 AHE nuclear loci from genome skimming reads of Chenunionetta molleuri. The raw recovered sequences had a mean length of 297.3 bp (median = 289 bp; range = 79–998 bp), which was reduced to a mean of 167.5 bp (median = 147 bp; range = 0–822 bp) after trimming (Supplementary Table S5). Following consistent in silico processing pipelines, 569 target nuclear loci were extracted from published AHE datasets for 10 parreysiine taxa and two Nodularia outgroups^[31] and combined with those recovered from the genome skimming reads of C. molleuri. The resulting concatenated supermatrix for the 13 taxa comprised 569 loci (95,518 bp), including 12,812 variable sites and 9,556 parsimony-informative sites.

To evaluate the reliability of the in silico recovered data, we compared sequence statistics across all 13 samples, including 12 samples generated by targeted AHE sequencing and C. molleuri recovered from genome skimming. Gene recovery rate was uniformly high across all samples (Fig. 5a). Result and trimmed sequence lengths both followed unimodal distributions (Fig. 5b, c), with trimming reducing sequence length in a systematic manner (Fig. 5d). Read counts were broadly comparable across samples on a log₁₀ scale (Fig. 5e), and result sequence length scaled positively with read depth (Fig. 5f). The mean result sequence length per sample was generally consistent across the 12 targeted-sequencing samples; the slightly shorter mean length observed for sample C. molleuri (Fig. 5g) is consistent with the lower and more stochastic coverage of genome skimming data compared to targeted enrichment. Recovered sequences tracked the reference locus set close to the 1:1 expectation (Fig. 5h), and the majority of AHE loci exhibited low length variability across samples (median CV < 0.15; Fig. 5i), confirming the suitability of these markers for phylogenomic analysis.

Figure 5. 

Summary statistics of in silico AHE gene capture across 13 Unionidae samples. (a) Gene recovery rate for each sample. (b) Frequency distribution of result sequence lengths with mean (red dashed line) and median (orange dotted line) indicated. (c) Frequency distribution of trimmed sequence lengths with mean and median indicated. (d) Length comparison between result and trimmed sequences. (e) Distribution of read counts on a log₁₀ scale. (f) Scatter plot showing the relationship between result sequence length and read depth (log₁₀ scale). (g) Mean result sequence length per sample. (h) Comparison of result sequence length with reference median length; the dashed line indicates the 1:1 relationship. (i) Distribution of the coefficient of variation (CV) of result sequence length across the 569 AHE loci.

Both IQ-TREE, FastTree and ASTRAL-III analyses of the AHE nuclear dataset recovered Chenunionetta as sister to Harmandia, with Scabies as sister to the (Chenunionetta+ Harmandia) clade. Quartet-based discordance analysis using Phytop revealed substantial gene tree conflict in node 1 (Fig. 6). Phytop analysis assigned an ILS index (ILS-i) of 62.0% and an IH index (IH-i) of 22.1% for the Chenunionetta–Harmandia relationship, with estimated proportions ILS-e = 41.3% and IH-e = 8.4%. The χ² test for quartet symmetry was significant (p = 0.039; significance assessed at p < 0.05), indicating that the two alternative quartet topologies are not distributed as expected under a strictly ILS-only model, and suggesting a secondary contribution from historical gene flow.

Figure 6. 

ASTRAL-III species tree inferred from 569 in silico-captured AHE nuclear loci. Branch support values (ASTRAL-III local posterior probabilities, IQ-TREE ultrafast bootstrap, and FastTree local support values) were uniformly maximal (100%) at all nodes and are therefore omitted for clarity. Insets display the quartet-based discordance analyses at the corresponding numbered nodes (1–10). Tip label of Chenunionetta is highlighted in red.

Systematic descriptions

Family Unionidae Rafinesque, 1820

Subfamily Parreysiinae Henderson, 1935

Tribe Indochinellini Bolotov, Pfeiffer, Vikhrev, & Konopleva, 2018

Genus Chenunionetta Xiang, He & Guo, 2026

Type species. Chenunionetta molleuri (Morlet, 1891) comb. nov.

Unio molleuri Morlet, 1891: 238, 242–243, pl. 7, fig. 4. [Basionym]

Chenunionetta xishuangbannaensis (Z.-G. Chen, Y.-T. Dai, H. Zheng, S. Oyang & X.-P. Wu, 2026): 57, fig. 2. [syn. nov.]

Zoobank: urn:lsid:zoobank.org:pub:DA019407-2B56-4CA6-9854-4B6400326A2B

Type locality. 'Muonngsun' (Muang Sing), Laos (sensu Haas, 1969^[9]); originally stated as 'Vallée du Mékong, sur les bords du fleuve de ce nom' by Morlet (1891)^[48].

Type material. Syntype: MNHN-IM-2000-1751 (Muséum National d'Histoire Naturelle, Paris); paratype: SMF 3599 (Senckenberg Museum, Frankfurt).

New material examined. CHINA: Yunnan Province: Xishuangbanna Dai Autonomous Prefecture: Luosuo River, 10 specimens, deposited in NCUMB (voucher codes: 21_NCU_XPWU_UF01–UF10). Shell small to medium-sized (length 36.2–41.7 mm, mean 37.8 mm; height 24.1–26.5 mm, mean 25.3 mm; width 15.8–19.9 mm, mean 18.3 mm) (Supplementary Table S6).

Distribution and habitat. Currently known only from the upper Lancang–Mekong River system, restricted to the Luosuo River (a first-order tributary near Menglun Town, Xishuangbanna, China) and historically reported from the Mekong valley (Muang Sing) in Laos. Inhabits shallow, slow-flowing reaches with substrates of mud and pebbles.

Taxonomic remarks. Historically, both Unio molleuri Morlet, 1891 and Unio broti Deshayes, 1876 were subsumed as junior synonyms of Unionetta fabagina (Deshayes, 1876) by subsequent authors (e.g. Haas, 1914^[8]; Haas, 1969^[9]; Brandt, 1974^[10]; Pfeiffer et al., 2018^[11]). This prolonged synonymy obscured the distinct validity of molleuri. Our integrative morphological and molecular analyses unequivocally demonstrate that molleuri represents a distinct evolutionary lineage, warranting its resurrection as a valid species. Direct comparison with the syntype (MNHN-IM-2000-1751) and paratype (SMF 3599) of U. molleuri confirms that the Xishuangbanna specimens match Morlet's original concept. Among our 10 voucher specimens (Supplementary Figs S8 and S9), several individuals share the shorter, more truncated posterior end comparable to the syntype, demonstrating that the observed shell shape variation falls within the intraspecific range. The recently described Chenunionetta xishuangbannaensis, which was based on this same distinct lineage, is here relegated to a junior synonym of C. molleuri (syn. nov.).

The resurrection of Unio molleuri Morlet, 1891 from synonymy under U. fabagina is validated by direct comparison with type specimens (MNHN-IM-2000-1751, SMF 3599). The genus Chenunionetta, recently established for this lineage^[15], together with the integrative evidence presented here, resolves a historical taxonomic puzzle within the Unionetta fabagina species complex. Notably, the phylogenetic topology derived from our three-gene dataset is highly consistent with the evolutionary relationships reconstructed by Pfeiffer et al.^[11] and Bolotov et al.^[12]. Although our phylogenetic analyses recovered a sister relationship between Chenunionetta molleuri and the specimen tentatively identified as ?Unionetta fabagina from the Tonle Sekong watershed, Cambodia, in Pfeiffer et al.^[11], this node received low statistical support. This lack of resolution is likely attributable to missing 16S sequences for several of these taxa, as well as the limited number of phylogenetically informative sites within the three-gene dataset for resolving deep intergeneric relationships. We refrain from addressing the formal taxonomic status of this lineage in the present study. Resolving its generic assignment will require the examination of additional specimens and the generation of comprehensive molecular data, such as complete mitochondrial genomes, in future investigations.

The recent description of Unionetta xishuangbannaensis by Chen et al.^[14], based on the same lineage, is therefore relegated to a junior synonym of Chenunionetta molleuri under the principle of nomenclatural priority established by the International Commission on Zoological Nomenclature (ICZN). While molecular systematics has undoubtedly propelled the discovery of cryptic species and refined morphology-based taxonomy, we emphasize that a new clade on a phylogenetic tree does not automatically equate to a new species. Such a genetically distinct lineage may merely represent a known species that has not yet been sequenced. Therefore, the formal description of a new species needs a comprehensive review of existing taxonomic literature, including historically synonymized names, alongside rigorous morphological comparisons with type specimens. Integrating phylogenomic data with classical taxonomy remains essential to avoid the redundant description of already established taxa.

Genome skimming for AHE nuclear loci recovery and gene tree discordance

This study demonstrates the practical feasibility of recovering anchored hybrid enrichment loci from genome skimming data using GeneMiner2, offering an economical pipeline for unionid phylogenomics. Traditional library preparation requires dedicated hybridization capture with custom probe sets, which is labor-intensive, and can be prohibitively expensive for rare taxa^[31,49,50]. By contrast, genome skimming generates shallow whole-genome coverage at a fraction of the cost, allowing off-target reads to be mined subsequently^[22,23,25].

We recovered all 569 target loci from the Unioverse probe set from a single genome skimming library, achieving recovery rates comparable to dedicated capture in unionids^[25,31]. The ability to integrate newly recovered loci with published AHE datasets is particularly valuable for systematic studies of non-model taxa. We merged our genome skimming-derived loci with the published Unioverse dataset^[31], enabling direct phylogenomic comparison across studies without re-sequencing reference taxa. This interoperability maximizes the cumulative value of phylogenomic data and reduces redundant sequencing, which is an important consideration for biodiversity research in resource-limited settings. Nevertheless, recovery rates from genome skimming are contingent upon genome size, repetitiveness, and sequencing depth. Optimizing sequencing depth and coverage evenness remains essential for maximizing locus recovery, particularly in freshwater mussels, which possess exceptionally large (2–3 Gb) and highly heterozygous genomes^[51].

A central finding of this study is the detection of substantial gene tree discordance surrounding the divergence of Chenunionetta and Harmandia. Quartet based discordance analysis revealed that incomplete lineage sorting predominates over historical introgression in shaping the observed conflict. Such heterogeneity is characteristic of rapid, ancient radiations where speciation events occur in close temporal succession, leaving insufficient time for ancestral polymorphism to sort completely^[52]. The significant quartet asymmetry nonetheless indicates a secondary contribution from historical gene flow, which is consistent with the biology of Unionidae. Freshwater mussels possess planktonic glochidial larvae that parasitize host fishes, potentially facilitating secondary contact and hybridization between divergent lineages. The complex drainage history of the Lancang and Mekong system, including river capture events and paleo-connections between previously isolated tributaries^[53−55], provides a plausible biogeographic setting for such historical gene flow.

Biogeographic implications and conservation significance

The Mekong River system acts as a significant biogeographic barrier for freshwater mussels, effectively separating the Indochinese and East Asian faunal regions. Xishuangbanna, situated at the upper terminal reach of the Lancang and Mekong basin, serves as a crucial hub and key node connecting the freshwater biotas of China and Southeast Asia. The discovery of Chenunionetta molleuri, a lineage with clear Southeast Asian affinities, in the Xishuangbanna region corroborates the pivotal role of this area in species evolution, transboundary dispersal, and ecological conservation. This genus-level endemism highlights that upper Lancang tributaries have functioned as isolated evolutionary arenas for the in-situ diversification of unionid lineages, driven by historical tectonic and climatic events such as the Tibetan Plateau uplift and the intensification of the East Asian monsoon^[20,56].

Despite its evolutionary significance, the restricted distribution of C. molleuri renders it highly vulnerable to ongoing anthropogenic pressures. Habitat degradation from sand mining, riparian deforestation, agricultural pollution, and hydropower development continuously threatens the limited range of this genus in small tributaries such as the Luosuo River. Furthermore, historical taxonomic confusion and political boundaries have long hindered comparative research on shared taxa among adjacent countries. This fragmentation of knowledge has resulted in an incomplete picture of regional biodiversity and protection needs. We strongly advocate for enhanced international communication and collaborative field surveys. Relying solely on localized populations for threat assessments without transnational data integration may lead to misdirected management efforts. Accurate preservation prioritization for this shared taxon requires a unified and systematic evaluation of its true distribution and ecological requirements across the entire Lancang and Mekong basin.

This study validates Chenunionetta molleuri (Morlet, 1891) comb. nov. as a freshwater mussel species endemic to the upper Lancang–Mekong basin. Beyond resolving a historical taxonomic puzzle, our findings validate the utility of recovering nuclear loci from genome skimming data, providing an economical approach particularly valuable for rare taxa.

Several methodological constraints point to productive avenues for subsequent investigation. First, the absence of certain key taxa from our dataset leaves open the possibility that the phylogenetic position of Chenunionetta may shift with expanded sampling. Second, the loci recovered from genome skimming data exhibited shorter mean lengths compared to anchored hybrid enrichment, which may impact the phylogenetic signal for deeper divergences within Unionidae. Third, because our interpretation of gene tree discordance relies on quartet concordance analysis, subsequent research should apply explicit coalescent modelling to quantify the relative contributions of incomplete lineage sorting and introgression. Ultimately, as political boundaries continue to fragment biodiversity exploration, integrative taxonomy anchored by type specimen verification offers a reliable path toward a comprehensive understanding of the exceptional freshwater diversity in the Indochinese region.