1. Introduction

Reproductive systems or modes are one of the most important and variable life-history traits of higher eukaryotes and are inherently linked to several fundamental population genetic parameters, including effective population size (N_e), selection efficacy, and mutation accumulation (Charlesworth and Charlesworth, 2010; Glémin and Galtier, 2012). Accordingly, reproductive systems are expected to play a key role in both species diversification and extinction (Cutter, 2019). Flowering plants provide outstanding opportunities to investigate the genetic and evolutionary consequences of variation in reproductive systems (Arunkumar et al., 2015; Zhang et al., 2022; Zeng et al., 2024). Numerous angiosperm lineages include many purely sexual species with contrasting mating systems, but also others that reproduce through both seed and clonal propagation, and some that rely largely on various mechanisms of asexual reproduction (Fryxell, 1957; Klimes et al., 1997; Silvertown, 2008; Vallejo-Marín et al., 2010; Barrett, 2015). This reproductive diversity if investigated in a phylogenetic and comparative framework enable studies on the influences of transitions from sexual to asexual reproduction on molecular evolutionary processes operating at the genome level.

The absence of recombination in plant genomes can vary widely in scale from a short segment at a hemizygous locus (Gutiérrez-Valencia et al., 2021), to large regions on the sex chromosomes (Charlesworth, 2017), to the extreme case of complete genomes in asexual organisms (Glémin and Galtier, 2012). There are two major genetic consequences of a lack of meiotic recombination that may be associated with clonal reproduction. First, the efficacy of natural selection is decreased for adaptive and deleterious mutations. Adaptive mutations at separate loci can interfere with one another in the absence of recombination thus reducing genomic rates of adaptive evolution (Hill and Robertson, 1966; Felsenstein, 1974; Keightley and Otto, 2006). Moreover, the proportion of mildly deleterious mutations is predicted to increase inexorably in clonal genomes leading to genomic degradation – Muller's Ratchet (Muller, 1964).

A second consequence of asexual reproduction is that levels of polymorphism should change at the population level (Balloux et al., 2003; Glémin and Galtier, 2012; Hartfield, 2016). In sexual populations, a unique mutation can be fixed more quickly by selection when in a homozygous state generated by segregation and recombination. However, heterozygosity (and thus polymorphism) is predicted to increase over time in asexual populations because the same mutation rarely occurs on a different copy of the chromosome and heterozygote mutations remain unrelated owing to a lack of segregation (Halkett et al., 2005; Balloux et al., 2003; Hartfield et al., 2016). Overall, the strongest theoretical prediction for the negative genomic consequences of asexuality is a decreased efficacy of selection compared to sexuality. Thus, asexuality will tend to be associated with an increased persistence and fixation of harmful mutations (Lynch et al., 1993), and this probably explains the higher extinction risk and rarity of asexual lineages compared with those that are sexual (Avise, 2008).

Investigating the contrasting genomic and evolutionary consequences of sexual versus asexual reproduction is difficult because many factors other than sex and recombination influence the effectiveness of selection, including the history and mode of asexuality, fluctuations in population size, and the potential effects of facultative sexuality (Hojsgaard and Hörandl, 2015; Hartfield, 2016; Brandt et al., 2017; Ho et al., 2019; Hodač et al., 2019; Cunha et al., 2022). One approach to testing the predictions is to conduct comparative investigations on the abundance of deleterious mutations in closely related sexual and asexual populations or lineages. Reduction in the efficacy of selection is predicted to increase the accumulation and fixation of weakly deleterious mutations. This should be evident in asexual species as a higher ratio of nonsynonymous to synonymous polymorphism (π_N/π_S) (equivalently, the ratio of 0- and 4-fold degenerate positions in protein-coding sequences, π₀/π₄) and/or divergence (d_N/d_S), representing the genomic consequences of asexuality over short- and long-time scales, respectively.

Most but not all comparative studies based on a few genes have found evidence for the reduced efficacy of selection in asexual compared to sexual populations or lineages (reviewed in Hartfield, 2016). However, mixed results have been reported even with genome-wide data, such as in Oenothera (Hollister et al., 2015), ribbon worms (Ament-Velásquez et al., 2016), hexapods (Brandt et al., 2019), and clonal fish (Kočí et al., 2020). In an extreme case a study on oribatid mites found the opposite pattern, with more efficient selection in asexual than sexual species, and large population size may have alleviated the negative consequences of mutation accumulation (Brandt et al., 2017). Given these conflicting results from previous studies, which have sometimes been limited to a small number of genes and/or taxa, it may be beneficial to conduct genome-wide analyses across numerous species in well-defined monophyletic lineages to better understand the role of reproductive mode in the evolution of genomes.

