Increasing the number of molecules per library

To test whether the number of (recombinant) molecules per library would be increased by increasing the DNA loading, we generated two additional 10X libraries from the same DNA pool by loading 0.40ng (P50L40) and 0.75ng (P50L75) into the Chromium Controller (Fig. 1; Methods). We sequenced these libraries with 104.8 and 212.3 million read pairs in anticipation of an increased number of molecules within both libraries (Fig. 2b–d). Following the same analysis as for the first library, the higher DNA loading drastically increased the number of recovered molecules. As compared to the 2.7 million molecules for P50L25, we now found 4.7 and 10.5 million molecules for P50L40 and P50L75 after filtering (Table 1). This also increased the number of recombinant molecules to 1012 in P50L40 and 2519 in P50L75 as compared to 558 in P50L25. This, however, had the cost of also increasing the FP rate by 5% in P50L40 and 10% in P50L75 (Table 1).

To investigate the effect of the FPs, we compared the distribution (i.e., recombination landscape) of all 2519 COs in the library with the highest FP rate, P50L75, with two other landscapes: one calculated only from 1874 true COs (those overlapping with the 400 benchmark COs) and one calculated only from the remaining 645 false recombinant molecules (Methods). The recombination landscape calculated for all COs (TP+FP, i.e., the set of TP and FP COs) was nearly identical to the landscape generated from true COs only (TP), while the distribution of the FPs was highly random (Fig. 3a; Supplementary Fig. ²: K–S test, p-value=9.6e-01), suggesting that FPs hardly obscure true recombination landscapes. In fact, correlating the sliding window-values of each of the three CO landscapes revealed that the landscape calculated from all CO sites was almost perfectly correlated to the one generated from real COs (Fig. 3b: Correlation test, Pearson’s r 0.97, p-value<2.2e-16), while the frequency of FP along the chromosomes was not correlated to the frequency of real COs and was only marginally correlated to landscape of all COs.

Open in a separate window

Fig. 3

Genome-wide CO landscape formation and feature association. a Sliding window-based (window size 1Mb, step size 50kb) recombination landscapes calculated for true positives (TP), false positives (FP) and the combined set of TP+FP recombinant molecules of P50L75. Heterochromatic regions⁴⁴ are indicated by rectangles in light blue. b Correlation tests of regional CO frequencies comparing TP vs. TP+FP (Pearson’s r 0.97), FP vs. TP+FP (Pearson’s r 0.25) and FP vs. TP (Pearson’s r 0.05). c Association of the TP and TP+FP sets with different genomic features in contrast to the random expectation (permutation test). d. Correlation tests of regional GC-content with CO frequency within the TP and TP+FP sets. e Association of TP or TP+FP sets with DNA methylation, compared with a random expectation (t-test; Methods). Note: in c and e, the expectations were obtained based on randomly sampled intervals within the reference genome excluding heterochromatic regions. The middle position of a CO interval was used for associations with genomic features in c, where the promoter of a gene was defined as the 1000bp region upstream of the transcription start site, and the gene start/end as the first/last 200bp of a gene. For each feature, a permutation test (oneway test in R coin package) was performed between the TP set (or TP+FP set) with 1,000 sets of randomly sampled 3000 intervals. The asterisks indicate the observed values (either from the TP or the TP+FP set) were significantly different from a random expectation. Source Data are provided as a Source Data file Source_Data_main_Figure_³.zip

Association of COs with genomic features

