Long-Read Sequencing Reveals Pre-Meiotic Gene Conversion in Sperm Cells

Research Ethics and Sample Governance

Research involving TwinsUK samples was conducted with informed consent under the TwinsUK BioBank ethics framework. The study was approved by the North West–Liverpool Central Research Ethics Committee (REC reference 19/NW/0187; IRAS ID 258513), as well as by earlier approvals granted to TwinsUK by the St Thomas’ Hospital Research Ethics Committee and subsequently the London–Westminster Research Ethics Committee (REC reference EC04/015). Sperm samples from SL were provided by Discovery Life Sciences as discarded medical waste under IRB protocol DLS-BB050. These samples were collected during standard-of-care testing requested by the patient’s physician, with excess material subsequently released for research.

Sperm Sample Extraction and PacBio Sequencing

We analysed 15 sperm samples from 13 donors of European ancestry. Eleven donors provided one sample, while two donors provided two samples collected at different ages (AA1-s1 and AA1-s2; AN-s1 and AN-s2). Donor ages at sample collection ranged from 24 to 74 years (Supplementary Table 1). Samples AA2-t1 and AA2-t2 were obtained from monozygotic twins. The TwinsUK cohort includes more than 14,000 volunteers, predominantly middle-aged female individuals, who have participated in a longitudinal study for approximately 30 years. The cohort includes lifestyle and health questionnaires, biomedical measurements, biological sample collection and multi-omics profiling—including genetics, metagenomics and metabolomics—across multiple visits.

For TwinsUK samples, we extracted high-molecular-weight (HMW) DNA from bulk sperm using the Circulomics Nanobind Tissue Kit (102-302-100) and the NEB Monarch HMW DNA Extraction Kit for Cells & Blood (T3050) UHMW protocol. We modified the protocol to account for the tighter packing of sperm chromatin. The lysis solution was supplemented with 150 mM 1,4-dithiothreitol to disrupt protamine disulfide bonds, and the proteinase K incubation was extended from 30 min to 2 h. All other steps of the Circulomics HMW DNA extraction protocol were unchanged. HMW DNA was sheared into 10–14 kb fragments using the Megaruptor 3 system (B06010003) at speed setting 30. Circular consensus sequencing (CCS) libraries were prepared using the standard CCS library preparation protocol 1.0 (100-222−300). Libraries were sequenced on Sequel IIe and Revio instruments at the Wellcome Sanger Institute using full-resolution base-quality scores.

For SL samples, we thawed specimens on ice and pelleted the cells in an isotonic sperm wash solution to remove debris and minimise somatic-cell contamination. DNA was extracted using the NEB Monarch HMW DNA Extraction Kit for Cells & Blood (T3050) UHMW protocol, with modifications to the cell-lysis step. Cells were digested at 56 °C for 1 h at 300 rpm in 100 µl Nuclei Prep Buffer, 100 µl Nuclei Lysis Buffer, 10 µl Proteinase K (NEB P8107, 20 mg ml⁻¹) and 10 µl 1 M 1,4-dithiothreitol (GoldBio in dH₂O; final concentration, approximately 50 mM). We then performed an additional 20-minute digestion with 5 µl RNase A (NEB T3018, 20 mg ml⁻¹). All remaining steps followed the original protocol. HudsonAlpha conducted HMW DNA quality control, CCS library preparation and sequencing on the PacBio Revio platform.

Blood Sample PacBio Sequencing Data Analysis

Raw PacBio CCS sequencing data and genome assemblies were obtained from the Platinum Pedigree dataset. We excluded flow cells containing multiple samples because demultiplexing errors caused cross-sample contamination. The four first-generation samples were also excluded because they originated from cell lines. The final analysis included 12 samples.

De Novo Haplotype-Resolved Genome Assembly