Dioscorea (Dioscoreaceae) provides a valuable comparative system for investigating the evolutionary and genomic consequences of variation in reproductive systems. The genus is most well-known because of the domestication of several clonally propagated species known as 'yams' (especially Dioscorea rotundata, D. alata, D. trifida, D. polystachya and D. esculenta), which are staple crops used as important sources of food and nutrition in tropical and subtropical regions (Sugihara et al., 2020; Bredeson et al., 2022; Li et al., 2022). All Dioscorea species are characterized by dioecy, and thus require obligate outcrossing, while they exhibit diverse reproductive systems, including sexual, asexual, and mixed reproduction. Although many species are indeed primary sexual outbreeders, a significant number also propagate to varying degrees by clonal propagation through bulbils, rhizomes, or tubers. Bulbils are a particularly effective means of clonal reproduction and are commonly formed in the axils of leaves or bracts with their size, shape and colour varying greatly among species. Bulbils function as a mechanism of clonal dispersal by gravity, water, and by foraging animals. Phylogenetic analyses and ancestral state reconstructions of Dioscorea indicate that asexual propagation by bulbil production is a derived character in most clades of the genus (Chen et al., 2022). However, a spectrum of reproductive conditions occurs in the genus (Fig. 1e–j), and these range from exclusive sexual reproduction through seeds (e.g. D. subcalva), through a mixed reproductive strategy of both seeds and bulbils (e.g. D. polystachya) to species that reproduce predominately by asexual means primarily involving only bulbils (e.g. D. melanophyma). Thus, in Dioscorea the transition from sexual to asexual reproduction with intermediate mixed sexual-asexual species allowed us to evaluate how these different reproductive modes might be associated with contrasting features of the genome and molecular evolution.

Here, we employ phylogenetic and comparative genomic approaches to investigate molecular signals associated with variation in reproductive modes of Dioscorea. Our study specifically addressed the following questions: (1) What are the patterns of divergence among reproductive modes? Specifically, we address the hypothesis that asexual and/or species with mixed reproductive modes will exhibit higher molecular evolutionary rates (d_N/d_S) than that of primarily sexual species. (2) What are the patterns of polymorphism and do nonsynonymous mutations accumulate at different rates among the reproductive modes? If so, what is the distribution patterns represented in the categories of nearly neutral (conservative) and mildly harmful (radical) mutations? To our knowledge, our study represents the first detailed examination of molecular divergence and polymorphism in a lineage of plant species with variable degrees of asexuality.

2. Materials and methods 2.1. Taxon sampling and sequencing

To evaluate the potential influence of variation in reproductive modes on genetic divergence parameters, we obtained samples or published transcriptomic data (43 and 14 samples, respectively; See Table S1) from 19 Dioscorea species comprising seven that reproduce by sexually produced seeds (D. deltoidea, D. nipponica, D. panthaica, D. zingiberensis, D. subcalva, D. esculenta, and D. composita), four in which reproduction is primarily asexual through clonally produced bulbils (D. melanophyma, D. pentaphylla, D. kamoonensis, and D. arachidna), and eight that regenerate by a mixture of both of these reproductive modes (D. alata, D. bulbifera, D. polystachya, D. fordii, D. rotundata, D. cirrhosa, D. japonica, and D. trifida). Our allocation of the 19 species into the three reproductive modes was based on our own field observations, experience with growing plants, and the literature on Dioscorea (Ting and Gilbert, 1996; Hang and Sun, 2020).

In addition, to comparing within-species polymorphism and the accumulation of deleterious mutations among plants of the three primary reproductive systems, we selected three focal species: sexual Dioscorea subcalva, asexual D. melanophyma and D. polystachya with mixed sexual and asexual reproduction for more in-depth genetic comparisons. We collected leaf tissue samples for RNA sequencing from populations of the three species, with 1–2 individuals per population (D. subcalva n = 14 individuals, D. melanophyma n = 9 individuals, D. polystachya n = 18 individuals). The samples for both divergence and within-species polymorphism involved seedlings either grown from seed or from bulbils in the glasshouse. Information on the location of samples used for divergence and polymorphism analyses are presented in Table S1 and Fig. S1, respectively.

We isolated total RNA using an RNAprep Pure Plant Kit (TIANGEN Biotechnologies Corporation, Beijing, China) and prepared the cDNA library for transcriptome sequencing using a cDNA Synthesis Kit (Illumina, San Diego, CA, USA) following the recommendations of manufacturer. The cDNA libraries were then sequenced on the Illumina HiSeq 4000 platform (Illumina, San Diego, CA, USA) to obtain PE reads (150 bp) of each cDNA. In total, we generated 57 transcriptomes with the RNA-seq data obtained from either current or previous studies (see Table S1).

2.2. Transcriptome assembly and data processing