To test the association of COs with genomic features, we checked all of the CO sites for their annotation in the Col-0 reference sequence (Fig. 3c). In comparison to randomly placed CO sites, permutation tests showed that the 1874 true COs sites were significantly enriched in promoter regions (p-value 5.0e-03) and intergenic regions (p-value 1.4e-03) and were significantly depleted in gene bodies (p-value 6.0e-04) and transposable elements (p-value 9.0e-04) recapitulating regional preferences, which have been described before^40,41. When additionally including the FP, the CO set showed the same significant regional associations with one exception: COs were slightly, but significantly enriched compared to random COs at gene ends. Gene ends represented only a minor fraction of all COs in both datasets and the difference between the two was marginal (5.16% within TP+FP CO sites vs. 4.90% within TP COs sites) and these percentages were much closer to each other than the percentage in the randomized CO dataset (4.40%). A later analysis with a larger CO dataset helped to avoid this spurious enrichment (see next section). To examine CO associations with genomic features at a finer scale, we analyzed the regional preferences of the CO sites in individual transposable element super families because these were found to have strong differences in CO rates⁴². Though this only included 331 putative CO sites, we found that COs were significantly depleted within LTR-Gypsy and En-Spm super families and enriched within Helitron and LINE elements, as has been previously shown⁴², whether or not FP were included (Supplementary Fig. ³). In addition, the presence of FP did not affect the previously reported relationships between COs and GC-content or DNA methylation⁴⁰. The two datasets showed nearly identical correlations between CO frequency and GC-content (Fig. 3d: Correlation test, Pearson’s r −0.46 and −0.45, both p-value<2.2e-16). Similarly, COs in both datasets were found in regions with low levels of methylation (Fig. 3e).

Together this suggests that our method is not only effective when facing an increasing amount of FP in libraries with a large number of molecules, but that it is also powerful enough to accurately identify chromosome-wide and local CO patterns.

Increasing the number of genomes per library

Many of the recombinant molecules that we found in the pool had identical CO breakpoints. This implied that we re-discovered some of the CO breakpoints in independent molecules, suggesting that there were more recombinant molecules in the library than independent CO breakpoints in the pooled genomes. For instance, there were only 363 distinct COs recovered by the 1874 recombinant molecules in P50L75 (Table 1). Even though identifying a single CO breakpoint multiple times can help to increase its resolution and reliability, it does reduce the overall number of distinct COs that can be found with one library.

The number of recombinant molecules that can be identified within one 10X library greatly depends on the number of molecules in the library. As the probability to overlap with a CO breakpoint is theoretically the same for each molecule, the number of CO molecules increases linearly with the number of analyzed molecules. The number of molecules heavily relies on the amount of DNA that is loaded in the 10X Chromium machine (while sequencing depth has only a marginal effect on molecule recovery). As a consequence, each 10X library includes an almost fixed number of recombinant molecules.

In turn, the inclusion of more plants in the pool does not increase the number of recombinant molecules. However, the inclusion of more plants does increase the number of independent recombinant molecules and thereby increases the number of distinct CO events found with one library (Supplementary Note ²). While this maximizes the number of distinct CO breakpoints that can be found with one library, it also implies that a majority of the breakpoints that are in the pool will not be found. Thus, pools with a small number of plants/CO breakpoints (smaller than the number of recombinant molecules) will reveal all CO breakpoints with high confidence; pools with a large number of plants/CO breakpoints will reveal a large number of distinct CO breakpoints, but cannot achieve the identification of all CO breakpoints. Hence, the false negative rate is mostly determined by the number of independent genomes in a pool and much less by the actual detection success of recombinant molecules.

To test the effect of genome number on the number of distinct COs in practice, we applied our method to pools of 200 and 1250 F₂ plants. We generated two individual 10X libraries each for a pool of 200 (P200R1, P200R2) and a pool of 1250 (P1250R1 and P1250R2) different F₂ plants (Fig. 1; Supplementary Figs. ¹ and ⁴; Methods). These libraries were sequenced with 168.2–194.4 million read pairs (Supplementary Data ¹) and revealed 12.3–16.0 million molecules. For P200R1 and P200R2, we found 2538 and 2334 recombinant molecules, which identified 1779 (69.5%) and 1606 (68.8%) distinct CO sites. For the pool of 1250 plants, there were 2880 and 3430 recombinant molecules, which uncovered 2386 (82.8%) and 2788 (81.3%) distinct CO sites. Comparison of these COs with 3320 COs determined using individual whole-genome sequencing for an independent set of 437 Col x Ler F₂s^20,21 showed consistent genome-wide patterns, including on identifying genomic regions with relatively high CO occurrences (Supplementary Fig. ⁵). Although we were unable to assess FP rate in these larger pools, the above consistencies in feature associations (even at gene ends with ~5% COs) suggested that FPs do not interfere with the proper identification of true CO patterns. While the larger pools revealed more distinct CO sites, there were still nearly 20% COs overlapping, but it does not necessarily imply that these are not unique CO events. As the genome size of A. thaliana is relatively small (~135 Mb⁴³), independent COs have an unneglectable probability to overlap (Supplementary Note ²; Supplementary Data ⁵). According to simulations, given a pool of 1250 F₂s, ~15% of the independent CO events are expected to overlap (Supplementary Data ⁵ and ⁶) making up for large parts of the overlapping recombinant molecules in the 1250 pools.

