Genomic consequences of residual recombination in a hybrid apomictic hickory complex

Domain: neuroscience
Published: 2026-02-05
Editorial outcome: minor_revision

Abstract
Apomixis, a form of clonal asexual reproduction in plants, is often accompanied by residual sex, yet its genomic consequences remain poorly understood. Here, we assembled a haplotype-resolved genome of Carya hunanensis and analyzed whole-genome sequencing data from 195 adults and 180 mature embryos across four hickory species, representing a hybrid apomictic complex with both sexual and asexual lineages. We find apomictic species exhibited genomic signatures of clonality, notably loss of heterozygosity (LOH), suggesting recombination induced by rare sexual events. Despite harboring more heterozygous deleterious variants, apomictic adults showed lower realized mutation loads, particularly in hybrid C. hunanensis, whose apomictic haplotype disproportionately carried deleterious alleles. Remarkably, rare embryos from apomicts underwent recombination-mediated LOH, exposing deleterious mutations to selection. These findings reveal the genetic cost of residual sex, while also indicating its role in generating novel genotypes, supported by close relatedness among adult apomicts. Our study provides a unique genomic snapshot of how residual sex and recombination mitigate mutation accumulation and potentially facilitate clonal maintenance in natural asexual systems.

Introduction
Sexual reproduction, characterized by meiosis and fertilization, dominates in animals and flowering plants. While facilitating genetic recombination and biparental inheritance, it also comes with substantial evolutionary costs, including energetic demands, disruption of co-adapted gene complexes, and the halving of parental genome transmission1,2. In contrast, obligate asexuality avoids these burdens3 but is associated with other disadvantages. Transitions between the two modes of reproduction are rarely abrupt and often involve intermediate reproductive modes that combine asexuality with some recombination4,5. This residual recombination can lead to loss of heterozygosity (LOH), exposing recessive deleterious alleles and may reducing long-term viability6,7. Despite this, some asexual lineages have evolved mechanisms to mitigate genetic load and maintain heterozygosity8–10. However, these mechanisms are often system-specific, underscoring the importance of investigating diverse asexual reproductive strategies.

Apomixis is a form of asexual reproduction in plants, in which seeds are produced without meiosis (or its avoidance, apomeiosis) and fertilization11–13. It is classified as either gametophytic or sporophytic, depending on whether unreduced embryo sacs originate from gametic or somatic cells14. Most apomictic taxa are facultative, combining sexual and asexual modes and often retaining features like partial meiosis or pollination-induced endosperm15,16. As a result, maternal plants can sometimes produce both sexually and asexually derived seeds within a single generation, even within a single fruit in multi-ovulate species5,17, highlighting the complexity and variability of apomictic reproduction. Polyploidy is prevalent among apomictic taxa of hybrid origin, a combination though to confer genome buffering and increase environmental resilience18–20. In contrast, diploid hybrid apomicts, which have been reported rarely in Borella21,22 and some Citrus species23,24, may be more vulnerable. Despite its adaptive advantages such as reproductive assurance, local genotype fixation, and increased colonization ability, apomixis remains taxonomically rare, occurring in fewer than 0.15% of flowering plant species (roughly 400 out of 300,000)16,25,26.

This paradox has led to hypotheses about evolutionary constraints on apomictic lineages27. One of the main hypotheses points to its developmental origin: apomixis may arise from incomplete or unstable reprogramming of the conserved sexual pathway, making it prone to disruption28–30. Theoretical models, including Muller’s ratchet and the Hill–Robertson effect, further predict that reduced recombination in asexuals weakens purifying selection, allowing mildly deleterious mutations to accumulate over time and ultimately reducing evolutionary potential31–33. However, empirical comparisons between apomictic and sexual lineages have yielded mixed results: some asexuals show increased genetic load22, while others do not or even exhibit genome homogenization34,35. Moreover, although “a bit of sex” might mitigate or prevent Muller’s ratchet36,37, the influence of facultative apomixis is often overlooked in genomic studies22,38. If occasional sexual reproduction helps maintain apomictic lineages, one would expect to find multiple distinct but closely related clonal lineages (e.g., parent-offspring or sibling-like structure) within a species, an observation that remains underexplored. There is therefore a need to clarify the genomic consequences of residual recombination in apomicts, including that arising from rare sexual reproduction or historical hybridization, and to assess how such processes influence lineage persistence.

A group of closely related hickory species, Carya cathayensis, C. dabieshanensis, C. hunanensis, and C. tonkinensis (Supplementary Fig. 1), offers a promising model for studying natural apomixis in diploid plants. All four species are endemic to southern China and valued for their pecan-like edible kernels, yet they exhibit allopatric distributions and remain largely confined to specific regions39–41. Carya cathayensis was first shown to exhibit apomixis through detailed embryological analyses showing nucellar embryogenesis42, and supported by controlled pollination experiments showing seed development without paternal contribution43. Polyembryony, a trait associated with sporophytic apomixis44,45, has been reported in both C. cathayensis and C. hunanensis46. Recent phylogenomic analyses indicate that C. hunanensis is a hybrid lineage derived from C. cathayensis, C. dabieshanensis, and C. tonkinensis47, a finding supported by its intermediate husk morphology (Supplementary Fig. 1). Two distinct morphotypes, TD and YL, have been identified in allopatric populations of C. hunanensis, suggesting potential cryptic lineage diversification. Although apomixis and hybrid origin appear to coexist in this complex, the facultative nature of apomixis remains unconfirmed due to the absence of population-scale genomic evidence.

In this study, we assembled a haplotype-resolved, chromosome-level reference genome for the diploid C. hunanensis (2n=2x=32; Fig. 1 and Supplementary Figs. 2–3), and conducted whole-genome resequencing of 131 newly collected leaf samples from adult trees, together with 64 publicly available genomes48, totaling 195 individuals across the four Carya species (Supplementary Data 1). We also generated a novel dataset by sequencing 180 mature embryos from three apomictic species (Supplementary Data 2). Using this integrated genomic dataset covering both vegetative and reproductive phases, we provide robust evidence for both interspecific hybridization and facultative apomixis within this species complex, and address the following key questions: 1) Does apomixis increase the accumulation of deleterious mutations relative to sexual relatives? 2) What is the relative contribution of apomictic and sexual parental lineages to genetic load in hybrid apomicts? 3) How does recombination-mediated LOH shape deleterious mutation patterns and affect lineage persistence?

Results
Haplotype-resolved genome assembly of Carya hunanensis and its hybrid origin supported by population structure
To generate a high-quality reference genome for the putative hybrid C. hunanensis, we produced 153.27 Gb (~218× coverage) of multi-platform sequencing data using Illumina, PacBio, and Hi-C technologies. The final assembly totaled 1.34 Gb (contig N50=11.31 Mb; scaffold N50=43.81 Mb) and exhibited high contiguity, coverage, and accuracy, with 97.13% of sequences anchored to 32 pseudochromosomes (Fig. 1; Supplementary Figs. 4–5 and Supplementary Table 1). Approximately 50.33% of the genome consisted of repetitive elements and a total of 78,799 protein-coding genes were predicted in the initial, unphased assembly (Supplementary Tables 1–2). This elevated gene count is consistent with the nearly doubled genome size (1.34 Gb) relative to the estimated haploid genome size (701.66 Mb) based on 19-mer analysis. Together with a high heterozygosity rate of 2.27% (Supplementary Fig. 3), these results indicated that both divergent haplotypes were represented. Mapping patterns revealed contrasting alignment preferences among the three putative parental species: C. cathayensis and C. dabieshanensis aligned predominantly to the first 16 pseudochromosomes (96.59% and 96.77%, vs. ~31% to the latter), whereas C. tonkinensis displayed the reverse (97.36% to the latter 16 vs. 28.52% to the first). These results enabled the partitioning of the C. hunanensis genome into haplotype E (Chr1A–16A) and haplotype W (Chr1B–16B), which was further confirmed by SubPhaser49 analysis (Fig. 1 and Supplementary Figs. 6–7). Gene prediction on the separated haplotypes yielded 38,635 and 39,070 genes (Supplementary Table 1), respectively, in line with the expected gene content of a diploid genome.

We analyzed whole-genome resequencing from 195 leaf samples (131 newly collected and 64 publicly available genomes) representing four Carya species in southern China: 55 C. cathayensis, 47 C. dabieshanensis, 57 C. hunanensis (TD and YL morphotypes), and 36 C. tonkinensis (Fig. 2a)47,48. The average sequencing depth ranged from 21.61× to 26.03×, with mapping to the C. hunanensis haplotype E genome (Supplementary Data 1). ADMIXTURE50 analysis at optimal K=6 identified six genetic clusters (Fig. 2b and Supplementary Fig. 8). Carya hunanensis separated into two morphotypes (TD and YL) along with two admixed individuals (Chudx01 and Chuzy01); C. cathayensis and C. tonkinensis each formed distinct clusters, while C. dabieshanensis comprised at least two primary subclusters. Notably at K=2, where C. cathayensis and C. dabieshanensis are lumped into one cluster and C. tonkinensis formed a separate group, C. hunanensis displayed ~47% ancestry from the former group and ~53% from the latter. The YL morphotype remained admixed at K=3–4 but formed an independent cluster by K=5 (Fig. 2b). PCA corroborated these findings, with PC1 (35.07% variance) positioning C. hunanensis intermediately between the three putative parental species, and PC2 separating the TD and YL morphotypes (Fig. 2c). These findings were further supported by significant pairwise divergence (FST and dxy) among the five taxa (Supplementary Fig. 9), and HyDe51 analysis confirming the hybrid origin of C. hunanensis (Z>344, P=0; ~41% ancestry derived from C. cathayensis and/or C. dabieshanensis, and ~59% from C. tonkinensis; Supplementary Data 3), with γ values close to 0.5, indicating limited hybridization events.