We performed quality checks on raw reads within each data set using FastQC v.0.11.2. The clean reads were filtered from the raw sequence with Trimmomatic v.0.32 (Bolger et al., 2014) after trimming adaptors and removing both ambiguous (N > 10 %) and low-quality reads (Phred score < 30). We de novo assembled transcriptomes of each species using the Trinity v.2.1.1 (Haas et al., 2013) and discarded contigs with length less than 300 bp. We merged the similar sequence stretches by using CD-HIT v.4.8.1 (Fu et al., 2012). The longest transcript was retained per gene/graph for subsequent analysis. We evaluated the quality of assemblies by aligning reads to the assemblies using HISAT2 v.2.1.0 (Kim et al., 2015) and scoring the completeness of single-copy conserved plant orthologs using the program BUSCO (Simão et al., 2015).

Open reading frames (ORF) were identified from the filtered transcriptome using TransDecoder v.2.1 (with strand-specific option). The longest ORF was retained when multiple ORFs were recovered from a single transcript. We scanned the predicted ORFs for homology to known proteins in two ways: a BLAST v.2.3.0 (Camacho et al., 2009) search against a local version of Swiss-Prot database with hits of e-value < 1 × 10⁻⁵, and by identifying hidden Markov model protein domains with HMMER v.3.3.1 (Eddy, 2011) against the local Pfam-A.hmm protein database with domain hits of FullSeqScore and FullDomainScore > 20. We retained the final coding regions by integrating the BLAST predictions and Pfam search results. In cases where multiple ORFs were predicted for a single transcript the top score ORF per transcript was selected. Detailed statistics for the transcriptome assemblies and annotations are provided in Table S2 and Table S3, respectively.

To identify orthologous sequences shared by the 19 species for interspecific divergence analysis, we ran OrthoFinder v.2.3 (Emms and Kelly, 2019) on the identified ORF protein sequences with default parameters. Finally, 362 single-copy orthologs that were fully represented in all of the 19 species were used for further analysis. We conducted sequence alignments for each orthologous group using MAFFT v.7.4.9.0 (Katoh and Standley, 2013) and nucleotide sequences were aligned based on the corresponding protein sequences using ParaAT v.2.0 (Zhang et al., 2012). We manually checked the sequence alignment of each orthologs set.

For our comparisons of intraspecific polymorphism among reproductive modes, we mapped reads using HISAT2 v.2.1.0 (Kim et al., 2015) with default parameters of each sample back to the transcripts of three focal species, respectively. We then used a specifically designed program, reads2snps (Gayral et al., 2013), to predict genotypes with the parameters "-par 1 -th1 0.99 -min 30 -nbth 10". The program employs a maximum likelihood framework to reduce spurious heterozygous genotype calls caused by paralogous mapping. We retained only biallelic SNP positions present in at least 80% of individuals and removed contigs with mapping sites less than 100. To ensure the reliability of intraspecific polymorphism analysis, we retained variants in shared orthologous contigs among three focal species.

2.3. Patterns of molecular divergence among reproductive modes

To investigate the association between reproductive mode and fixation of weakly deleterious mutations, we estimated the ratio of nonsynonymous to synonymous divergence (d_N/d_S = ω) among species with sexual, asexual and mixed reproductive modes. We used the back-translated nucleotide alignments of the 362 19-species orthologs as the input for ω ratio estimates. Using the concatenated data sets, we first reconstructed the phylogenetic tree of the 19 species (Fig. 1a) by RAxML v.8.2.12 (Stamatakis, 2014) calculated with 1000 rapid bootstraps and the GTRGAMMAIX model estimated by ModelTest-NG v.0.1.6 (Darriba et al., 2020). Each alignment and the fixed topology of species tree were then parsed as input for the CodeML, which is included in the PAML v.4.9 package (Yang, 2007). We evaluated the fit of following four hierarchical branch models (Fig. 1a–d) to the data: 1) the one-ratio model (M₁), which assumes the same ω value for all branches; 2) a two-ratio model (M₂), which estimated two separate ω values for internal (ω_int) and external (ω_ext) branches; 3) a three-ratio model (M₃), in which three reproductive mode branches were allowed to have different ω with internal branches all treated as sexual branch; and 4) a four-ratio model (M₄) similar to M₃ but with the internal branches allowed to have a separate ω. We compared the fit of models by using likelihood ratio tests under specific degrees of freedom (df). We applied different models to both concatenated genes and individual alignments. We used the branch-site model to test if there any loci in asexual reproductive species present the signatures of positive selection.

2.4. Patterns of polymorphism among reproductive modes

Using our three focal species with contrasting reproductive modes and counts of the number of homozygous and heterozygous genotypes, we calculated the observed heterozygosity by dividing the number of heterozygotes by the total number of genotypes called. We averaged statistical parameters across orthologs weighted by sequence length.