Estimating relative recombination frequency

Since the probability to identify a CO breakpoint within a molecule depends on the average recombination frequency among the pooled genomes, it might be possible to estimate the average recombination frequency from the fraction of observed recombinant molecules. However, as the probability to identify a CO breakpoint also depends on the length and sequencing coverage of the molecules within a library, it is only meaningful to calculate an average recombination frequency that is relative to the actual molecule characteristics (or relative recombination frequency). To test for significant differences in relative recombination frequencies between libraries, it is therefore essential that they have close-to-identical molecule characteristics. For libraries with differences in molecule length and sequencing coverage distributions, subsampling can be used to generate such identical distributions. Once the distributions are similar, the relative recombination frequency can be determined as the number of recombinant molecules per million molecules (C^M). Repeated subsampling even allows for the calculation of confidence intervals for C^M values and thereby allows for testing significant differences between the average recombination frequencies of different pools.

First, we compared the CO frequencies of P200R1 and P200R2, which were two independent libraries generated from the same DNA (Fig. 1; Fig. 2b–d; Methods). After subsampling (Supplementary Fig. ⁶), there was no significant difference in C^M values (26.7±0.4 and 26.8±0.4), as expected for these libraries (Table 2). We repeated this comparison for P1250R1 and P1250R2, which differed greatly in their original molecule characteristics (Fig. 2b–d). After subsampling to the same molecule distributions as for the smaller pools (Supplementary Fig. ⁶), the C^M values of 27.2±0.6 and 27.1±0.7 were also not significantly different from each other (Table 2). Moreover, when comparing the 200 and 1250 recombinant pools, we also found no significant difference between any of the pools (all confidence intervals overlapped), which is expected, given that all libraries were generated from individuals of the same F₂ population. Together, these results show that C^M values are stable against differences in the molecule characteristics and pool sizes and allow for determining and comparing average recombination frequencies between samples.

Table 2

Relative recombination frequency in F₂ and pollen pools

	Pool	Moleculesa	CO moleculesb	C M
F2	P200R1	4,845,588±1,556	129±7	26.7±0.4
F2	P200R2	4,856,886±1,371	130±7	26.8±0.4
F2	P1250R1	4,650,420±1,662	126±9	27.2±0.6
F2	P1250R2	4,309,936±1,454	117±10	27.1±0.7
Pollen	P8000	4,902,994±1,722	212±10	43.2±0.6

Open in a separate window

^{a, b}The values are μ±σ, where μ is the mean and σ is the standard deviation of the molecule and CO numbers in 50 random sub-samplings. For C^M, i.e., the ratio of CO molecules to total molecules scaled by a factor of 10⁶, the values are μ±s giving 95% confidence intervals for the mean μ, where s is 2.0096×σ/50^0.5 with σ being the standard deviation. Source Data are provided as a Source Data file Source_Data_main_Table_².zip

Pollen DNA extraction and library preparation