The chloroplast phylogeny of 196 individuals (including C. illinoinensis as the outgroup) resolved five clades corresponding to the three focal species and two C. hunanensis morphotypes (Supplementary Fig. 10). The TD morphotype and two admixed individuals clustered with C. dabieshanensis, while the YL morphotype formed a sister lineage within the same clade. Haplotype analysis revealed a single haplotype in all 55 sampled C. cathayensis individuals, likely reflecting its apomictic reproduction42, two haplotypes in C. dabieshanensis, four in C. hunanensis, and greater diversity in C. tonkinensis (Fig. 2d). These results suggest that C. hunanensis likely inherited its maternal lineage from a variant of C. dabieshanensis (8–10 bp divergence), with C. tonkinensis as the paternal lineage (Fig. 2d and Supplementary Fig. 10).

Polyembryony and genomic signatures of clonality and residual recombination in three apomictic hickory species
Our germination experiments with 120 mature seeds per taxa revealed polyembryony (a hallmark of sporophytic apomixis44,45) in C. cathayensis (34.26%), C. dabieshanensis (24.77%), and both morphotypes of C. hunanensis (>98%), whereas C. tonkinensis produced only single seedlings (Supplementary Fig. 11 and Supplementary Table 3). The widespread occurrence of polyembryony, particularly in the hybrid C. hunanensis where seeds yielded up to eight seedlings, strongly supports apomixis in these taxa.

Kinship analysis was performed on genome-wide SNP data from adult individuals comprehensively sampled to represent the full geographic distribution of each species. All 55 individuals of C. cathayensis formed a single clonal lineage, potentially reflecting obligate apomixis in all adult trees of this species (Fig. 3a). Similarly, 44 of the 47 individuals of C. dabieshanensis grouped into five clonal lineages (comprising 9, 21, 9, 2, and 3 individuals; Fig. 3b), and 55 of 57 individuals of C. hunanensis (including 44 of the TD morphotype and 11 of the YL morphotype; Fig. 3c) were grouped into two clonal lineages, all possibly derived from apomictic reproduction. Pairwise comparisons within C. dabieshanensis, involving the five clonal lineages and the three remaining individuals (Cdabmj02, Cdazf01, and Cdazf02), as well as within C. hunanensis, involving the two clonal lineages and the two remaining individuals (Chudx01 and Chuzy01), revealed first- to third-degree relationships, suggesting occasional sexual reproduction (Fig. 3b, c). In contrast, no clonality was detected in sexually reproducing C. tonkinensis, where most pairwise relationships (602/630) were more distant than third-degree (Supplementary Fig. 12).

The folded site frequency spectrum (SFS) analysis revealed distinct genetic signatures. In C. tonkinensis, the SFS exhibited a unimodal distribution characteristic of sexually reproducing species (Supplementary Fig. 13a). In contrast, C. cathayensis showed a bimodal distribution with a pronounced peak in the first and last bin (Supplementary Fig. 13b). Likewise, C. dabieshanensis and C. hunanensis exhibited multimodal SFS patterns, with the three largest clonal groups of C. dabieshanensis (21, 9, and 9 individuals) and both TD and YL morphotypes of C. hunanensis closely resembling the SFS of C. cathayensis (Supplementary Fig. 13c, d). These patterns indicate that apomixis leads to reduced genetic diversity and skewed allele frequency spectra.

Sexually reproducing individuals typically show variable and stochastic LOH patterns due to recurrent recombination and random mating. In contrast, several asexual lineages may exhibit shared LOH profiles within clonal groups, reflecting non-random and recurrent patterns along chromosomes. In principle, such patterns could arise from mitotic DNA double-strand break repair, which is generally expected to generate only very short and localized LOH tracts and is unlikely to occur simultaneously across multiple chromosomes52, or from automixis, a modified form of meiosis predicted to produce genome-wide but structured LOH with positional biases around centromeres or telomeres, although this process has not yet been rigorously confirmed in angiosperms4,53. Residual sexual reproduction within apomictic clones could also generate such patterns, as it would be expected to produce extended LOH tracts at variable chromosomal positions, and subsequent selfing or mating among closely related clonal lineages could give rise to novel LOH profiles distinct from those of the parental clones (illustrated schematically in Supplementary Fig. 14). To assess genomic LOH patterns, we analyzed heterozygosity in 50-kb sliding windows across the 16 chromosomes of 195 adults from four Carya species. Using a cutoff of <0.001, benchmarked against the sexual control C. tonkinensis (Supplementary Figs. 15–16), we found nearly identical LOH blocks in all C. cathayensis individuals and eight distinct LOH profiles in C. dabieshanensis that matched to its clonal groups (Supplementary Figs. 17–18), both may reflect recombination-derived LOH tracts subsequently fixed through clonal propagation. In C. hunanensis, LOH distributions were largely absent, likely due to its hybrid origin (Supplementary Fig. 19). Carya tonkinensis exhibited no shared LOH patterns, supporting its obligately sexual nature (Supplementary Fig. 20).

To assess whether new LOH could arise from residual sexual recombination, we performed low-coverage (10×) whole-genome sequencing of 180 mature embryos from three apomictic taxa, with 60 embryos per species sampled from three trees (Supplementary Data 2). Under the same threshold as above, only 8 (C. cathayensis), 11 (C. dabieshanensis), and 3 (C. hunanensis, YL morphotype) embryos displayed large LOH segments that differed from those in most other embryos (Fig. 3d–f), indicative of rare recombination events during embryo development and plausibly arising from residual sexual reproduction. Applying a more stringent threshold of heterozygosity <0.0005 yielded consistent LOH patterns in both adult and embryo datasets (Methods; Supplementary Figs. 21–26).

Comparison of genetic load between apomictic and sexual hickory species
To evaluate genetic diversity and selection efficacy among the four hickory species, we analyzed genome-wide SNP data from 195 adult individuals. As expected, the apomictic hybrid C. hunanensis exhibited the highest genome-wide nucleotide diversity (π=1.24 × 10⁻2) and individual heterozygosity (h=1.48 × 10⁻2), with similar values across its TD and YL morphotypes. This was followed by the sexually reproducing C. tonkinensis (π=9.62 × 10⁻3; h=4.17 × 10⁻3), while the apomictic species C. dabieshanensis (π=4.48 × 10⁻3; h=3.84 × 10⁻3) and C. cathayensis (π=3.40 × 10⁻3; h=3.73 × 10⁻3) had substantially lower values. Purifying selection efficacy was assessed using the ratio of heterozygosity at 0-fold (nonsynonymous) to 4-fold (synonymous) degenerate coding sites. Among the four species, C. hunanensis exhibited the highest absolute heterozygosity at both site classes (all P<2.2×10–16; Fig. 4a, b), likely reflecting the retention of divergent parental haplotypes. However, it showed the lowest 0-fold/4-fold ratio, followed by the other two non-hybrid apomictic species, C. cathayensis and C. dabieshanensis, while the sexually reproducing C. tonkinensis had the highest ratio (all P<2.2×10–16; Fig. 4c).

We further quantified mutation load in four functional categories, including synonymous (SYN), tolerated (TOL), deleterious (DEL), and loss-of-function (LOF) mutations, based on SNP data from the same 195 individuals and two outgroups (C. illinoinensis and C. kweichowensis). For this analysis, we used SIFT4G predictions based on C. hunanensis haplotype E reference genome (derived from apomictic lineages) as the primary reference. Consistent with its hybrid origin, C. hunanensis showed markedly higher heterozygosity and lower homozygosity across the four mutation types than the other three species (all P<2.2×10–16; Fig. 4d and Supplementary Fig. 27). Accordingly, this corresponds to a lower realized mutation load, which reflects the contribution of homozygous variants, whereas the excess heterozygous variants mainly contribute to the masked load. Some differences were also observed between its TD and YL morphotypes (Supplementary Fig. 28). Excluding the hybrid, C. tonkinensis displayed significantly lower heterozygosity and higher derived homozygosity than the two non-hybrid apomictic species (all P<2.2×10–16; Fig. 4d). Specifically, derived homozygous variants accounted for 98.13% of LOF and 97.58% of DEL mutations in C. tonkinensis, compared to 95.31% and 92.75% in C. cathayensis, and 93.64% and 91.20% in C. dabieshanensis), while haplotype W was inherited from the sexual parent (C. tonkinensis). All analyses are based on 57 adult trees of C. hunanensis. Notably, C. cathayensis carried more homozygous-derived and fewer heterozygous deleterious variants than C. dabieshanensis (P<2.2×10⁻16 to P=2.4×10⁻5; Fig. 4d). Gene Ontology (GO) enrichment analysis revealed that homozygous DEL and LOF mutations in the three apomictic species were enriched in multiple biological processes, including functions related to pollen recognition and reproductive development such as embryo, seed, and fruit development (Supplementary Data 4). Analyses based on Grantham scores54 and Haplotype W (derived from sexual lineages) produced concordant overall patterns, confirming the robustness of these results (Methods; Supplementary Figs. 29–30).