To determine whether deleterious mutations accumulate at different rates in species with asexual and mixed reproduction compared with sexual species, we used three complementary approaches. First, we followed well-established precedence by assuming that synonymous mutations are selectively neutral (Kimura, 1968) and estimated the proportion of nonsynonymous mutations relative to silent mutations. Specifically, we extracted 0-fold degenerate sites as nonsynonymous sites and 4-fold degenerate sites as synonymous sites and inspected whether the frequency-dependent and independent measure of diversity (π and θ_W defined by Nei and Li, 1979; Watterson, 1975, respectively) of synonymous (π_S, θ_S), nonsynonymous (π_N, θ_N) and their ratios (π_N/π_S, θ_N/θ_S) differed among reproductive modes. We conducted McDonald-Kreitman (MK) tests (McDonald and Kreitman, 1991) to determine whether the ratio of nonsynonymous to synonymous (P_N/P_S) polymorphism differed from the ratio of nonsynonymous to synonymous divergence (D_N/D_S) with D. esculenta, D. rotundata and D. trifida used as outgroups for sexual, mixed and asexual ingroups, respectively. We quantified the direction and degree of departure from neutrality using the neutrality index (NI), the odds ratio of the MK contingency table, defined as (P_N/D_N)/(P_S/D_S) (Rand and Kann, 1996).

Second, we tested whether harmful mutations were represented in the categories of nearly neutral and mildly harmful mutations. To perform this analysis, we kept only 0-fold degenerate sites and classified nonsynonymous mutations into two categories: conservative and radical amino acid changes using the Conservative-Radical Index (CRI) (Sharbrough et al., 2018). CRI is an overall index of the degree of amino acid changes averaging seven classification schemes, including considerations of charge, polarity, volume (Zhang, 2000), correlation with K_A/K_S, charge and aromaticity, charge and polarity (Hanada et al., 2007) and charge and hydrophobicity (Grantham, 1974). Specifically, we modified the original codon table (for invertebrate mitochondrial genome) to fit our study system, and classified mutations as conservative and radical amino acid changes with CRI ≤ 0.5 and CRI > 0.5, respectively. Analogously, we calculated π and θ for conservative (π_C and θ_C), radical (π_R and θ_R) mutational changes and weighted the amount of the two categories of mutations by dividing with synonymous diversity (π_S, θ_S). We also extend this estimate of θ using only the variants of singletons (hereafter θ_U). Because instances of mutation were expected to occur more frequently and preserve longer in populations of species with asexual or mixed reproduction, this approach can be used to estimate the difference in mutation accumulation rate among reproductive modes.

Finally, as a complementary analysis, we conducted assessment of selection efficacy using the distribution of fitness effects (DFE) models (Eyre-Walker and Keightley, 2009). Although the assumption of random mating is likely to be violated in primarily asexual populations, previous studies have validated its effectiveness for sexual versus asexual comparisons (e.g. Hollister et al., 2015; Koufopanou et al., 2020), and also outcrossing versus selfing comparisons (Arunkumar et al., 2015; Wang et al., 2021). Therefore, we applied the DFE method to contrast variation in selection efficacy among our three focal species. We estimated the genome-wide DFE for new nonsynonymous mutations, the rate of adaptive evolution (α), and the proportion of adaptive (ω_A) and non-adaptive (ω_NA) nonsynonymous mutations by using DFE-α v.2.16 (Eyre-Walker and Keightley, 2009). We computed the locus-specific synonymous and nonsynonymous unfolded site frequency spectra (SFS) based on predicted genotypes of each species; two outgroup species (D. esculenta and D. rotundata for sexual D. subcalva; D. fordii and D. rotundata for mixed D. polystachya; D. trifida and D. pentaphylla for asexual D. melanophyma) were independently and jointly used to polarize the variant loci of our focal species. We obtained the estimates of DFE and α for each of three species with the assumption of neutral evolution of synonymous sites. Finally, we summarized the proportion of mutations from the DFE fall in three discrete categories of N_es (N_es < −10, −10 < N_es < −1, and −1 < N_es < 0) with s being the estimated selective disadvantage of the mutation. We performed 200 bootstrap replicates for each species datasets by randomly resample the loci, and generated mean values and 95% CIs.

3. Results 3.1. Patterns of divergence in clades with contrasting reproductive systems

The reconstructed phylogeny of 19 Dioscorea species indicated that sexual species were located at the basal position of the tree topology, with asexual species strongly supported as monophyletic (ML bootstrap value: 100%) and derived. Species with mixed reproduction were generally intermediate between those with sexual and asexual reproduction in their positions on the tree (Fig. 1a).

In the concatenated data set, the ω value under the one-ratio model (hereafter M₁; Table S4) was 0.239, indicating that most Dioscorea nuclear genes evolve under purifying selection. The two–ratio model (M₂) provided a significantly better fit than M₁ (χ² = 327.72, df = 1, P < 0.001; Table S4), with a higher ω value in external (ω_ext = 0.265) than internal (ω_int = 0.210) branches. This result is consistent with the hypothesis that most mildly deleterious mutations are removed over long-time scales but slightly deleterious mutations are relatively more persistent over shorter time scales. M₃ exhibited a significantly better fit than both M₁ (χ² = 175.80, df = 2, P < 0.001; Table S4) and M₂ (χ² = 151.92, df = 1, P < 0.001; Table S4). In this model, asexual branches had higher ω values (0.338) relative to branches of the other two types of reproductive mode (0.233 and 0.239 for sexual and mixed reproduction branches, respectively). This result suggests that asexual species experience a higher rate of nonsynonymous mutation accumulation and less effective purifying selection compared with their sexual or mixed reproductive counterparts. The M₄ model (Fig. S2) provided a similar pattern of ω variation in internal-external branches and among three external branches, a result concordant with the results for M₂ and M₃.