We used hifiasm³⁰ (v0.19.5-r592; default parameters; HiFi-only mode) to generate a haplotype-resolved de novo assembly for each sample. The hifiasm assembler produces a partially phased assembly graph representing two haplotypes. We converted these graphs into two FASTA sequence files per sample.

Each haplotype was scaffolded against the T2T-CHM13 reference genome⁵² using RagTag⁵³ (v2.1.0; arguments “-u -w –aligner minimap2”). RagTag orients and positions contigs along the T2T reference while introducing gaps. We manually expanded gaps between adjacent contigs to 30,000N to prevent reads from mapping across a gap to two contigs. Such mappings could reflect a phase-switch error and lead to false identification of crossover (CO) events.

PRDM9 Allele Genotyping

We obtained zinc-finger array sequences for 74 previously reported PRDM9 alleles²⁷ and mapped them to every assembled haplotype from each donor. For each haplotype, we assigned the PRDM9 allele with the fewest mismatches, insertions or deletions relative to the allele sequence and assembly. All A alleles, as well as the B and D alleles, matched perfectly. We did not identify a perfect match for the newly reported PRDM9 allele in AN-s1/2. We also confirmed that the flanking sequences on both sides matched the expected regions reported previously²⁷, with up to 2 bp mismatches. Finally, we visually inspected the assembly surrounding the PRDM9 region to confirm that allele variation was not caused by assembly errors.

CCS Read Alignment and Filtering

We identified CCS reads containing residual adapter sequences using HiFiAdapterFilt⁵⁴ and excluded them from downstream analyses. Reads were aligned separately to each haplotype with minimap2^55,56 (v2.26-r1175; arguments “-ax map-hifi –cs=short –eqx –MD”). We filtered each BAM file to retain only primary alignments using the “-F 0×900” flag in samtools view⁵⁷, and retained reads with a mapping quality of MAPQ = 60. We removed reads that mapped to only one haplotype, mapped to different chromosomes or aligned to opposite strands of the two haplotypes.

Identification of Candidate Recombinant Reads

Comparison of Haplotype Alignments

An alignment of a read to a haplotype partitions the read into segments, with each segment representing one alignment operation: (1) a ‘match’, which aligns perfectly to the reference haplotype; (2) a ‘mismatch’, which aligns but differs from the reference; (3) an ‘insertion’; (4) a ‘deletion’ relative to the reference, where the corresponding read segment has length 0; or (5) ‘soft clipping’.

We compared the two alignments of each read and generated a refined read partition based on the intersection of both alignments. Each segment in this joint partition therefore has two alignment operations, one for each haplotype. For example, the read segment [0, 40) may match both haplotypes; [40, 41) may match haplotype 1 but mismatch haplotype 2; [41, 41) may represent a deletion relative to haplotype 1 but an empty match to haplotype 2; and [41, 50) may match haplotype 1 but represent an insertion relative to haplotype 2. This joint partition formed the basis of subsequent recombinant-read analysis.

We excluded reads with more than 10 bp of soft clipping to either haplotype because these reads were frequently chimeric. We also removed reads containing more than 100 sequence errors, defined as segments—usually 1 bp—that mismatched both haplotypes.

SNP Detection and Filtering

We identified candidate single-nucleotide polymorphisms (SNPs) along each read using the joint partition. Candidate SNPs were positions that matched one haplotype but differed from the other. Many single-nucleotide mismatches occurred in regions associated with potential assembly errors, including low-complexity sequences, tandem repeats, read ends and haplotype regions with low coverage. We therefore applied multiple quality filters.