Col CEN3 qrt 420 and Ler seeds were stratified for 4–7 days at 4°C before sowing on soil in 18-pot trays and growing under 16h light, 8h dark cycles at 20°C until flowering. Pollen from a single Ler plant was used to pollinate a single Col CEN3 qrt 420 stigma to generate F₁ plants. To obtain pollen and extract DNA, we adapted a method from Drouaud and Mezard⁴⁹ as follows: Inflorescences were collected from a single F₁ plant and ground in 1mL of 10% sucrose using a mortar and pestle. 9mL of 10% sucrose were added to the mortar slurry and the resulting homogenate was mixed by pipetting up and down with a wide bore 1mL pipet tip and filtered through an 80-μM nylon mesh before centrifugation at 350×g for 10min at 4°C. The supernatant was discarded and the pollen pellet was washed two times with 10% sucrose. The pellet was resuspended in four volumes of lysis buffer (100mM NaCl, 50mM Tris–HCl (pH 8)), 1mM EDTA, 1% SDS and proteinase K was added to achieve a concentration of 20μg/mL. Five to ten 2-mm glass beads were added to the sample and it was vortexed at full speed for 30s. To check for pollen disruption, a 1-μL sample was removed and mixed with 10μL of lysis buffer and examined under a microscope. During this time, we verified that the pollen sample was free of large amounts of cell debris. The sample was then vortexed for 30s and checked for pollen disruption. An equal volume of Tris-saturated phenol was added and the sample was placed on a rotating wheel at room temperature for 30min. After centrifuging at 15,000×g for 10min, the supernatant was transferred to a new tube and mixed with an equal volume of 24:1 chloroform:isoamylalcohol and homogenized by shaking the tube. The tube was centrifuged again at 15,000×g for 10min and the supernatant was mixed with 0.7 volumes isopropanol in a new tube and inverted gently. After another centrifugation step at 15,000×g for 10min, the supernatant was discarded and the pellet was washed with 1mL of 70% EtOH. After a final centrifugation at 15,000×g for 2min, the supernatant was discarded and the DNA pellet was allowed to dry at room temperature. The pellet was resuspended in 50μL of 10mM Tris-HCl, pH 8 with 0.1mM EDTA and stored at 4°C.

DNA was analyzed by field inversion gel electrophoresis to confirm that most molecules were around 48kb (Supplementary Fig. ^1f) before preparing the linked-read library (load: 1.00ng) using the Chromium Genome Reagent Kit, the Chromium Genome Library Kit & Gel Bead Kit, Chromium Genome Chip Kit v2, and Chromium i7 Multiplex Kit according to the manufacturer’s instructions. As 1.00ng is equivalent to 10⁶Mb, the haploid genome size is ~135 Mb⁴³, and each pollen grain is composed of three cells, there were ~3000 pollen grains expected with unique recombinant genomes.

In total, eight 10X and 50 standard Illumina libraries were sequenced on HiSeq 3000/4000 with 151bp paired-end reads, where the individual F₂s were at ~5x and the pools of F₂s and pollen nuclei were at 173–659×(Supplementary Data ¹).

Molecule recovery and genotyping using DrLink

Col-0 and Ler reference genomes were indexed using the function mkref of longranger (v2.2.2 10X Genomics). Linked-reads of each sample were aligned against the Col-0 and Ler reference genomes separately using longranger wgs with options:–id=ALIGN-ID–reference=MKREF-INDX-FOLDER–fastqs=READ-PATH–sample=READ-ID–localcores=40–localmem=192–noloupe–sex=male–vcmode=freebayes–library=LIBRARY-ID and all other options under default settings (note that option–sex was set as male, which is required by the tool but not affecting the alignment here). Read alignments with the same read barcodes were clustered to recover molecules by DrLink molecule function according to their genomic location ensuring that neighboring read alignments were not more than 25kb away from each other. The resultant molecules (over 1kb) were further filtered for falsely merged molecules by removing very long or densely covered molecules (Supplementary Note ¹). Next, with a set of SNP markers and the barcoded read alignments (in VCF) given by longranger, molecules were genotyped by DrLink recombis function to identify recombinant molecules, with options for filtering less confident predictions (Supplementary Data ³).

Generating the CO benchmark set