Out of the 362 single-copy orthologs, we removed 109 orthologs that either had ω > 1 or lacked any synonymous or non-synonymous substitutions. This procedure left 253 orthologs that we analyzed independently. For total genes, ω values for asexual and external branches were significantly higher resulting in a positive difference in (ω_ext – ω_int) (Wilcoxon-signed rank test, V = 37, 895, mean = 0.040, P < 0.001) (Fig. 2a and b) and (ω_asexual – ω_sexual) (V = 37, 726, mean = 0.080, P < 0.001) (Fig. 2c and d). The ω values for asexual branches were also significantly higher than branches for species with mixed reproduction (V = 35, 894, mean = 0.055, P = 0.018) (Fig. 2e and f), whereas there was no significant difference between ω values of sexual and mixed branches (V = 34, 158, mean = 0.025, P = 0.191, Fig. 2g and h). Overall, average ω across the 253 loci was significantly higher in clades with asexual (mean = 0.292) versus sexual (0.213) or mixed (0.237) reproductive modes (P < 0.001 and P = 0.006, respectively), and importantly, ω was never elevated in clades with mixed reproduction compared with those with sexual reproduction, either including (P = 0.146) or excluding (P = 0.745) the orthologs (31) with ω > 1 (Fig. 3a). Furthermore, the branch-site model test indicated that in asexual species only 1.58 % of the fast-evolving genes exhibited a signature of positive selection (Table S5).

3.2. Patterns of polymorphism and efficacy of selection in species with contrasting reproductive systems

The predominantly asexual reproduction Dioscorea species (D. melanophyma) was characterized by elevated heterozygosity, but reduced sequence diversity (Table 1). Asexual individuals exhibited a significantly higher level of heterozygosity (mean ± SD: 0.285 ± 0.040, P < 0.001) than D. polystachya, with mixed reproduction (0.142 ± 0.045) and D. subcalva, that reproduced primarily by sexually produced seed (0.190 ± 0.080), with the latter two species not significantly different in heterozygosity (P = 0.069). Average genome wide nucleotide diversity was significantly reduced (π_Total = 0.0021) in the asexual species compared with the sexual species (π_Total = 0.0058, P < 0.001) and the species with mixed reproduction (π_Total = 0.0047, P < 0.001). Similar trends were evident at nonsynonymous sites (π_N), although the relative reduction in diversity of synonymous sites (π_S) was less severe in the asexual species. These patterns resulted in the highest average ratio of nonsynonymous to synonymous nucleotide diversity (π_N/π_S = 0.378) in the asexual species, intermediate values for the species with mixed reproduction (π_N/π_S = 0.344) and the lowest average values for the sexual species (π_N/π_S = 0.297) (Fig. 3b). The distributions of θ_W variation for synonymous (θ_S) and nonsynonymous (θ_N) sites were similar to those of π, and thus a similar trend of θ_N/θ_S variation to π_N/π_S among three reproductive modes (Table 1).