Inheritance of deleterious mutations in the hybrid apomict Carya hunanensis from sexual and apomictic parental lineages
To further explore the distribution of deleterious mutations inherited from sexual and apomictic parents, we compared the two haplotype-resolved genomes of the apomictic hybrid C. hunanensis: haplotype E, derived from the apomictic parents C. cathayensis and C. dabieshanensis, and haplotype W, derived from the sexual parent C. tonkinensis (Figs. 1–2 and Supplementary Data 3). Across 57 adult individuals of C. hunanensis, individual-level heterozygosity at 0-fold and 4-fold degenerate sites, as well as the 0-fold/4-fold ratio, were all significantly higher in haplotype E than in haplotype W (all P<2.2×10–16; Fig. 5a–c), indicating greater retention of variation in the apomictic-derived genome. Consistently, the ratios of deleterious to synonymous variants (DEL/SYN and LOF/SYN) at heterozygous sites were also significantly elevated in haplotype E (all P<2.2×10–16; Fig. 5d). In contrast, no significant differences were observed at homozygous sites between the two haplotypes (P=0.58 and 0.054; Fig. 5e).

Mutation load in LOH regions reveals genomic costs of occasional recombination events in apomictic hickory species
Loss of heterozygosity (LOH) can unmask recessive deleterious mutations, thereby reducing individual fitness and potentially driving selection against recombination in asexual lineages55,56. To assess these genomic consequences, we compared deleterious mutation loads between LOH and non-LOH regions (heterozygosity cutoff of 0.001) in adults of the non-hybrid apomictic C. cathayensis and C. dabieshanensis, both exhibiting extensive genome-wide LOH (Supplementary Figs. 17–18, 21–22). On average, C. cathayensis exhibited more LOH blocks than C. dabieshanensis (1973 vs. 1805 with P=1.1×10–6; Fig. 6a), consistent with its lower genome-wide heterozygosity. In both species, heterozygous DEL and LOF mutations were significantly reduced in LOH regions, typically with 0–6 sites. Notably, C. cathayensis harbored more heterozygous deleterious mutations than C. dabieshanensis (average 1.60 vs. 0.77 for DEL with P=2.3×10–7; 3.36 vs. 0.60 for LOF with P<2.2×10–16; Fig. 6b), reflecting differences in mutation accumulation. At homozygous sites, the ratio of deleterious to synonymous mutations (DEL/SYN) was significantly lower in LOH regions compared to non-LOH regions in both species (C. cathayensis: P=1.4×10⁻5 and C. dabieshanensis: P=2.4 × 10⁻5; Fig. 6c). However, homozygous LOF/SYN ratio showed divergent patterns: C. cathayensis had significantly higher ratios in LOH regions (P=3.6 × 10–9), while C. dabieshanensis had lower ratios (P=2.2 × 10–7; Fig. 6d), suggesting species-specific effects of recombination or its frequency.

To assess the genetic cost of recombination in apomictic hickories, we examined genome-wide mutation loads in 180 mature embryos of C. cathayensis, C. dabieshanensis, and the YL morphotype of C. hunanensis. Most embryos were genetically identical to their maternal plants, consistent with clonal inheritance, but a small subset showed genome-wide LOH patterns indicative of residual sexual reproduction (Fig. 3d–f). These sexual embryos displayed markedly reduced genome-wide heterozygosity (Supplementary Fig. 31) and increased homozygosity at deleterious mutation sites, including both DEL and LOF variants (P=1.5 × 10–8 to 2.8 × 10–6; Fig. 6e and Supplementary Figs. 32–33).

We next focused on LOH regions to test whether mutation load is further amplified within these vulnerable genomic segments. We compared the number of homozygous DEL and LOF mutations within ancestral LOH regions, which are shared by both sexual embryos and their apomictic siblings, and in novel LOH regions that occur only in sexually derived individuals. In most cases, sexual embryos carried substantially more homozygous deleterious variants within both ancestral and novel LOH regions (Fig. 6f and Supplementary Figs. 34–35). On average, sexual individuals of C. cathayensis exhibited 2.2–3.3-fold increases in homozygous DEL and 1.1–3.0-fold increases in LOF burden relative to apomictic individuals within LOH regions (Fig. 6f). The corresponding increases in C. dabieshanensis were 1.7–2.5-fold (DEL) and 1.2–2.1-fold (LOF) (Supplementary Fig. 34). The most pronounced effects were observed in the hybrid taxon C. hunanensis, with 3.6–45.6-fold (DEL) and 4.3–84.1-fold (LOF) increases in homozygous mutations within LOH regions (Supplementary Fig. 35). Nonetheless, inter-individual variability was observed: in one sexual embryo each of C. cathayensis and C. dabieshanensis, no novel LOF mutations were detected, and in one C. hunanensis embryo, the deleterious burden was only slightly elevated compared to its clonal counterparts (Fig. 6f and Supplementary Figs. 34–35).

Methods
Genome assembly, annotation, and haplotype phasing
To identify ploidy levels of the presumed hybrid taxon of C. hunanensis, fluorescence in situ hybridization (FISH) analysis was performed on root tip meristem tissues of the TD morphotype (Supplementary Figs. 1–2). Fresh young leaves of the TD morphotype of C. hunanensis (Tongdao County, Hunan Province, China; 26°6′ 53.63″N, 109°41′39.61″E) were collected from an adult tree. High-quality genomic DNA was extracted using the CTAB method. Illumina short-read (350 bp) and PacBio long-read (20 kb) libraries were sequenced on the NovaSeq 6000 and PacBio Sequel II platforms, respectively. Hi-C libraries were prepared and sequenced on the NovaSeq 6000. The quality-filtered Illumina short reads were used to estimate genome size, heterozygosity, and the repeat content of the genome by analyzing 19-mer frequencies using Jellyfish v2.370 and Genomescope v2.071 (Supplementary Fig. 3). High-accuracy CCS reads were used to assemble the draft genome with hifiasm v0.1672. Chromosome-level assembly was then performed using Hi-C data, with contigs clustered, ordered, and oriented using LACHESIS73. Misplacements and orientation errors were manually corrected based on chromatin interaction patterns, resulting in an assembly anchored to 32 chromosomes (Supplementary Fig. 4 and Supplementary Table 1). Genome assembly completeness was assessed by mapping Illumina whole-genome sequencing reads using BWA v0.7.1274 and by BUSCO v5.2.275 with the “embryophyta_odb10” database (Supplementary Fig. 5). Gene prediction was conducted through a combination of ab initio prediction, homology-based inference, and transcriptome data from RNA sequencing. Repetitive sequences, including transposable elements (TEs) and tandem repeats, were also identified using a combination of homology-based and ab initio approaches. See details of this section in Supplementary Methods.

The initial assembly obtained was twice the size of the anticipated genome size (Supplementary Fig. 3), suggesting that it contained two haplotypes. Since C. hunanensis is thought to be a hybrid between C. tonkinensis and the ancestor of C. cathayensis and C. dabieshanensis47, we separated these haplotypes using a method based on individual genome coverage. Illumina resequencing reads from five representative individuals of C. cathayensis, C. dabieshanensis, the TD and YL morphotypes of C. hunanensis, and C. tonkinensis were aligned to the 16 pseudochromosome pairs. The highest coverage from C. cathayensis and C. dabieshanensis was used to define homologous chromosomes as ChrA, with regions of lower coverage designated as ChrB (opposite coverage patterns were observed in C. tonkinensis). To validate the accuracy of haplotype phasing, haplotype-specific repetitive DNA sequences were identified using SubPhaser49 with default parameters. This analysis confirmed the haplotypes identified by the genome coverage-based method. Consequently, the C. hunanensis genome was divided into two haplotypes, referred to as C. hunanensis haplotype E (East) and C. hunanensis haplotype W (West) (Fig. 1), reflecting their inferred parental geographic origins.

Sample collection for leaf and fruit tissues and genome sequencing
To explore genetic details associated with apomixis, we conducted whole-genome sequencing on leaf tissue from 131 individuals across four Carya species. This included 41 C. cathayensis, 26 C. dabieshanensis, 36 C. hunanensis (comprising both TD and YL morphotypes), and 28 C. tonkinensis (Supplementary Data 1). Samples were collected from the natural distribution ranges of these species, ensuring that all sampled trees were mature (over 50 years old) and located in their native habitats. Field identification was primarily based on the morphological characteristics of leaves, fruits, and trunks. Genomic DNA was sequenced on the Illumina NovaSeq 6000 platform, generating 150 bp paired-end reads with an average depth of 30× for leaf samples. In addition to the newly sequenced data, we integrated previously published resequencing data from Zhang, Ding47 and Zhang, Yang48, which included 14 individuals of C. cathayensis, 21 C. dabieshanensis, 21 C. hunanensis and 8 C. tonkinensis. This resulted in a total sample size of 195 individuals: 55 C. cathayensis, 47 C. dabieshanensis, 57 C. hunanensis, and 36 C. tonkinensis (Supplementary Data 1). Furthermore, to detect potential recombination events mediated by residual sexuality in putatively apomictic taxa, we collected 180 mature embryo samples from three trees per taxon, with approximately 20 fruits per tree, for C. cathayensis, C. dabieshanensis, and the YL morphotype of C. hunanensis (Supplementary Data 2). These fully developed embryos from the fruit samples were sequenced using the same platform and protocol as the leaf samples but with an average depth of 10×.