First, each candidate SNP had to be flanked by at least 10 bp of matched alignment on both sides. Second, we excluded SNPs within 1,500 bp of either read end for Sequel II data and within 400 bp for Revio data (Supplementary Methods and Supplementary Figs. 1 and 2). Third, we removed SNPs with base-quality (BQ) scores below 60 for Sequel II data, or below 40 or 50 for Revio data, depending on the largest BQ bin (Supplementary Methods and Supplementary Figs. 3 and 4). Fourth, we used sdust⁵⁸ (v0.1-r2; default parameters), an implementation of the dustmasker algorithm⁵⁹, to identify low-complexity regions in both haplotypes and removed SNPs located in either region. Fifth, we used Tandem Repeat Finder⁶⁰ (v4.09.1; arguments “2 6 6 80 10 50 500 -ngs -h -l 10”) to identify tandem repeats in both haplotypes and excluded SNPs mapping to either repeat. SNPs that passed all filters were classified as high-confidence SNPs. Reads without high-confidence SNPs were excluded from subsequent analyses.

For Sequel II data, we also identified a less stringent set of ‘classification SNPs’ using the same criteria, but with a minimum BQ score of 30 rather than 60 and a minimum distance of 500 bp from read ends rather than 1,500 bp. These SNPs were used for event classification and analysis but not for event detection. For Revio data, the corresponding thresholds were a minimum BQ score of 30 rather than 40 and a minimum distance of 200 bp rather than 400 bp.

Haplotype Assignment

We scanned each read for high-confidence SNPs and recorded the haplotype matched by each SNP. More than 99.9% of reads contained SNPs that were fully consistent with at least one haplotype and were therefore not candidate recombinant reads. These reads were nevertheless useful for identifying and filtering assembly errors. For each read, we recorded the percentage of SNPs consistent with haplotype 1 and haplotype 2.

Additional SNP Coverage Filters

Assembled haplotype regions with insufficient coverage from consistent reads were more susceptible to errors and false-positive calls. We therefore required every SNP to be overlapped by at least three reads containing more than 95% SNPs consistent with haplotype 1, and by at least three reads meeting the same criterion for haplotype 2.

Candidate Recombinant SNPs and Reads

After coverage filtering, we recalculated for every read the percentage of high-quality SNPs consistent with each haplotype. Reads whose SNPs were not fully consistent with either haplotype were retained as candidate recombinant reads.

Classification of Candidate Recombinant Reads

Additional Read-Level Filtering

We applied additional filters to classified reads and their annotations to reduce false-positive recombinant events. A transition pair was defined as two adjacent SNPs on a read that mapped to different haplotypes. Every candidate recombinant read contained at least one transition pair. We considered two transition pairs from different reads that mapped to identical haplotype coordinates as potential phasing or assembly errors, because each recombination event was assumed to occur in exactly one read. Reads containing transition pairs observed in other reads were therefore discarded. Given the sequencing depth of approximately 20×–160× (mean, 78×), the probability of observing the same recombinant molecule more than once by chance was low: 4.4 × 10⁻⁷ to 2 × 10⁻⁶ for CO events and 1.6 × 10⁻⁶ to 7.8 × 10⁻⁴ for non-crossover (NCO) events. Repeated events were therefore more likely to represent technical artefacts.

We also discarded reads if the genomic interval between a transition pair was not covered by at least three reads in either haplotype, as described under ‘Further SNP filtering measures’. Low-coverage regions may indicate a phasing error in the haplotype assembly and could produce false CO calls.

Finally, we removed reads suspected of cross-sample contamination. A read was classified as contaminated if it contained a single-base mismatch to both haplotypes at a site of known population variation, suggesting that it originated from another genome. We first aligned all reads to the GRCh38 reference genome and compared 1 bp mismatch positions with known population-variation sites from the 1000 Genomes Project⁶¹. We then aligned all previously defined high-quality SNPs to GRCh38 to generate a dataset-specific genetic-variation callset and removed reads containing mismatches at these variable sites.

Recombination Event Classification

We classified candidate reads as CO when they contained one transition pair flanked by at least one SNP on each side—for example, ‘1122’, representing two SNPs from haplotype 1 followed by two SNPs from haplotype 2. Reads were classified as NCO when they contained two transitions, such as ‘121’. Reads with one transition but without flanking SNPs on at least one side—for example, ‘12’, ‘112’ or ‘122’—were classified as ‘ambiguous’ because CO and NCO events could not be distinguished. Reads containing more than two switches were classified as ‘complex’.