For each individually sequenced F₂, read alignments against the Arabidopsis Col-0 reference sequence³⁸ and variant calling were performed using Bowtie2⁵⁰ (version 2.2.8) and SAMtools/BCFtools⁵¹ (version 1.3.1). TIGER³³ was used to reconstruct the parental haplotypes using SNPs in co-linear regions between the Col-0 and Ler³⁹ genomes. We determined an average of 8.3 COs per diploid genome, totaling 415 COs with a median breakpoint resolution of 0.5kb (Supplementary Data ²). Of those, 15 COs were located in regions with an inter-marker distance of more than 10kb. As local CO breakpoint identification relies on markers near the breakpoint, we excluded these COs, generating a final benchmark set of 400 COs.

Estimating chromosomal CO landscapes

Recombination landscapes were estimated using sliding windows along each chromosome (window size 1Mb, step size 50kb). Within a given window, the CO frequency f was calculated by f=n / t_i, where n is the number of COs in the window and t_i is the total CO number within the respective chromosome i across each of the focal libraries. This ensures that different CO landscapes can be compared to each other within the given chromosome.

Comparison of CO distribution in pollen and F₂ populations

Local CO occurrence in pollen (P8000) and F₂ populations (two libraries: P1250R1/2) were compared based on sliding windows (window size 1Mb, step size 50kb). Specifically, for each window, we counted both recombinant molecules R^x and non-recombinant molecules N^x, where x was either P1250R1/2 or P8000. Then for each window, two Fisher’s exact tests were performed between (R^P8000, N^P8000) and (R^P1250R1/2, N^P1250R1/2). Finally, the windows with both p-values<0.05 were identified as different between P8000 and P1250R1/2, and neighboring windows were connected into larger regions.

Relative CO frequency estimation

For frequency estimation, we first subsampled molecules from the pools of interest by intersecting their molecule characteristic distributions. Using the molecule length distributions as the basis (Supplementary Fig. ^6a), the molecules were randomly sampled down within 1kb large bins. For each bin, the number of randomly selected molecules per pool was equal to the lowest number of molecules across all pools within this specific bin. To increase the subsampling rate, 80% of the sampled molecules were further randomly selected within each bin. After this second subsampling, the read and base coverage distributions were also highly similar (Supplementary Fig. ⁶b, c). Final molecule filtering was applied following the same strategy as used for the complete sets of molecules, by applying 30kb as size and 24 as read number thresholds.

Then, COs per million molecules (C^M) was calculated for each pool. The process of molecule subsampling, CO identification and C^M calculation was repeated fifty times for each pool. The means (μ) of all 50 C^M values and their 95% confidence intervals (μ±s, where s is 2.0096×σ/50^0.5 with σ being the standard deviation) were used to assess significant differences between samples.

DNA methylation level estimation

Within a recent study, we have assessed DNA methylation for A. thaliana Col-0⁵². Methylation level M for a CO or a random interval was calculated by M=N_met / N, where N_met is the number of reads supporting methylated cytosines at all C and G sites, while N is the total number of reads at these sites.

The authors would like to thank Ian R. Henderson (Department of Plant Sciences, University of Cambridge) for providing the CO breakpoint lists, Erik Wijnker (Wageningen University) for providing seeds and Ulrike Hümann, Manish Goel, Wen-Biao Jiao, Vidya Oruganti, and Onur Dogan (Max Planck Institute for Plant Breeding Research) for their help in the greenhouse. We also would like to thank 10X Genomics for their help on setting up the longranger software, advice on DNA extraction, and their kind donation of library reagents to support the development of recombination identification. We thank Lutz Froenicke and the DNA Technologies Core at UC Davis for 10X library sequencing support and discussions. We also acknowledge helpful discussions with Detlef Weigel at the Max Planck Institute for Developmental Biology and Kyle Fletcher, William Palmer, and Sebastian Reyes-Chin-Wo at UC Davis. We thank Felicity Jones and Frank Chan (Friedrich Miescher Laboratory of the Max Planck Society) for inspiring the extension of this work to gametes. This work was supported by the Max Planck Society postdoctoral fellowship, and a combined grant by the Deutsche Forschungsgemeinschaft (DFG) and the Agence Nationale de la Recherche (ANR) under grant number SCHN1257/8–1 (KS), and a UC Davis Genome Center Pilot Project grant (BAR).