Read mapping and variant calling
The raw reads from leaf and mature embryo samples were trimmed for adapters, bases with a quality score below 20, and any reads shorter than 50 bp using Trimmomatic v0.3876. The remaining high-quality paired-end reads were aligned to the C. hunanensis haplotype E reference genome using the BWA-MEM algorithm in BWA v0.7.1274 with default parameters, unless otherwise specified for mapping to haplotype W (see relevant sections below). SAMtools v0.1.1977 was employed to retain uniquely and properly paired reads, and to convert the Sequence Alignment Map (SAM) format to Binary Alignment Map (BAM) format. For both sample types, the SENTIEON DNAseq software package v202112.0578 was used to remove duplicate reads, realign indels, and perform single nucleotide polymorphism (SNP) calling for each individual. Joint SNP calling was conducted only for the leaf samples.

For both leaf and embryo samples, SNPs were filtered following the criteria established by Zhang, Cao58. The filtering process involved: (1) retaining only biallelic sites; (2) excluding sites with a mapping depth below one-third or above twice the individual’s average depth (with a minimum depth cutoff of 6× for leaf samples and 3× for embryo samples), while applying Q20 quality score thresholds and heterozygosity adjustments79; (3) retaining all SNPs regardless of missing rate by default, although in certain analyses, more stringent filters were applied, such as excluding SNPs with >20% missing data or removing all SNPs with any missing genotypes; (4) eliminating singleton variants; (5) excluding SNPs located within coding regions or their 20-kb flanking sequences to ensure neutrality; and (6) thinning SNPs by retaining a single site within each 20-kb window. For specific downstream analyses, custom SNP filtering steps were further applied based on the requirements of the software tools used and computational constraints.

Population structure and genetic differentiation analysis
To assess the genetic relatedness among individuals and taxa from four Carya species, Bayesian clustering was performed using ADMIXTURE50 based on 4,828 pruned SNPs from 195 individuals (filtered using steps 1–6, excluding SNPs with >20% missing data). ADMIXTURE employs a maximum likelihood-based model to estimate ancestral components in unrelated individuals. The number of genetic clusters (K) was predefined from 1 to 8, and cross-validation (CV) error was used to determine the optimal K value. Additionally, principal component analysis (PCA) was also conducted on the same SNP dataset using the R package SNPRelate v1.34.180 with default parameters. To quantify genetic divergence between species and within-species lineages (as inferred by ADMIXTURE and PCA), we estimated interspecific FST and dxy in non-overlapping 20-kb windows using pixy81. Since pixy requires both polymorphic and monomorphic (invariant) sites for unbiased estimation of π and dxy, we generated a separate VCF dataset for all 195 individuals using the GATK tools CombineGVCFs and GenotypeGVCFs tools in GATK v.4.4.0.082 with the “-all-sites” option enabled. After applying filtering steps 1–3 while retaining sites with <20% missing data, this pipeline produced a VCF dataset with invariant sites used for the estimation of genetic variation in pixy81.

Testing for hybridization events
Hybridization events were detected using HyDe software51, which employs Paterson’s D statistics based on ABBA-BABA site patterns. The parameters γ and 1-γ denote the inheritance probabilities from two parental sources. Generally, γ values near 0.5 indicate a recent or infrequent hybridization event, while values closer to 0 or 1 suggest an ancient hybridization event that has persisted in present-day species. The analysis was conducted in two phases using a SNP dataset of 23,395 sites from 195 individuals and the two aforementioned outgroups. Filtering steps 1–4 and 6 were applied, and sites with missing data were removed due to computational constraints. Initially, based on phylogenetic results from Zhang, Ding47, C. cathayensis and C. dabieshanensis were treated as a species complex, C. tonkinensis remained unchanged, and C. kweichowensis and C. illinoinensis served as outgroups to focus on hybridization within C. hunanensis. Subsequently, we extended the analysis to the aforementioned four species, using the same outgroups, to explorer reticulation among these lineages.