For CO reads, we calculated the mean CO recombination rate in centimorgans per megabase (cM Mb⁻¹) across the interval between the two SNPs spanning the CO event (‘12’ or ‘21’). Methods for comparing CO and NCO genetic-length distributions are provided in the Supplementary Methods and Supplementary Fig. 5.

Probability of Detecting the Same NCO Event More Than Once

We modelled recombinant events as a Poisson process to estimate the probability of observing the same gene-conversion event more than once. Approximately 30 to 50 CO events and 200 to 1,000 NCO events occur per meiosis, with 25,000 to 50,000 meiotic recombination hotspots across the genome. The expected CO and NCO event rates per hotspot per meiosis (λ) were therefore 0.0006–0.002 and 0.004–0.04, respectively. Under a Poisson distribution, the probability of observing two or more CO or NCO events in the same hotspot was calculated as P(Poisson X ≥ 2) = 1 − e^−λ (1 + λ).

Read Mapping to Reference Genomes

For genetic-distance analyses, we used the European male-specific refined genetic map³², provided in GRCh37 coordinates. We aligned reads to the GRCh37 reference genome⁶² with minimap2^55,56 (v2.26-r1175; arguments “-ax map-hifi –cs=short –eqx –MD”) to analyse CO-rate information. The DSB map^37,38 was provided in GRCh38 coordinates, so reads were similarly aligned to the GRCh38 reference genome⁶³ for DSB analyses. To calculate distances from telomeres, we also aligned reads to the T2T reference genome⁵².

Positional Bias of Converted Markers Relative to the DSB–PRDM9 Motif

To evaluate positional bias relative to the DSB–PRDM9 motif, we calculated the distance between converted markers in NCO reads overlapping a DSB–PRDM9 motif and the motif centre, accounting for the strand on which the motif occurred. We tested for asymmetry using a permutation test. First, we calculated the fraction of reads with a negative distance. We then randomly reversed the sign of each distance, equivalent to randomising motif strands, and recalculated the fraction of negative distances 10,000 times. The permutation test P value was the fraction of randomisations in which the negative-distance fraction was equal to or greater than the observed fraction.

PRDM9 Motif Detection in Recombinant Reads

We detected PRDM9 motifs using FIMO⁶⁴ with the PRDM9 A/A position weight matrix provided by Hinch et al.³⁸. FIMO was run on recombinant reads and on a control set of 100,000 randomly selected CCS reads. Only motif matches with a q value < 0.01 were retained for downstream analysis.

GC-Biased Gene Conversion Analysis

To calculate GC-biased gene conversion (gBGC) among converted SNPs in NCO reads, we included only SNPs in which one allele was G or C and the other was A or T. gBGC was calculated as the proportion of converted SNPs carrying G or C within this subset.

Inference of Non-Crossover Tract-Length Distributions

We inferred the NCO tract-length distribution using an approximate Bayesian computation-like approach. Tracts were modelled either with a single geometric distribution having mean length L, or with a mixture of two geometric distributions with proportions m:1 − m and mean lengths L₁ and L₂.

The fitting procedure used three statistics derived from detected NCOs: (1) the number of converted SNVs; (2) the distance between the first and last converted SNVs, representing the lower bound; and (3) the distance between the immediately flanking SNVs upstream and downstream of the conversion, representing the upper bound.

For each proposed parameter set—L for the single-geometric model, or m, L₁ and L₂ for the mixture model—we simulated the expected distributions of these statistics. We randomly selected 10.7 million reads. For each read, we simulated an NCO tract by (1) drawing a random start position along the read; (2) selecting a geometric component with probability m or 1 − m; and (3) sampling a tract length from the selected distribution. We then determined which SNVs were converted, assessed whether the event was detectable and recorded the number of converted SNVs, lower bound and upper bound. Undetectable NCO events were excluded.