In our comparisons of the three focal species, we found a significant effect of reproductive mode on the neutrality index (NI). The value in the sexual species was close to 1 (NI = 0.931; Fisher's exact test, P = 0.070), which fits the null expectation of equal ratios of polymorphism and divergence. In contrast, both the species with asexual and mixed reproduction exhibited significant deviations from neutrality (asexual: NI = 1.488, P < 0.001; mixed: NI = 1.265, P < 0.001), thus providing evidence for the accumulation of harmful mutations.

To further investigate how different levels of harmful mutations accumulate within each of the reproductive modes, we partitioned nonsynonymous mutations into conservative and radical amino acid changes. We observed consistent patterns in conservative and radical amino acid polymorphisms π across the three species with contrasting reproductive modes (Table 2); asexual Dioscorea melanophyma exhibited the highest level (π_C/π_S = 0.229, π_R/π_S = 0.187), followed by D. polystachya, with mixed reproduction (π_C/π_S = 0.183, π_R/π_S = 0.142), whereas sexual D. subcalva had the lowest levels (π_C/π_S = 0.152, π_R/π_S = 0.114) (Fig. 4a and b).

We further compared the levels of nucleotide diversity of singletons across reproductive modes and mutational types (Table 2). This analysis revealed that asexual D. melanophyma exhibited significantly higher levels of both private conservative amino acid polymorphism (θ_U/θ_U = 0.279) and private radical amino acid polymorphism (θ_U/θ_U = 0.241) than the species with mixed (θ_U/θ_U = 0.252, θ_U/θ_U = 0.213) and sexual reproduction (θ_U/θ_U = 0.229, θ_U/θ_U = 0.189) (Fig. 4c). Together, these results indicate that both conservative and radical nonsynonymous changes are eliminated from the genomes of the sexual species more rapidly than in the species with asexual and mixed reproduction. D. polystachya, the species with mixed reproduction, also provided evidence for a reduction in selection intensity against nonsynonymous polymorphism, although the degree of reduction was smaller than in asexual D. melanophyma.

Our estimates of the DFE for nonsynonymous mutations using the unfolded synonymous and nonsynonymous SFS (Fig. S3) revealed striking differences in the distribution of fitness effects among our three focal species with different reproduction modes (Fig. 5). Sexual D. subcalva accumulated relatively fewer effectively neutral nonsynonymous mutations and weakly deleterious mutations with 10.92–12.75% and 7.62–7.92% falling into the range −1 < N_es < 0 and −10 < N_es < −1, respectively (Fig. 5a). In contrast, a significantly higher proportion of sites were classified as nearly neutral and weakly deleterious mutations in both the species with asexual and mixed reproduction (−1 < N_es < 0: 16.22–16.29% – asexual, 15.44–15.99% – mixed; −10 < N_es < −1: 12.13–12.29% – asexual, 11.14–11.53% – mixed) (Fig. 5a). These patterns were consistently recovered by using the two outgroup species for each of three ingroup species.

Using two outgroups for polarization, we found qualitatively similar results, for both the species with asexual and mixed reproduction, with higher frequencies in −1 < N_es < 0 and −10 < N_es < −1 bins (Fig. 5b), consistent with a genome-wide accumulation of slightly deleterious mutations in these species. The DFE inferred from folded SFS was largely consistent with the unfolded inference except that there were no significant differences in the proportion of −1 < N_es < 0 between sexual and mixed reproducing species (Fig. S4). Together, these results suggest that the transition from sexuality to asexuality has been accompanied by relaxed purifying selection.

We examined whether the efficacy of positive selection was influenced by the evolution of asexuality using three indices of adaptive substitution rate implemented in DFE-α v.2.16 (Eyre-Walker and Keightley, 2009). The proportion of nonsynonymous substitutions that were adaptive (α) was greater in sexual Dioscorea subcalva than in asexual D. melanophyma and mixed reproductive D. polystachya (Table 3). By splitting nonsynonymous mutations into adaptive (ω_A) and non-adaptive proportions (ω_NA), we observed a lower ω_A and a higher ω_NA in both the species with asexual and mixed reproduction compared with the species that reproduced primarily by sexual means (Table 3).

4. Discussion

We examined associations between reproductive mode and between-species divergence and within-species polymorphism in a sample of wild yam species with contrasting degrees of sexual and clonal reproduction. We found evidence that species that were primarily asexual exhibited a reduced efficacy of selection and the accumulation of deleterious mutations for both divergence and polymorphism. In contrast, species with mixed reproductive strategies involving both seed and clonal reproduction showed no evidence of an increased fixation of harmful mutations at divergence level, although an increase nonsynonymous mutation was evident at polymorphism level. Below, we compare our results with other recent studies of sexual and asexual taxa and evaluate hypotheses that might account for the associations we observed.

4.1. Influence of predominant asexuality on molecular divergence

Our phylogenetic analysis provided evidence that the predominantly asexual clade of Dioscorea species that we investigated evolved from sexual ancestors with the mixed reproductive clade as an intermediate state (Fig. 1a). We detected an elevated signature in the ratio of nonsynonymous to synonymous mutations (d_N/d_S) in the asexual clade (0.338 and 0.292 for concatenated and individual ortholog, respectively) compared to its sexual counterpart (0.233, 0.213) (Figs. 2 and a). This result supports the prediction of less effective purifying selection in asexual than sexual species and higher rates of deleterious mutation accumulation over long time scales.

In our study, the extent of molecular evolutionary divergence between sexual and asexual species of Dioscorea (increase of asexual over sexuals: 45.06% and 37.09% for concatenated and individual ortholog, respectively) was greater than had been previously reported in comparative studies of plant taxa, including Oenothera (range 4.1–20.6% in different clades; Hollister et al., 2015), Boechera (0–3.1% depending on populations and sequence site categories; Lovell et al., 2017) and the Ranunculus auricomus complex (only outlier genes involved in sporogenesis and gametogenesis with elevated d_N/d_S ratios; Pellino et al., 2013). The sexual species of both Oenothera and Boechera exhibit partial self-fertilization, which may potentially decrease the extent of divergence between sexual and asexual species, because selfing can tend to weaken the effectiveness of purifying selection and result in genetic consequences similar to those observed in asexual reproduction (Sicard and Lenhard, 2011; Cutter, 2019; Tsuchimatsu and Fujii, 2022). In our study, the ancestral state of Dioscorea is sexual reproduction and significantly all species in the genus are dioecious thus enforcing obligate outbreeding. The dioecious sexual system may therefore contribute towards the contrasting evolutionary divergence between sexual versus asexual clades of Dioscorea.

The occurrence of time-dependent increases in the proportion of fixed deleterious mutations within asexual species, as predicted under Muller's ratchet (Muller, 1964; Lynch et al., 1993), probably explains the wide differences in d_N/d_S values reported in earlier studies. These ranged from relatively large versus small differences to no signature of elevated divergence. Comparative studies conducted at the genus level, e.g. Oenothera (Hollister et al., 2015) and Dioscorea (this study) have revealed significantly larger increases in ω values in these asexual clades than in studies of asexual populations/lineages that were conducted at the species level, such as in Boechera spatifolia (Lovell et al., 2017) and Ranunculus auricomus complex (Pellino et al., 2013). This suggests that the extent of molecular divergence associated with asexuality may largely depend on the time since differentiation between sexual and asexual taxa.

It is notable that the earlier plant studies of molecular divergence mentioned above involved species with asexual seed production via apomixis (Boechera, Ranunculus) or permanent translocation heterozygosity (Oenothera), rather than species with prolific clonal propagation, as in some of the Dioscorea species we investigated. The maintenance of meiotic gamete production in these other asexual species might allow purifying haploid selection to eliminate deleterious mutations expressed during the gametophytic stage (Immler, 2019; Beaudry et al., 2020). Asexual reproduction by vegetative propagation circumvents meiosis and gametophyte development and it is possible that this could result in higher levels of deleterious mutation accumulation in clonally reproducing species. A comparative study of molecular divergence in apomictic versus clonally propagated species would be worthwhile to better understand the consequences of different modes of asexuality on plant genomes.

4.2. Effects of partial asexuality on molecular evolution

In contrast to our results for primarily asexual species, our polymorphism data indicated significant amounts of nonsynonymous variants segregating in populations with mixed reproductive systems (Fig. 3b and Fig. 4; Tables 1 and 2), but with no overall increase in the fixation of these mutations over long-time scales (Fig. 2 and Fig. 3a). Several explanations may account for this unexpected pattern.

It is possible that asexuality in the mixed reproductive clade is evolutionarily too young to change polymorphic mutations into fixed divergence and hence d_N/d_S ratios are still similar to those in the sexual clade. However, this explanation seems unlikely because our phylogenetic analysis of Dioscorea indicated that the mixed reproductive clade evolved prior to the asexual clade (Fig. 1a) and the latter exhibited substantial amounts of divergence compared to the sexual clade. An alternative explanation assumes that species with mixed reproductive systems do not accumulate significant amounts of genome-wide deleterious mutations because facultative sexuality prevents pronounced levels of mutation accumulation. Theoretical studies indicate that novel mutations in facultatively asexual populations accumulate more slowly than in populations of obligate asexuals in which each mutation is added to the mutational load (Hojsgaard and Hörandl, 2015; Hartfield, 2016). This idea of "a little bit of sex serves to avoid Muller's ratchet" has received some support in a polyploid species (Hodač et al., 2019), although ineffective selection has been reported in facultatively asexual Spirodela polyrhiza (Ho et al., 2019; Wang et al., 2024), and clearly the long-term effects of facultative asexuality still remain to be confirmed empirically in plants. A recent study of the water flea Daphnia, a cyclical parthenogen with prolific asexual reproduction, documented a reduced efficacy of selection associated with its mixed reproductive strategy (Fields et al., 2022).

Another mechanism that may contribute to the avoidance of mutation accumulation in species with mixed reproduction involve demographic factors associated with large census population size. Reproduction and dispersal by both seeds and bulbils in Dioscorea species with mixed reproduction may be an important factor in promoting the larger geographical distribution and population density of most species (Sugihara et al., 2021; Chen et al., 2022), and these ecological factors may help to counteract the negative genetic impacts of asexuality. Thus, in species with mixed reproduction several factors operating in concert may limit the pace of divergence and genomic decay, including the amount and history of facultative sexuality in lineages and the demographic characteristics of individual species.

4.3. Patterns of genetic diversity and mutation accumulation accompanying asexuality

Some asexual taxa are characterized by recent polyploid-hybrid origins (e.g. Oenothera, Hollister et al., 2015; Boechera spatifolia, Lovell et al., 2017) and also the introgression of genomic regions (e.g. Pennisetum/Cenchrus, Ozias-Akins et al., 2003) or entire chromosomes (e.g. Daphnia pulex, Tucker et al., 2013) from related species, which may account for the finding of increased heterozygosity in asexuals. If the origins of asexual and mixed reproductive clades of Dioscorea also involved these processes, we might expect similar trends in patterns of variation between nucleotide diversity and heterozygosity. However, we observed a sequential increase in heterozygosity from mixed reproduction to predominantly asexual species relative to their sexual counterparts, but the opposite trends for nucleotide diversity. The increase in heterozygosity may be due to fixation of inherited nucleotide differences from sexual progenitors and the subsequent accumulation of mutations associated with a lack of recombination and segregation across asexual genomes. Our results therefore provide evidence for the hypothesis that transitions between sexual and asexual reproduction in Dioscorea directly influence levels of genetic diversity and individual heterozygosity in natural populations and are likely to be less influenced by the confounding effects of polyploidy and hybridity.

For predominately asexual and mixed reproductive species, the accumulation of nonsynonymous mutations we detected involved both nearly neutral and moderate harmful variants. This pattern was supported independently by the analysis of the biochemical properties of amino acid changes (Fig. 4; Table 2) and DFE estimates (Fig. 5). For example, the DFE results for focal species indicated that there were significantly more sites in the N_es: −1 to 0 and −10 to −1 categories in both the asexual species and the species with mixed reproduction compared to the sexual species (Fig. 5). This result matched closely the comparative results of genetic load quantified for conservative and radical genetic burden (Fig. 4) and could be indicative of less effective selection occurring on mutations with a wide range of fitness effects, including nearly neutral mutations associated with reduced N_e. The accumulation of mildly deleterious mutations is more likely to be associated with Muller's ratchet because of the absence of recombination in asexual genomes. Interestingly, a previous study of the freshwater snail (Potamopyrgus antipodarum) detected the accumulation of radical amino acid changes but not conservative variants in an asexual lineage (Sharbrough et al., 2018). Insufficient genomic sampling may have accounted for the earlier finding of no signal of elevated conservative mutation accumulation as a limited range of variants from the mitochondrial genome were involved (Sharbrough et al., 2018).

Finally, the patterns of mutation accumulation during transitions in reproductive mode from sexual to asexual reproduction are likely to differ from those involving transitions from outcrossing to selfing. Previous molecular studies reported that relaxed purifying selection and increased purging of deleterious genes in selfers were prominent features of genome evolution (Arunkumar et al., 2015; Wang et al., 2021). These results were explained as the consequence of the reduced N_e and stronger selection efficacy against recessive mutations in the highly homozygous genomes of selfers. As would be expected theoretically, we did not detect the effects of purging in asexuals (Fig. 5) presumably because of the absence of recombination and masking of deleterious recessive mutations in heterozygotes.

5. Conclusions

Our study represents the first investigation of divergence and polymorphism at the genomic level in plant species with primarily sexual reproduction via seeds, predominantly clonal propagation and mixed reproductive modes. Our analyses confirmed the prediction of relaxed purifying selection and accumulation of deleterious mutations following the transition from sexuality to asexuality. Transitions from sexuality to asexuality are often associated with the irreversible loss of sexual trait function (reviewed in Eckert, 2002) and also an associated risk of extinction (Lynch et al., 1993). Significantly, in our focal asexual species – Dioscorea melanophyma – we have observed widespread pollen sterility (W. Zhou unpubl. observations) suggesting that this species may be relying exclusively on bulbil production for the maintenance of populations and colonization of new sites. In some ways this form of clonal propagation 'mimics' sexual seed production and may in the long-term buffer species against extinction risk. Indeed, Dioscorea species with mixed reproductive strategies involving bulbil production could experience the ecological benefits of dual reproductive modes without the inexorable accumulation of mutations associated with clonal reproduction over long-time scales. Mixed reproductive strategies involving seed and clonal propagation predominate in herbaceous perennial plants and their prevalence may be, in part, because they offer not only ecological advantages but also genetic benefits.

Acknowledgements

We thank Spencer C.H. Barrett for helpful comments on the article. This research was supported by the Top-notch Young Talents Project of Yunnan Province (YNWR-QNBJ-2019-203), the National Natural Science Foundation of China (32470394), and the Key Basic Research Program of Yunnan Province (202201AS070057, 202101BC070003, and 202103AC100003).

Data availability

RNA-seq data are deposited in the National Genomic Data Center of the China National Center for Bioinformation (GSA: CRA012445) with accession numbers provided in Table S1.

CRediT authorship contribution statement

Xin Wang: Investigation, Formal analysis, Data curation. Qing-Hong Feng: Investigation, Formal analysis, Data curation. Zhi-Hua Zeng: Formal analysis. Zhi-Qiang Zhang: Investigation. Jie Cai: Investigation. Gao Chen: Investigation, Conceptualization. De-Zhu Li: Investigation. Hong Wang: Conceptualization. Wei Zhou: Writing – review & editing, Writing – original draft, Investigation, Formal analysis, Data curation, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.The author Wei Zhou is an Editor for Plant Diversity and was not involved in the editorial review or the decision to publish this article.

Appendix A. Supplementary data

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