Maternally inherited plastid genome analysis
The high-quality paired-end reads from the 195 individuals, along with one individual of outgroup C. illinoinensis, were mapped to their respective species’ chloroplast genomes (http://cmb.bnu.edu.cn/juglans) using the BWA-MEM algorithm in BWA v.0.7.1274. After generating the realigned BAM files as described above, the ‘mpileup’ and ‘call’ functions in BCFtools v.1.1583, along with the ‘vcf2fq’ command from vcfutils.pl, were used to convert these BAM files into VCF format and derive consensus chloroplast genome sequences. From these sequences, 80 shared protein-coding genes were extracted across all individuals based on their corresponding chloroplast genome annotations. Phylogenetic trees were reconstructed from the concatenated alignments of the 80 chloroplast genes using IQ-TREE v.2.1.384 with the ‘-m MFP’ model and 1000 rapid bootstrap replicates. Additionally, haplotypes for all individuals were identified from the concatenated sequences using DnaSP v685, and a median-joining (MJ) haplotype network was constructed using PopART v.1.786.

Germination test for polyembryony identification
To examine the occurrence of the multi-seedling phenomenon, an indicator of polyembryony and a potential sign of sporophytic apomixis44,45, a germination experiment was conducted with four Carya species, including the newly identified TD and YL morphotypes of C. hunanensis (Supplementary Fig. 1). Fresh mature fruits were collected from six randomly selected trees per taxon, with 20 fruits harvested from each tree, resulting in 120 fruits per taxon. The fruits were stored at 4°C for 24 h, then soaked in water for 2–3 h at room temperature. Following this, the fruits were buried in peat soil at room temperature under a 1–2 cm soil layer. Germination was monitored every two days, recording the number of seeds germinated until no further germination occurred after 30 days. The total number of seeds exhibiting the multi-seedling phenomenon and its occurrence rate were calculated for each taxon.

Kinship inference and loss of heterozygosity analysis
We utilized the KING v.2.1 software87 to infer family relationships and assess inbreeding levels within each of the four hickory species. This software estimates kinship coefficients independently of sample composition or population structure, based on the difference between shared heterozygosity and shared homozygosity. For accurate kinship inference, it is recommended to use genome-wide SNP data without prior pruning or filtering of “good” SNPs that pass quality control87. To maximize mapping rates and SNP recovery, all adult individuals of C. cathayensis, C. dabieshanensis, and C. hunanensis were mapped to the C. hunanensis haplotype E reference genome, while all C. tonkinensis individuals were mapped to the C. hunanensis haplotype W reference genome (Supplementary Data 1). The initial filtering steps 1–3 were applied, retaining sites with missing data, resulting in datasets of 6,515,619 SNPs for C. cathayensis, 9,506,975 SNPs for C. dabieshanensis, 23,681,021 SNPs for C. hunanensis, and 20,859,321 SNPs for C. tonkinensis. Kinship coefficients were calculated using the ‘-kinship’ option with default settings in KING. Relationship classifications were assigned based on the ranges specified in the KING user manual, where kinship coefficients greater than 0.354 correspond to duplicate/MZ twins, 1st-degree relationships range from 0.177 to 0.354, 2nd-degree relationships range from 0.0884 to 0.177, and 3rd-degree relationships range from 0.0442 to 0.0884.

For genome-wide heterozygosity analyses, heterozygosity was calculated in 50-kb sliding windows for each of 195 adults and 180 embryos, using gVCF files as input to a custom Perl script. All samples were mapped to the C. hunanensis haplotype E reference genome to provide a consistent analytical framework, allow direct comparisons across tissues and reproductive modes, and minimize biases associated with using different references. Within each window, heterozygosity was calculated as the number of heterozygous sites divided by the total number of non-missing sites. To minimize technical artifacts, sites were excluded if sequencing depth fell outside an individual-specific range defined as one-third to twice the average coverage, but with a minimum cutoff of 6× for adults (~30× coverage) and 3× for embryos (~10× coverage). Multiallelic and indel sites were removed, and heterozygous genotypes were further filtered by allele balance79, requiring allelic ratios of 0.2–0.8 at higher depths (≥20× for adults, ≥6× for embryos) and 0.1–0.9 at moderate depths (≥10× for adults, ≥3× for embryos), otherwise a homozygous call was assigned. Finally, gVCF block coordinates were expanded per base to ensure alignment, and windows with >50% missing sites were reported as NA. Regions of loss of heterozygosity (LOH) were defined as 50-kb windows with heterozygosity <0.001, a cutoff benchmarked against the sexual control C. tonkinensis, where only 2.01% of 345,979 windows fell below this value (Supplementary Figs. 15–16). This sexual baseline provides a conservative threshold, minimizing false positives (rare even under continuous meiotic recombination) and avoiding circularity that could result from defining cutoffs on apomictic genomes enriched in homozygous tracts. After the stringent filtering described above, the rare windows with heterozygosity <0.001 observed in sexual likely reflect genuine regions of low diversity rather than technical artifacts.

Following the rationale of Weir, Capewell88 and Simion, Narayan89, which examined LOH to investigate mitotic recombination in asexual reproduction, we adapted this approach to assess meiotic recombination in apomictic hickories. In principle, sexual reproduction is expected to generate variable LOH patterns among individuals through recurrent recombination, whereas clonal reproduction should yield highly consistent LOH profiles within lineages. To evaluate whether the LOH detected in apomictic species could instead reflect residual sex, we contrasted our data with expectations from alternative mechanisms: mitotic repair, which typically generates short and localized tracts unlikely to occur simultaneously across multiple chromosomes52, and automixis, a modified meiosis predicted to produce genome-wide but structured LOH with positional biases4,53 (schematic comparison in Supplementary Fig. 14). For the 195 adults, species- and lineage-level LOH profiles were compared, with lineage-shared large blocks spanning all 16 chromosomes interpreted as potential recombination events that had subsequently become clonally fixed. To further test whether LOH could arise from residual sexual reproduction, we analyzed 180 embryos (60 per apomictic species). In this framework, the occurrence of a minority of embryos with LOH profiles distinct from the predominant clonal background was interpreted as evidence of rare recombination events arising during sexual reproduction. The genomic distribution of LOH blocks for each individual was visualized using a custom R script.

Folded site frequency spectrum and genetic diversity analysis
After identifying the reproductive modes and clonal lineages of the four hickory species, we assessed and compared their genomic variation. First, intraspecific nucleotide diversity (π) for adults from all four species (including the TD and YL morphotypes of C. hunanensis) was calculated using pixy81, incorporating both invariant and variant sites within a 20-kb sliding window. Second, individual heterozygosity (h) was calculated as the number of polymorphic sites divided by the total length of the C. hunanensis haplotype E reference genome, based on a filtered SNP dataset comprising 39,934,221 sites from 195 individuals (filtering steps 1–3). Third, we calculated the folded site frequency spectrum (SFS) for each taxon using a Perl script from Zhang, Cao58, based on 6,132,224 SNPs from the same set of individuals. SNPs with missing data were removed during filtering to minimize biases in allele frequency estimates and ensure accurate SFS reconstruction.

Evaluation of genetic load across sexual and asexual reproductive types
We calculated two metrics to quantify and compare the selection efficacy and deleterious genetic load across species. First, we determined the genome-wide ratio of heterozygosity (heterozygotes/total genotypes) for 0 to 4-fold degenerate sites within coding regions, based on the C. hunanensis haplotype E reference genome annotation. We then calculated the ratio of 0-fold nonsynonymous to 4-fold synonymous degenerate sites for each individual. Second, we annotated the effects of SNP variants (56,274,446 SNPs from 195 individuals and two outgroups, filtered using steps 1–3) on protein-coding genes, classifying them as synonymous, missense, or loss-of-function (LOF) variants using SnpEff v.5.090. LOF variants were defined as those involving a gain and/or loss of stop codons, as well as the loss of start codons. Missense SNPs were further classified as deleterious (score ≤ 0.05) or tolerated (score > 0.05) using SIFT4G program with custom databases built from the UniRef90 protein dataset91. To avoid reference bias in identifying derived alleles and deleterious variants, we determined the polarity of variants only when the two outgroups, C. illinoinensis and C. kweichowensis, shared identical homozygous states. For each SNP position in the SNP dataset, we determined the ancestral and derived allelic states using Est-SFS software92, applying a probability cutoff of 0.95 and assuming a Rate-6 substitution model. The ancestral alleles were assumed to be non-deleterious, while the derived alleles were considered potentially deleterious. We calculated the total number of derived alleles by counting heterozygous genotypes once and homozygous-derived genotypes twice. To verify robustness, we additionally estimated deleterious burden with Grantham scores54, which quantify the degree of amino acid change, designating missense mutations with scores ≥ 150 as deleterious and < 150 as tolerated93,94. We further repeated SIFT4G predictions on the same 197 individuals mapped to the alternative haplotype W reference genome (58,940,457 SNPs after identical filtering). Both approaches yielded results consistent with the SIFT4G predictions based on the Haplotype E reference; therefore, we retained the Haplotype E–based SIFT4G results in the main text and used them for subsequent analyses. Gene ontology (GO) enrichment analysis for genes associated with DEL and LOF mutations in apomictic hickory species was performed using TBtools95 with default parameters.

Comparison of inherited genetic loads between two genome haplotypes in Carya hunanensis
Building on our identification of hybrid diploid apomixis in C. hunanensis, we performed a comparative analysis of mutation burdens inherited from the two parental haplotypes. All 57 sampled adult C. hunanensis individuals were separately aligned to the haplotype E (derived from the apomictic parent) and haplotype W (derived from the sexual parent) reference genomes. After applying SNP filtering steps 1–3, we obtained 43,020,882 and 42,202,704 SNPs from 59 individuals (including two outgroups) for haplotype E and haplotype W, respectively. For each haplotype, we calculated the genome-wide nonsynonymous (0-fold) and synonymous (4-fold) site counts within coding regions and determined the nonsynonymous/synonymous (0-fold/4-fold) mutation ratio, following the same approach as previously described. Additionally, to better compare the difference in deleterious sites between the two haplotype genomes, we focused on the ratios of deleterious to synonymous mutations (DEL/SYN) and loss-of-function to synonymous mutations (LOF/SYN), considering both heterozygous and homozygous variants.

Quantification of mutation loads in LOH regions of adult trees and mature embryos
Loss of heterozygosity (LOH) can reduce individual fitness by unmasking recessive deleterious mutations in diploid genomes. This mechanism has been proposed as a contributing factor to the rarity of asexual organisms in nature and may drive selection for reduced recombination rates in asexual lineages55,56. Given the extremely low distribution of LOH regions in the hybrid apomictic species C. hunanensis, our analysis focused on C. cathayensis and C. dabieshanensis to investigate patterns of deleterious variation associated with LOH. Using whole-genome sequencing data from adult individuals, we applied a 50-kb sliding window with a heterozygosity cutoff of 0.001 applied to classify genomic regions into LOH and non-LOH categories. To quantify the deleterious mutation burden, we analyzed homozygous and heterozygous sites separately. For homozygous sites, we counted the number of predicted deleterious (DEL), loss-of-function (LOF), and synonymous (SYN) mutations per window, and calculated mutation load as the ratios DEL/SYN and LOF/SYN for each region type. For heterozygous sites, because heterozygous variants were frequently absent in many windows, which precluded reliable ratio-based normalization, we quantified mutation burden by directly comparing the raw counts of heterozygous DEL and LOF mutations between LOH and non-LOH regions.

To assess the extent of new deleterious mutation accumulation in sexual embryos relative to their maternal plants, we further analyzed SNP data from 180 mature embryos (60 per species) across the three apomictic species. Apomictic embryos were treated as genetic equivalents of their maternal trees, as their genotypes are expected to be identical due to the absence of meiosis and fertilization. Accordingly, we generated 26,760,972, 26,722,247, and 34,948,656 SNPs (filter steps 1–3) from 60 individuals and two outgroup samples for C. cathayensis, C. dabieshanensis, and the YL morphotype of C. hunanensis, respectively. First, we compared the overall mutation burden between apomictic and sexual embryos to evaluate the potential genomic costs associated with sexual reproduction. The sexual embryos were identified in the above-mentioned loss of heterozygosity analysis, and their genetic load was assessed using the same methods as applied to adult individuals. Considering sample size imbalance in the YL morphotype of C. hunanensis (57 apomictic vs. 3 sexual embryos), we used the Wilcoxon rank-sum test for this comparison, whereas Welch’s t-test was applied to adult comparisons (as described above) as well as to the other embryo contrasts (C. cathayensis: 52 vs. 8; C. dabieshanensis: 49 vs. 11). In addition, since LOH typically exposes recessive deleterious mutations to selection, we further classified LOH regions into two categories: (i) ancestral LOH, defined as regions shared by all embryos within the apomictic species, and (ii) novel LOH, referring to regions newly formed through recombination in sexual embryos. Within each LOH type, we quantified the number of homozygous deleterious (DEL) and loss-of-function (LOF) mutations. We hypothesized that sexual embryos exhibiting no significant increase in deleterious mutations within new LOH regions may represent novel genotypes with a higher likelihood of being retained by selection.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript investigates the genomic consequences of residual recombination in a hybrid apomictic hickory complex, examining clonal maintenance, mutation load, and the role of loss of heterozygosity (LOH). The study employs a combination of high-quality genome assembly, population genomic analyses of adult individuals, and embryo sequencing to address fundamental questions regarding apomixis and its long-term genetic implications. While the study presents a technically sound genomic analysis of Carya apomixis, it fails to meet the novelty threshold required for Nature Communications. The core findings represent incremental advances over prior work, and methodological limitations undermine key claims. I recommend rejection without resubmission for the reasons detailed below.\n1.Insufficient Novelty and Conceptual Advance. The central premise—that residual sexuality mitigates mutation accumulation in apomicts—is well-established across multiple systems. For instance, studies in polyploid apomicts like Ranunculus auricomus (Pellino et al., 2013, Mol. Ecol.) demonstrated purging via LOH; research in diploid systems such as Boechera spp. (Lovell et al., 2017, PLoS Genet.) quantified reduced mutation loads in apomicts through heterozygosity masking; and work on hybrid taxa in Citrus (Wang et al., 2017, Nat. Genet.) linked heterosis to deleterious mutation buffering. This study applies these established concepts to a new taxonomic group (Carya) without proposing novel evolutionary mechanisms, challenging existing paradigms, or providing cross-system insights that redefine apomixis theory.\n2.Methodological concerns undermine key conclusions. (a) LOH Definition and Validation: The classification of LOH as 50-kb windows with heterozygosity <0.001 is arbitrary and lacks justification. This threshold is not compared to sexual controls (e.g., C. tonkinensis) or validated against independent methods (e.g., allele frequency spectra), raising doubts about whether reported LOH regions truly reflect recombination or technical artifacts. (b) Sample Size Limitations: The inference of \"residual sexual reproduction\" relies on extremely small numbers of embryos with LOH (8 in C. cathayensis, 11 in C. dabieshanensis, 3 in C. hunanensis). This paucity of data makes it impossible to robustly quantify recombination rates or assess its evolutionary significance, weakening claims about \"costs and benefits of residual sex.\" (c) Mutation Load Calculation: The use of SIFT scores to classify deleterious variants without accounting for lineage-specific functional constraints (e.g., divergence between hickory haplotypes) may misclassify neutral or adaptive mutations as deleterious, skewing comparisons between apomictic and sexual lineages.\n3.Inconsistent Interpretation of Results. (a) It claims apomicts have \"lower realized mutation loads\" but acknowledges they harbor \"more heterozygous deleterious variants.\" This tension is attributed to \"hybrid buffering\" in C. hunanensis, but no mechanistic evidence (e.g., gene expression data showing masked recessive effects) supports this assertion. (b) LOH is proposed to \"expose deleterious mutations to selection,\" yet sexual embryos with LOH show higher homozygous deleterious loads. The suggestion that \"selection eliminates harmful alleles early\" is speculative and unsupported by data on embryo viability or fitness.", "reviewer_2": "The manuscript uses advanced kinship analysis of Carya species with dense SNP data and the KING algorithm to study clonality and possible apomixis. The authors apply the KING algorithm, which allows them to analyze unfiltered SNP datasets and look for duplicated (clonal), first, second, and third-degree kin relationships. The study claims that there is a lot of clonality, sometimes almost complete, within these species.\nBy adding resequencing data from adult tissues and sequencing 180 embryos, the authors try to answer old questions about hybridization, persistent lineages, and the effect of apomictic versus sexual reproduction in this complicated hickory group.\nThe study is relevant, the methodology is solid, and the evolutionary and ecological questions are important. Still, some points in the manuscript need to be explained better to make the conclusions clearer and reliable. Studying clonality on a genome-wide level in woody apomictic plants is new and important, including for the ecology and agriculture of hickory. The use of site frequency spectra (SFS) and loss-of-heterozygosity (LOH) analysis makes the argument for clonality even stronger and gives different views on the genetic structure.\nThe analysis shows clonality at the population level very well. Including 180 mature embryos from three apomictic species is a very valuable part of the study. The authors say that rare loss-of-heterozygosity (LOH) events in embryos are probably caused by leftover recombination from occasional sexual reproduction, not from mitotic mistakes. This explanation is reasonable and interesting, and it fits with the finding that there is no genome-wide LOH.\nHowever, it would be helpful to better explain the difference between LOH caused by mitosis and LOH caused by automixis. For example, this could be done by showing allele imbalance of LOH patterns on chromosomes. My main doubts comes from the fact that Carya is a sexually outcrossing (allogamous) species. Generally, in diploid species, the level of heterozygosity is roughly halved every generation when selfing occurs. In contrast, cross-pollination keeps or even increases heterozygosity in the population, preserving genetic variation.\nIn outcrossers, selfing causes higher LOH with a fast increase of homozygosity and a decrease of heterozygosity. Cross-pollination results in less LOH and keeps heterozygosity and genetic diversity high. This is why allogamous species that do repeated selfing experience a quick loss of heterozygosity and often suffer inbreeding depression.\nMy personal idea to explain the results (embryos from sexual reproduction showing LOH), together with the authors’ conclusion, is that in the crosses, both pollen and egg come from two plants of the same clone. But this should be studied more carefully. Still, seeing such LOH in only one generation (embryos from a mother plant) with crossing is unexpected. The authors should maybe sequence three biological replicates only for the samples that they think come from sexual reproduction, just to confirm that the LOH is consistent. I also suggest explaining this part more clearly. This is important because the embryo genotypes, especially related to LOH and recombination, are a large part of the study’s analysis and interpretation. So, this section should be better described better in the main methods or detailed more clearly in the Supplementary Information.\nAlso, I suggest adding more discussion about other possible reasons, like mitotic recombination or technical problems.\nIn conclusion, this paper gives valuable information about apomixis, clonality, and how hybrid species in long-living trees stay stable through time. It just needs some data/revisions, especially to explain the embryo part better."}, "author_response": "Reviewer #1 (Remarks to the Author):\nThis manuscript investigates the genomic consequences of residual recombination in a\nhybrid apomictic hickory complex, examining clonal maintenance, mutation load, and\nthe role of loss of heterozygosity (LOH). The study employs a combination of high-\nquality genome assembly, population genomic analyses of adult individuals, and\nembryo sequencing to address fundamental questions regarding apomixis and its long-\nterm genetic implications. While the study presents a technically sound genomic\nanalysis of Carya apomixis, it fails to meet the novelty threshold required for Nature\nCommunications. The core findings represent incremental advances over prior work,\nand methodological limitations undermine key claims. I recommend rejection without\nresubmission for the reasons detailed below.\nAnswer: We thank the reviewer for acknowledging that our study presents a technically\nsound genomic analysis of Carya apomixis. However, we respectfully disagree with\nthe assessment that the work lacks sufficient novelty. As outlined in the detailed point-\nby-point responses below, our analyses provide the first genome-wide evidence of\nrecombination-mediated LOH in a diploid woody apomictic complex and demonstrate\nits consequences for deleterious variation. We believe these findings have both\nevolutionary and applied significance. Methodological issues and questions regarding\ninterpretation of results have also been addressed through clarifications and\nsupplementary analyses, with corresponding revisions in the manuscript.\n1.Insufficient Novelty and Conceptual Advance. The central premise—that residual\nsexuality mitigates mutation accumulation in apomicts—is well-established across\nmultiple systems. For instance, studies in polyploid apomicts like Ranunculus\nauricomus (Pellino et al., 2013, Mol. Ecol.) demonstrated purging via LOH; research\nin diploid systems such as Boechera spp. (Lovell et al., 2017, PLoS Genet.) quantified\nreduced mutation loads in apomicts through heterozygosity masking; and work on\nhybrid taxa in Citrus (Wang et al., 2017, Nat. Genet.) linked heterosis to deleterious\nmutation buffering. This study applies these established concepts to a new taxonomic\ngroup (Carya) without proposing novel evolutionary mechanisms, challenging existing\nparadigms, or providing cross-system insights that redefine apomixis theory.\nAnswer: We note that the reviewer’s interpretation of the three cited studies does not\nreflect the conclusions of the original articles, all of which were already discussed in\nour manuscript (Introduction, lines 62 and 72–75; Discussion, lines 373–376).\n1. Pellino et al. 2013 (Mol. Ecol.) did not demonstrate “purging via LOH.”\nThe study analyzed RNA-seq data from the Ranunculus auricomus complex (2\ndiploid and 1 tetraploid sexual species, plus 2 apomictic hexaploid lineages) and\ntherefore could not assess LOH from transcriptomic SNPs. More importantly,\ntheir main conclusion was that there is no genome-wide difference in\ndeleterious mutation accumulation between apomictic and sexual lineages,\nbased on dN/dS comparisons of 1,231 genes from five representative\nindividuals.\n2. Lovell et al. 2017 (PLoS Genet.) reports the opposite of the reviewer’s claim.\nIn diploid/diploid-hybrid Boechera, Lovell et al. (2017) found that apomicts\nharbor more derived mutations than sympatric sexuals, and at\nphylogenetically constrained sites the ancestral (conserved) allele is more likely\nretained in sexuals—consistent with stronger purifying selection in sexual\nlineages. Thus, the paper supports greater (not reduced) mutation load in\napomicts. Notably, the study did not explicitly control for facultative apomixis\n(variable proportions of sexual/asexual seeds were present in 11 of 13 apomictic\nsamples, as shown in their Table S1), which can complicate inference about\npurely asexual dynamics.\n3. Wang et al. 2017 (Nat. Genet.) does not analyze deleterious load or\n“heterosis buffering.”\nThe paper focuses on the genetic basis of polyembryony in citrus using\ngenome and transcriptome data. Although some population-genomic contrasts\namong primitive, wild, and cultivated groups of genetic diversity were\npresented, the study did not quantify deleterious mutation burden nor link\nheterosis to buffering of deleterious variants.\nTherefore, as we outlined in the manuscript, the broader empirical picture remains\nunclear. The notion that occasional or residual sex may reduce mutation accumulation\nhas been proposed in theoretical models (e.g., Green & Noakes 1995, J. Theor. Biol.).\nHowever, empirical support from apomictic systems remains limited and not well-\nestablished. The only case study to date, in the R. auricomus complex (Hodač et al.\n2019, BMC Evol. Biol.), genotyped three F progeny arrays using six microsatellite\nmarkers and combined these data with mathematical modeling in a polyploid apomictic\ncontext—without genome-wide validation. As highlighted by Hörandl et al. (2020, Crit.\nRev. Plant Sci., “Genome Evolution of Asexual Organisms and the Paradox of Sex in\nEukaryotes”), the interplay between facultative sexuality, selfing, and apomixis\nintroduces further complexity into mutation dynamics in asexual plants—an area that\nremains poorly understood, largely due to the scarcity of genomic resources.\nBy contrast, our work provides—to our knowledge—the first genome-scale\nempirical validation in apomicts that rare recombination-mediated LOH occurs\nand shapes deleterious variation. We analyze a diploid, woody, wild system (Carya)\nthat has been largely overlooked relative to previously studied apomictic systems,\nwhich are predominantly polyploid herbaceous plants or crops. Specifically, our study\nfocuses on two non-hybrid apomicts (C. cathayensis and C. dabieshanensis), a sexual\nspecies (C. tonkinensis), and their hybrid derivative, the apomict C. hunanensis. Our\ndataset comprises a new phased reference genome for C. hunanensis, 195 resequenced\nadult trees and 180 mature embryos (60 per apomictic species). We validate\nreproductive modes and conduct multi-level comparisons (sexual vs apomicts;\nhaplotype contrasts within the hybrid apomict; apomictic parents compared to\neach other; apomictic vs sexual embryos). Critically, we introduce genome-wide\nLOH analysis in apomicts and show that embryos from apomictic mothers sometimes\nundergo recombination-mediated LOH that unmasks recessive deleterious alleles (a\ngenetic cost), while adult relatedness patterns indicate that occasional sex can generate\nnovel genotypes that persist. Remarkably, despite harboring more heterozygous\ndeleterious variants, apomictic adults showed lower realized mutation loads than their\nclosely related sexual species. Thus, our study does not simply transpose an old concept\nonto a new taxon; rather, it provides the first genome-wide evidence in a diploid\nwoody apomictic complex for how rare recombination influences mutation dynamics,\nwith clear ecological and agricultural relevance.\n2.Methodological concerns undermine key conclusions.\n(a) LOH Definition and Validation: The classification of LOH as 50-kb windows with\nheterozygosity <0.001 is arbitrary and lacks justification. This threshold is not\ncompared to sexual controls (e.g., C. tonkinensis) or validated against independent\nmethods (e.g., allele frequency spectra), raising doubts about whether reported LOH\nregions truly reflect recombination or technical artifacts.\nAnswer: The issues regarding the LOH threshold, the use of sexual controls, and the\npossibility of technical artifacts were carefully considered in our analyses. To ensure\nclarity, we have now substantially expanded the description of these safeguards in the\nrevised manuscript (Lines 628–654).\nTo define LOH, we used a cutoff of heterozygosity <0.001, which was not arbitrary\nbut benchmarked against the sexual control C. tonkinensis, as the reviewer rightly\nsuggested. In this species, only 2.01% of ~346k 50-kb windows fell below this value,\nconfirming that 0.001 lies in the extreme left tail of the sexual heterozygosity\ndistribution (Fig. S15). A genome-wide heterozygosity plot (Fig. S16) further illustrates\na generally stable background level punctuated by rare LOH segments where\nheterozygosity approaches zero. To test robustness, we also repeated the analyses using\na more stringent threshold of heterozygosity <0.0005 (met by only 0.67% of windows\nin the sexual control) and obtained consistent LOH patterns (Figs. S20 and S23).\nTo minimize potential technical noise, heterozygosity was computed from per-\nbase gVCF information with stringent filters: sites outside individual-specific depth\nranges, indels/multiallelics, and heterozygotes failing allele-balance filters were\nexcluded, while windows with >50% missing sites were reported as NA. These\nsafeguards support the use of 0.001 as a conservative and biologically justified criterion\nfor LOH detection. In the revised manuscript, the 0.001-based results are presented in\nthe main text, with the 0.0005 cutoff included as supporting evidence (Lines 213–232).\nFigure S15A. Distribution of window heterozygosity across sexual control Carya\ntonkinensis.\nHistograms show heterozygosity values calculated in 50-kb windows for 36\nresequenced adult individuals of C. tonkinensis. The dashed vertical line marks the\npredefined cutoff of 0.001 used to define loss of heterozygosity (LOH). Only 2.01% of\n345,979 windows fall below this value, confirming that it represents the extreme left\ntail of the sexual heterozygosity distribution and provides a conservative threshold for\nLOH detection. More details can be found in the Methods.\n(b) Sample Size Limitations: The inference of \"residual sexual reproduction\" relies on\nextremely small numbers of embryos with LOH (8 in C. cathayensis, 11 in C.\ndabieshanensis, 3 in C. hunanensis). This paucity of data makes it impossible to\nrobustly quantify recombination rates or assess its evolutionary significance, weakening\nclaims about \"costs and benefits of residual sex.\"\nAnswer: Our study does not attempt to estimate recombination rates per se; rather, to\nestablish the occurrence of recombination-mediated LOH in apomictic embryos and to\nevaluate its genetic consequences. Although the numbers of sexually derived embryos\nare modest (8, 11, and 3 out of 60 sampled per species), their detection across the three\napomictic hickory species provides clear evidence of sexual reproduction–mediated\nrecombination. In fact, a single embryo with LOH would have been sufficient to show\nthat recombination occurs. In addition, since an individual apomictic adult hickory tree\ntypically produces hundreds to thousands of fruits each year, the absolute number of\nfruits generated through sexual reproduction is still considerable. Therefore, once their\noccurrence is established, it allows us to reveal whole-genome patterns of LOH and\nmutation load in these 22 (3+8+11) sexually produced embryos. These patterns directly\ninform the genetic costs and potential benefits of residual sex, thereby making it\nmeaningful to explore their evolutionary significance.\nTo avoid over-interpretation, we have revised the text to remove strong qualifiers\nsuch as high and adjusted the Abstract (Lines 27–31) to read: “Remarkably, rare\nembryos from apomicts underwent recombination-mediated LOH, exposing deleterious\nmutations to selection. These findings reveal the genetic cost of residual sex, while also\nindicating its role in generating novel genotypes, supported by close relatedness among\nadult apomicts.” We have also revised the relevant sections of the Results and\nDiscussion for clarity.\n(c) Mutation Load Calculation: The use of SIFT scores to classify deleterious variants\nwithout accounting for lineage-specific functional constraints (e.g., divergence between\nhickory haplotypes) may misclassify neutral or adaptive mutations as deleterious,\nskewing comparisons between apomictic and sexual lineages.\nAnswer: We thank the reviewer for this constructive point. We acknowledge that any\nsingle predictor may, in principle, misclassify some missense variants, particularly in\ncross-species contexts. Nevertheless, SIFT4G (Vaser et al. 2016, Nat. Protoc.) has been\nwidely applied in comparable studies of non-model trees and other plants using a single\nreference genome—for example, in Ostrya rehderiana and congeners (Yang et al. 2018,\nNat. Commun.) and in multiple Populus species (Liu et al. 2022, Mol. Biol. Evol.).\nTo further assess robustness, we evaluated our data using two complementary\napproaches: (i) applying an alternative predictor, Grantham scores (Grantham 1974,\nScience), to the same dataset mapped to the C. hunanensis Haplotype E genome\n(derived from apomictic lineages), and (ii) repeating SIFT4G analyses with the\nalternative haplotype, Haplotype W (derived from the sexual lineage). Both Grantham-\nbased and haplotype-W–based predictions closely mirrored the patterns obtained with\nSIFT4G on Haplotype E (Figs. S4D, S29–S30), with the key conclusions unchanged:\napomictic adults harbor more heterozygous deleterious variants but lower realized\n(homozygous) loads than the sexual relative C. tonkinensis. Therefore, we consider our\nresults robust to both predictor choice and divergence between hickory haplotypes.\nAccordingly, we present SIFT4G with Haplotype E as our primary predictor and\nreference, and Grantham scores together with Haplotype W as orthogonal checks for\ncross-species comparisons, with revisions added in the main text (Lines 252–274) and\nSupplementary Information (Figs. S29–S30). For within-species analyses (haplotype-\nlevel contrasts in C. hunanensis and embryo data from the three apomictic species),\nlineage divergence is not a confounding factor, so SIFT4G results based on Haplotype\nE were retained.\n3.Inconsistent Interpretation of Results.\n(a) It claims apomicts have \"lower realized mutation loads\" but acknowledges they\nharbor \"more heterozygous deleterious variants.\" This tension is attributed to \"hybrid\nbuffering\" in C. hunanensis, but no mechanistic evidence (e.g., gene expression data\nshowing masked recessive effects) supports this assertion.\nAnswer: We agree that our study does not provide direct mechanistic evidence (e.g.,\ngene expression data) for the masking of recessive deleterious alleles in the hybrid\napomictic species C. hunanensis. We have therefore revised the Discussion (Lines 381–\n387) to clarify this, adding the following text:\n“Among them, the hybrid apomict C. hunanensis is particularly notable, showing\ngenomic patterns consistent with expectations under the hybrid buffering effect64,\nwhereby admixture between highly divergent parental genomes promotes fixed\nheterozygosity 65. While our observations align with this theoretical framework, we note\nthat direct functional validation (e.g., gene expression assays) will be required to\nconfirm the underlying mechanism.”\n(b) LOH is proposed to \"expose deleterious mutations to selection,\" yet sexual embryos\nwith LOH show higher homozygous deleterious loads. The suggestion that \"selection\neliminates harmful alleles early\" is speculative and unsupported by data on embryo\nviability or fitness.\nAnswer: We thank the reviewer for this comment and acknowledge that our original\nwording was unclear. We have revised the sentence (Lines 410–417) to read:\n“While some sexually derived mature embryos showed elevated mutation loads,\nnewly formed LOH regions in two such embryos—one from C. cathayensis and one\nfrom C. dabieshanensis—contained no homozygous LOF mutations. This pattern is\nsuggestive of selection acting at multiple stages, with some genotypes carrying highly\ndeleterious alleles possibly purged early, whereas other sexual embryos still retained\nhomozygous LOF mutations, indicating that selection may also act later in development\nor during seedling establishment.”\nWe agree that this remains indirect evidence. Direct confirmation would require\nanalyses of fruits and seedlings at different developmental stages, which would provide\na clearer picture of how LOH and deleterious alleles interact during embryogenesis.\nReviewer #2 (Remarks to the Author):\nThe manuscript uses advanced kinship analysis of Carya species with dense SNP data\nand the KING algorithm to study clonality and possible apomixis. The authors apply\nthe KING algorithm, which allows them to analyze unfiltered SNP datasets and look\nfor duplicated (clonal), first, second, and third-degree kin relationships. The study\nclaims that there is a lot of clonality, sometimes almost complete, within these species.\nBy adding resequencing data from adult tissues and sequencing 180 embryos, the\nauthors try to answer old questions about hybridization, persistent lineages, and the\neffect of apomictic versus sexual reproduction in this complicated hickory group. The\nstudy is relevant, the methodology is solid, and the evolutionary and ecological\nquestions are important. Still, some points in the manuscript need to be explained better\nto make the conclusions clearer and reliable. Studying clonality on a genome-wide level\nin woody apomictic plants is new and important, including for the ecology and\nagriculture of hickory. The use of site frequency spectra (SFS) and loss-of-\nheterozygosity (LOH) analysis makes the argument for clonality even stronger and\ngives different views on the genetic structure. The analysis shows clonality at the\npopulation level very well. Including 180 mature embryos from three apomictic species\nis a very valuable part of the study.\nAnswer: We sincerely thank the reviewer for the positive evaluation of our study,\nparticularly for recognizing the novelty and importance of applying genome-wide\nkinship, SFS, and LOH analyses in a clonal or asexual system, and for highlighting our\nefforts to address longstanding questions about the genomic consequences of rare\nsexual reproduction for lineage persistence in apomicts.\nThe authors say that rare loss-of-heterozygosity (LOH) events in embryos are probably\ncaused by leftover recombination from occasional sexual reproduction, not from\nmitotic mistakes. This explanation is reasonable and interesting, and it fits with the\nfinding that there is no genome-wide LOH. However, it would be helpful to better\nexplain the difference between LOH caused by mitosis and LOH caused by automixis.\nFor example, this could be done by showing allele imbalance of LOH patterns on\nchromosomes.\nAnswer: We thank the reviewer for this constructive suggestion. We agree that\nclarifying the mechanistic differences between LOH generated by mitotic repair,\nautomixis, and residual sexual recombination improves interpretation.\nMitotic repair of DNA double-strand breaks typically generates very short and\nhighly localized LOH tracts through gene conversion or break-induced replication—\nusually hundreds of bp to a few kb (Symington et al., 2014, Genetics). Such events are\nsporadic and locus-restricted, and LOH tracts are unlikely to be generated across\nmultiple chromosomes.\nAutomixis, by contrast, produces genome-wide LOH because nearly all\nchromosomes are affected. Crucially, these tracts follow structured and predictable\npatterns: terminal fusion yields LOH around centromeres, while central fusion\ngenerates LOH at chromosome ends while retaining heterozygosity near centromeres\n(Lenormand et al., 2016, Genetics). If automixis were the prevailing mechanism, such\ncentromere- or telomere-anchored biases would recur consistently across embryos.\nHowever, our observed LOH tracts occur at variable chromosomal positions and do not\nmatch the expected automixis signatures.\nResidual sexual recombination offers the best explanation: meiotic crossing-over\nduring rare sexual events produces stochastic, embryo-specific LOH tracts that vary in\nlocation among individuals. This scenario fits our data from three apomictic Carya\nspecies, where 8/60, 11/60, and 3/60 embryos carried large LOH tracts.\nAs you suggested, we have added a schematic illustrating the characteristic LOH\nsignatures expected under mitotic repair, automixis, and residual sexual recombination\n(Fig. S14), and we have revised the corresponding Results and Methods sections in the\nmain text accordingly (Lines 200–232, 655–673).\nFigure S14. Schematic comparison of loss-of-heterozygosity (LOH) signatures\nexpected under different mechanisms.\nA Mitotic repair typically produces very short, localized LOH tracts (blue) arising from\ngene conversion or break-induced replication, usually restricted to a few hundred bp to\nkb. B Automixis generates genome-wide LOH with structured patterns: terminal fusion\nyields LOH around centromeres, whereas central fusion produces LOH at chromosome\nends while retaining heterozygosity near centromeres. C Residual sexual recombination\nresults in stochastic LOH tracts of variable chromosomal positions across embryos,\nreflecting random meiotic crossing-over. Blue regions indicate LOH, grey regions\nindicate heterozygous segments.\nMy main doubts come from the fact that Carya is a sexually outcrossing (allogamous)\nspecies. Generally, in diploid species, the level of heterozygosity is roughly halved\nevery generation when selfing occurs. In contrast, cross-pollination keeps or even\nincreases heterozygosity in the population, preserving genetic variation. In outcrossers,\nselfing causes higher LOH with a fast increase of homozygosity and a decrease of\nheterozygosity. Cross-pollination results in less LOH and keeps heterozygosity and\ngenetic diversity high. This is why allogamous species that do repeated selfing\nexperience a quick loss of heterozygosity and often suffer inbreeding depression. My\npersonal idea to explain the results (embryos from sexual reproduction showing LOH),\ntogether with the authors’ conclusion, is that in the crosses, both pollen and egg come\nfrom two plants of the same clone. But this should be studied more carefully. Still,\nseeing such LOH in only one generation (embryos from a mother plant) with crossing\nis unexpected.\nAnswer: We fully agree with the reviewer’s point. In our Carya system, the sexual\nembryos showing LOH are most plausibly explained by selfing or by intra-clone/close-\nkin inbreeding during rare sexual events. This interpretation is supported by (i) an ~50%\nreduction in heterozygosity relative to apomictic mothers (Fig. S31), and (ii) the\npresence of first–third-degree relatedness among adults, which is compatible with\nmating within clonal lineages as well as inbreeding among closely related genets under\nthe wind-pollinated condition. We have clarified this interpretation in the manuscript\nand revised the Results accordingly (Lines 208–213).\nThe authors should maybe sequence three biological replicates only for the samples that\nthey think come from sexual reproduction, just to confirm that the LOH is consistent.\nAnswer: We agree that biological replication would be the most direct way to confirm\nthat LOH patterns in sexual embryos are consistent and not sequencing artifacts.\nHowever, as the embryo sequencing was performed in 2023 and DNA extracts were not\npreserved at sufficient quality, additional replicates cannot be generated for this dataset.\nImportantly, several independent lines of evidence argue against artifacts. All\napomictic embryos showed fully consistent, heterozygous genomes, whereas the\nputative sexual embryos displayed large-scale LOH across nearly every chromosome.\nGiven our stringent quality controls, these patterns are unlikely to result from\nsequencing errors or misidentification of LOH.\nI also suggest explaining this part more clearly. This is important because the embryo\ngenotypes, especially related to LOH and recombination, are a large part of the study’s\nanalysis and interpretation. So, this section should be better described better in the main\nmethods or detailed more clearly in the Supplementary Information. Also, I suggest\nadding more discussion about other possible reasons, like mitotic recombination or\ntechnical problems. In conclusion, this paper gives valuable information about\napomixis, clonality, and how hybrid species in long-living trees stay stable through time.\nIt just needs some data/revisions, especially to explain the embryo part better.\nAnswer: We are grateful for the reviewer’s concluding remarks highlighting the value\nof our study. We also thank the reviewer for the helpful suggestion regarding the LOH\nanalyses in sexual embryos. As noted above, we have clarified and elaborated these\nanalyses, providing more detailed descriptions of the procedures used to minimize\npotential technical artifacts together with a clearer comparison of possible mechanisms\n(including mitotic and automixis recombination). These revisions have been\nincorporated into the main text and the Supplementary Information."}, {"round": 2, "reviewer_comments": {"reviewer_1": "This study presents the first genome-wide evidence of recombination-mediated loss of heterozygosity (LOH) in a diploid woody apomictic complex (Carya), demonstrating how rare sexual events unmask deleterious alleles while contributing to genotypic diversity. The authors show that apomicts maintain lower realized mutation loads than sexual relatives, despite higher heterozygosity for deleterious variants, suggesting hybrid buffering and purging. Methodological concerns about LOH detection and mutation load estimation were robustly addressed through sexual control benchmarking, sensitivity analyses, and orthogonal validation. While the number of recombinant embryos is limited, their detection across three species confirms residual sexuality. The study fills a notable gap in apomixis research, which has largely focused on polyploid herbs. The manuscript is now substantially strengthened through revisions and supplementary analyses, and I support its acceptance pending minor textual clarifications.", "reviewer_2": "Authors have answered to most of my suggestions and the manuscript is clearer now almost all my comments were addressed."}, "author_response": ""}]

Decision letter