We compared simulated and observed distributions by calculating the two-sample Kolmogorov–Smirnov statistic for each of the three summary statistics and summing the results. Parameters were selected by minimising this summed statistic. For the mixture model, we initially performed a coarse grid search using 11 values of m (0.97–1), 11 values of L₁ (10–100 bp) and 21 values of L₂ (200–500 bp). We then refined the search using narrower ranges: 7 values of m (0.988–0.994), 31 values of L₁ (20–50 bp) and 31 values of L₂ (650–1,550 bp). For the single-geometric model, we fixed m = 1 and evaluated 500 values of L between 10 and 1,000 bp. Confidence intervals were estimated using 100 nonparametric bootstrap samples, repeating the complete fitting procedure for each sample. For each parameter, the confidence interval was defined by the 2.5th and 97.5th percentiles of the bootstrap estimates.

Statistical Tests for Equality of Distributions

To test whether two samples had equal distributions, we used the two-sample Anderson–Darling statistic⁶⁵. This statistic is based on the integrated squared difference between the empirical cumulative distribution functions of two measurement sets, such as recombination rates. We assessed significance using permutation testing. Measurements were randomly assigned to the two samples while preserving sample sizes, and the Anderson–Darling statistic was recalculated for 10,000 permutations. The P value was the fraction of permutations with an Anderson–Darling statistic greater than the observed value.

To compare one focal sample with multiple other samples, we calculated the Anderson–Darling statistic between the focal sample and each comparison sample and summed the resulting statistics. For each of 10,000 iterations, measurements were independently permuted between the focal sample and each comparison sample while preserving the relevant sample sizes. We then calculated the summed Anderson–Darling statistic. The P value was the fraction of paired permutation sets producing a statistic greater than the observed statistic.

Age-Related Effects in Blood Samples

To test whether recombination-event rates were associated with sample age, we modelled the number of events per sample using a Poisson distribution, with the expected value proportional to sequencing coverage. We compared a null model, in which the event rate λ was constant across samples and the expected count was λ × coverage, with an alternative model in which the rate increased linearly with age and the expected count was λ × coverage × sample_age. We optimised the Poisson log likelihood for each model and calculated a likelihood-ratio statistic to assess whether age improved model fit. Statistical significance was evaluated by randomly permuting sample ages among individuals and refitting the age-dependent Poisson model. The permutation P value was the fraction of permutations in which the fitted slope exceeded the slope observed in the original data.

Estimation of Non-Meiotic NCOs

We estimated the proportion of non-meiotic NCOs using three complementary approaches. First, we modelled the genetic-length distribution of NCOs (Fig. 2b) as a mixture of the genetic-length distributions of all reads and COs. We maximised the likelihood of the observed NCO lengths on a logarithmic scale under a p:(1 − p) mixture of these distributions. Each distribution was smoothed using kernel density estimation to enable continuous probability-density evaluation. The maximum-likelihood estimate suggested that 45% of NCOs followed the non-meiotic length distribution, with a 95% confidence interval of 39.4–51.2%. However, the fit was imperfect, potentially because of model misspecification even if all NCOs originated during meiosis. This may also reflect variable resolution of DSBs into COs and NCOs.

Second, because 45% of NCO reads overlapped a DSB–PRDM9 motif, compared with 67% of COs and 13.5% of all reads, we estimated that approximately 41% of NCO reads were non-meiotic. Third, using GC bias, we estimated the non-meiotic fraction at approximately 32%. Overall NCOs showed a GC-conversion rate of 60.7%, compared with 65.7% among NCO reads overlapping a PRDM9 motif—representing putative meiotic NCOs—and a 50% background rate. This approach may be affected by classification errors because NCOs may not be accurately assigned as meiotic or non-meiotic.

Reporting Summary

Additional information about the research design is available in the Nature Portfolio Reporting Summary linked to this article.