
==== Front
Mol Biol Evol
Mol Biol Evol
molbev
Molecular Biology and Evolution
0737-4038
1537-1719
Oxford University Press UK

10.1093/molbev/msae183
msae183
Discoveries
AcademicSubjects/SCI01130
AcademicSubjects/SCI01180
Diversity in Recombination Hotspot Characteristics and Gene Structure Shape Fine-Scale Recombination Patterns in Plant Genomes
https://orcid.org/0000-0001-5990-7545
Brazier Thomas Unité Mixte de Recherche (UMR) 6553 - ECOBIO (Ecosystems, Biodiversity, Evolution), University of Rennes, CNRS, Rennes, France

https://orcid.org/0000-0001-7260-4573
Glémin Sylvain Unité Mixte de Recherche (UMR) 6553 - ECOBIO (Ecosystems, Biodiversity, Evolution), University of Rennes, CNRS, Rennes, France
Department of Ecology and Genetics, Evolutionary Biology Center and Science for Life Laboratory, Uppsala University, Uppsala, Sweden

Lee Grace Yuh Chwen Associate Editor
Corresponding author: E-mail: thomas.brazier@univ-rennes.fr.
Conflict of Interest None declared.

9 2024
20 9 2024
20 9 2024
41 9 msae18327 6 2024
20 8 2024
20 9 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of Society for Molecular Biology and Evolution.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact reprints@oup.com for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site—for further information please contact journals.permissions@oup.com.

Abstract

During the meiosis of many eukaryote species, crossovers tend to occur within narrow regions called recombination hotspots. In plants, it is generally thought that gene regulatory sequences, especially promoters and 5′ to 3′ untranslated regions, are enriched in hotspots, but this has been characterized in a handful of species only. We also lack a clear description of fine-scale variation in recombination rates within genic regions and little is known about hotspot position and intensity in plants. To address this question, we constructed fine-scale recombination maps from genetic polymorphism data and inferred recombination hotspots in 11 plant species. We detected gradients of recombination in genic regions in most species, yet gradients varied in intensity and shape depending on specific hotspot locations and gene structure. To further characterize recombination gradients, we decomposed them according to gene structure by rank and number of exons. We generalized the previously observed pattern that recombination hotspots are organized around the boundaries of coding sequences, especially 5′ promoters. However, our results also provided new insight into the relative importance of the 3′ end of genes in some species and the possible location of hotspots away from genic regions in some species. Variation among species seemed driven more by hotspot location among and within genes than by differences in size or intensity among species. Our results shed light on the variation in recombination rates at a very fine scale, revealing the diversity and complexity of genic recombination gradients emerging from the interaction between hotspot location and gene structure.

plant genomics
recombination hotspot
gene structure
recombination
crossovers
Agence Nationale de la Recherche 10.13039/501100001665 ANR-19-CE12472 0019
==== Body
pmcIntroduction

Meiotic recombination is a general feature of sexually reproducing species. During meiosis, crossovers (COs) ensure proper segregation of homologous chromosomes and reshuffle alleles from the two parental genomes. The distribution of COs is not homogeneous along chromosomes and varies at large and short scales (Haenel et al. 2018; Tock and Henderson 2018; Brazier and Glémin 2022; Lloyd 2022). Characterizing and understanding where recombination occurs have many implications because recombination has both local molecular effects (e.g. mutagenic effects, gene conversion) and indirect effects through the breakdown of linkage disequilibrium (LD) (Webster and Hurst 2012).

In many species described so far, recombination tends to occur within narrow regions called recombination hotspots (Kauppi et al. 2004; Mézard 2006; de Massy 2013; Lam and Keeney 2015), but in some species, such as Drosophila and Caenorhabditis species, proper hotspots are lacking. When they exist, the general view is that two main mechanisms determine where hotspots occur across the genome. In many mammals, the location of recombination hotspots is directed by the PRDM9 protein (Baudat et al. 2010; Grey et al. 2018). The zinc-finger domain that binds to specific DNA motifs is evolving fast, thus changing recombination landscapes over short evolutionary times (Boulton et al. 1997; Latrille et al. 2017). The PRDM9 gene is believed to be conserved across several vertebrate lineages because it redirects recombination away from the genic region to avoid deleterious mutations and chromosomal rearrangements (Baudat et al. 2010; Smagulova et al. 2011; Baker et al. 2017). Recently, PRDM9-driven hotspots have also been identified in nonmammalian species (Raynaud et al. 2024). In contrast, many animals but also plants and fungi, do not have a PRDM9-like system (Baker et al. 2017). In such species where PRDM9 is lacking, recombination was found to occur preferentially in gene regulatory sequences, especially in promoters and 5′ to 3′ untranslated regions (UTRs) (Myers 2005; Brick et al. 2012; Auton et al. 2013a; Choi et al. 2013; Choi and Henderson 2015; Li et al. 2015). Because this was found in very phylogenetically divergent species such as plants, yeasts and birds (Pan et al. 2011; Choi et al. 2013; Singhal et al. 2015) and when PRDM9 is experimentally knocked-out (Brick et al. 2012), it is thought to be the ancestral mechanism for locating recombination hotspots. Mechanistically, recombination hotspots seem globally associated with nucleosome occupancy and methylation patterns that define flanking regions of genes (Choi et al. 2013; Shilo et al. 2015; Tock and Henderson 2018). Although recombination hotspots can overlap coding sequences (CDSs), both mechanisms tend to lower recombination rates within genes compared to flanking regions and/or nongenic regions.

When hotspots are targeted to promoter like features it can generate recombination gradients within genes, as observed in a few plant species (Dooner and Martínez-Férez 1997; Fridman et al. 2000; Hellsten et al. 2013; He et al. 2017; Marand et al. 2019a). In rice, recombination hotspots overlap genes between the transcription starting site (TSS) and the transcription termination site (TTS), potentially increasing recombination in the first and last exons (Marand et al. 2019a). In Mimulus guttatus, recombination is higher in the first exon than in 5′ noncoding sequences (Hellsten et al. 2013). Even if recombination hotspots strictly remain in flanking regions, recombination tracts span hundreds of base pairs and are likely to end within genes. Indeed, 5′ ends of genes experience gene conversion gradients (Schultes and Szostak 1990; Detloff et al. 1992; Malone et al. 1992) and the combination of GC-biased gene conversion (gBGC) and 5′ to 3′ recombination gradients would provide an explanation for gradients of GC nucleotide content along genes commonly observed among plant species (Serres-Giardi et al. 2012; Glémin et al. 2014; Ressayre et al. 2015; Clément et al. 2017). Therefore, actual GC gradients are indirect evidence that recombination gradients could be widespread and organized as a function of the number of exons in addition to the distance to TSS and TTS, as observed for GC content (Glémin et al. 2014; Ressayre et al. 2015).

However, this simple dichotomous view has been recently challenged by the exploration of recombination landscapes in various non model species. In mammals, human and mouse, where PRDM9 mechanism was initially identified and characterized in detail, appear as extremes in a more continuous distribution where PRDM9 hotspot could co-occur with PRDM9-independent hotspots directed to promoter regions (Joseph et al. 2024). Such co-occurrence of PRDM9 and promoter hotspots has also be observed in a snake species (Hoge et al. 2024). Similarly, we do not know whether the precise location of hotspots around promoter features is the same and conserve across kingdoms. For example, in some animal species such as dog and snake, hotspots are preferentially located around TSS and CpG islands (Auton et al. 2013a; Hoge et al. 2024), whereas in some birds and plants they are associated with both TSS and TTS (Choi et al. 2013; Singhal et al. 2015; Marand et al. 2019a). In plants, the fine-scale recombination landscape has been characterized in a handful of species only, and the implicit view that plants have recombination hotspots mainly targeted to TSS and partly to TTS as in Arabidopsis thaliana deserves further investigation. In addition, the simple observation of recombination gradient within genes can be misleading if the gene structure is not properly taken into account, as exemplified with gradients of GC content (Ressayre et al. 2015).

Here, we extended the characterization of fine-scale recombination landscapes in diverse flowering plants, specially focusing on genic regions. We assessed (i) whether recombination hotspots are common in plants or whether some species lack them, such as Drosophila and Caenorhabidtis; (ii) whether gene promoters and terminators are enriched in COs or whether hotspots can also be targeted outside genic regions in some species; and (iii) how hotspot characteristics and gene structure (exons/introns) interact to shape the recombination gradient within genes. To do so, we constructed fine-scale recombination landscapes in 11 plant species, including nine new species (plus A. thaliana and rice), and we reanalyzed human data to compare with a species with PRDM9. We inferred fine-scale recombination rates from patterns of LD in polymorphism data (Auton and McVean 2007). Owing to tiny distances between genetic markers (1 to 2 kb) and the large number of meiosis observed in the genealogy of a population, variation in population-scaled recombination rates can be measured within genes and can be used to infer the precise hotspot locations (Myers 2005; Choi et al. 2013; Auton et al. 2014; Stukenbrock and Dutheil 2018).

We found signature of hotspots in all 11 species but we uncovered a diversity of hotspot characteristics with preferential location in TSS (as assumed to be common) but also more frequent in TTS for some species, or without clear location around genic features. Despite this diversity, we propose that a simple model that takes gene structure (size, number of exons) and hotspot position into account can explain the variety of recombination gradients we observed within genes.

Results

LD-Based Recombination Landscapes

To achieve fine-scale LD-based recombination maps at a gene scale, we gathered high-density polymorphism datasets in 11 flowering plant species (Table 1, supplementary table S1, Supplementary Material online). We identified species of interest, representing the diversity of plant genomes (small/large chromosomes, eudicots/monocots) and broad-scale recombination patterns (low/high genome-wide recombination rates), from a previous study (Brazier and Glémin 2022). We also used a human dataset (Sudmant et al. 2015) to compare the results with a species with hotspots targeted outside genic regions due to the PRDM9 mechanism. To ensure the comparability of the results we run the same pipeline on this human dataset instead of directly using published recombination maps. After filtering SNPs (minor allele frequency >0.05, missing data per site <0.1, biallelic SNPs, genotype quality score ≥30), we kept only datasets with at least 0.5 SNP/kb (range 0.6 to 9.7 SNP/kb, Table 1).

Table 1 Summary statistics of 12 datasets (11 plants and 1 human)

Species	Reference	# SNPs	Mb	Map masked (%)	SNP/kb	# Chr.	# Diploid genomes	Mean ρ/kb	Median ρ/kb	θ^π	ρ/θ	
A. thaliana	The 1001 Genomes Consortium (2016)	641,453	119	4.7	5.39	5	40	2.68	1.21	0.0014	1.97	
C. sinensis	Zhang et al. (2021)	1,808,673	2,985	42.3	0.6	15	40	4.64	4.60	0.0002	25.75	
C. lanatus	Guo et al. (2019)	366,856	362	11.7	1.0	11	40	1.53	0.07	0.0004	4.12	
G. max	Yang et al. (2021)	4,337,396	949	12.2	4.57	20	40	12.34	1.51	0.0012	10.38	
H. sapiens	Sudmant et al. (2015)	6,584,165	2,796	4.7	2.3	22	40	3.37	0.57	0.0008	3.99	
M. sieversii	Sun et al. (2020)	6,295,986	651	1	9.66	17	37	8.06	3.25	0.0036	2.21	
O. sativa indica subsp.	Wang et al. (2018)	1,710,624	372	8.7	4.59	12	40	2.77	2.85	0.0021	1.34	
P. vulgaris	Wu et al. (2020)	1,840,346	515	31.7	3.6	11	40	7.11	2.23	0.0011	6.69	
P. tremula	Liu et al. (2022)	1,663,613	361	11.2	4.6	19	36	4.67	2.70	0.0011	4.24	
S. bicolor	Lozano et al. (2021)	422,937	683	30.9	0.6	10	32	0.44	0.03	0.0002	2.77	
S. oleracea	Cai et al. (2021)	3,283,574	879	0.5	3.7	6	40	2.68	1.52	0.0011	2.53	
T. aestivum B subgenome	Zhou et al. (2020)	15,614,566	5,177	27.8	3.0	7	40	3.12	2.14	0.0004	7.69	
Species, reference of the original polymorphism data, number of SNPs after filtering, genome length (Mb) and part of the LD map masked (%), SNP density (number of SNPs per kb), number of chromosomes of the genome, number of diploid genomes sampled, mean and median ρ/kb, estimated θ^π per bp and ratio mean ρ/θ.

Fine-scale recombination landscapes have been estimated with LDhat (Auton and McVean 2007) using a custom pipeline (https://github.com/ThomasBrazier/ldhat-recombination-pipeline.git v1.1). We checked for hidden population structure and sampled individuals within a single consistent genetic group with the highest polymorphism level, accordingly. For selfing species we considered haploid genome by sampling randomly one allele at the few heterozygote SNPs. The population size parameter (θπ) was estimated from the sampled population. The historical demography of the population was estimated with SMC++ (Terhorst et al. 2017) and used to generate a demography-aware look-up table with LDpop (Kamm et al. 2016). For computational limits, a maximum of forty diploid genomes per species (80 haploid genomes for selfing species, autosomes only) were sampled for LDhat analyses since the quality of estimates does not improve much over twenty individuals (Raynaud et al. 2023). We checked the reliability of LD-based recombination landscapes by visual comparison (supplementary fig. S1, Supplementary Material online) between LD-based landscapes and broad-scale pedigree-based landscapes (Marey maps) of Brazier and Glémin (Brazier and Glémin 2022). We filtered all LD-based recombination maps to remove large genomic segments (>100 kb) with a constant recombination rate (see the percentage masked in Table 1).

The mean genome-wide population-scaled recombination rates ranged from 0.44 to 12.34 ρ/kb across species (median range 0.03 to 4.6 ρ/kb). Most of the recombination rate estimates were lower than 10, yet the tails of the distribution were consistently skewed towards high values across species (supplementary fig. S2, Supplementary Material online). The ratio ρ/θ were globally in the same order of magnitude as in human (Table 1), except for Camellia sinensis (ratio ∼25.7). LD-based methods have been mostly evaluated on a range of parameters close to human, with a ρ/θ between 0.1 and 10 not affecting much the accuracy of estimates except for low mutation rates (μ=10−9) (Raynaud et al. 2023).

Recombination Hotspots Are Common but with Different Characteristics among Species

Based on the LDhot statistical framework we inferred LD-based recombination hotspots in 11 plant species as well as in human in a standardized and comparable manner (Table 2, Fig. 1a). We applied two levels of filtering to the hotspot call set. As recombination hotspots are generally defined as narrow peaks of recombination over a short genomic range, we first filtered out all hotspots larger than 10 kb without making assumptions about their intensity (soft filtering). In order to improve our power to detect fine-scale patterns and associations with short targets (e.g. TSS, TTS), we did a harder filtering removing every hotspot larger than 10 kb and with an intensity lower than 4 and higher than 200 (putatively false hotspots due to variance in LDhat estimates). The hotspot intensity was measured as the peak rate divided by the background recombination rate (the mean background rate in a 50 kb window around the hotspot center). The two levels of filtering drastically reduced the number of hotspots (Table 2). Hotspot filtering had the same impact on hotspot shape for all species. Hotspot size was reduced to approximately 6 to 8 kb and the hotspot intensity increased, as expected by the two criteria we chose (Fig. 1a).

Table 2 Number of hotspots (density of hotspots per Mb of the masked genome length)

Species	# Hotspots	Density/Mb	# (Soft)	% (Soft)	Density/Mb (soft)	# (Hard)	% (Hard)	Density/Mb (hard)	# (Literature)	Reference	
A. thaliana	5,249	46.3	2,685	51.1	23.7	889	16.9	7.83	8,448	Choi et al. (2013)	
C. sinensis	41,843	24.3	14,098	33.7	8.18	5,299	12.7	3.08	…	…	
C. lanatus	9,068	28.3	4,511	49.7	14.1	2,726	19.0	5.4	…	…	
G. max	69,499	83.4	66,287	95.4	79.5	23,467	33.8	28.2	…	…	
H. sapiens	104,573	39.3	57,952	55.4	21.7	16,952	16.2	6.4	25,000 to 50,000	Myers (2005)	
M. sieversii	26,286	43.8	13,814	48.8	21.4	4,054	14.3	6.3	…	…	
O. sativa	12,443	36.6	6,134	49.3	18.0	3,251	26.1	9.5	14,125	Marand et al. (2019a)	
P. vulgaris	13,550	38.5	7,573	55.9	21.5	2,886	21.3	8.2	…	…	
P. tremula	13,804	43.0	6,799	49.2	21.2	1,620	11.7	5.1	6,000	Slavov et al. (2012) in P. trichocarpa	
S. bicolor	12,350	26.1	6,324	51.2	13.4	2,859	23.1	6.0	…	…	
S. oleracea	75,376	86.2	73,431	97.4	84.0	30,413	40.4	34.8	…	…	
T. aestivum	140,067	37.5	71,892	51.3	19.2	28,180	20.1	7.5	…	…	
Statistics are given for the unfiltered, soft and hard filtering datasets. The number of hotspots already known in the literature is given at the end.

Fig. 1. The distribution of recombination is heterogeneous at a fine scale and concentrated within recombination hotspots. a) The recombination rate (ρ/kb) as a function of the genomic distance to the CO hotspot center (kb) with three different filtering strategies. Soft filtering removed hotspots larger than 10 kb. Hard filtering removed hotspots larger than 10 kb and hotspots with an intensity lower than 4 or higher than 200. Fully colored lines are the recombination rates around hotspot center (midpoint of the hotspot interval) and shaded lines are the control by randomly resampling hotspot intervals within a ±50 kb neighboring region. The dashed horizontal line is the mean genome-wide recombination rate. b) The cumulative ordered recombination fraction as a function of the cumulative genome fraction for each chromosomes and species. We represented dashed vertical lines at 20% of the genome and dashed horizontal lines at 80% of the total recombination. c) Normalized (centered-reduced) recombination rates as a function of the distance to the hotspot center for the soft-filtered dataset. Species in the legend ordered by descending peak height. Human is in black.

The genome-wide hotspot density varied by a factor of 3.5 among species, independently of the filtering strategy (Table 2) and recombination rates varied over many orders of magnitude at a fine scale and were concentrated in a short fraction of the genome (Fig. 1b, supplementary fig. S2, Supplementary Material online). For most species, 80% of the total recombination was concentrated within approximately 20% of the genome. The species-specific patterns were roughly the same between plants and human. Our results for human, rice and A. thaliana were globally similar to state-of-the-art previous studies (Myers 2005; Choi et al. 2013; Auton et al. 2014; Marand et al. 2019a), both for the shape of hotspots, the cumulative distribution function (Fig. 1a,b) and the number of hotspots (Table 2). Overall our results showed that recombination was globally concentrated in a short fraction of the genome across 11 plant species, suggesting that recombination hotspots may be a general feature in plants, similar to human.

Despite qualitatively similar genomic sizes, hotspots varied in intensity among species (Fig. 1a,c). Three species exhibited particularly low levels of hotspot intensity compared to the background rate (Citrullus lanatus, Malus sieversii and Phaseolus vulgaris) while Spinacia oleracea and Populus tremula were the most intense plant species. Human presented an average pattern. Comparing the background recombination rate (±50 kb) and the genome-wide rate, showed that hotspot detection was sometimes biased in favor of genomic regions recombining on average slightly less or more than the genome-wide average (see the difference between soft lines and the dashed line in Fig. 1a). For the unfiltered and soft-filtered datasets, the bias was absent in three species and weak in eight other species (except C. lanatus), suggesting that we were not strongly limited to detect hotspots in regions experiencing extremely high or low recombination rates. However, for the hard-filtered dataset, the detection of hotspots was particularly biased towards lowly recombining regions in many species (see the difference between the hard filtering and the dashed line), probably because intense hotspots are more easily detected on a low recombining background. In general, the soft-filtered dataset seemed relatively more reliable for the exhaustive search of hotspot location along the genome since hard filtering was biased towards low recombining regions in many species. Subsequently, we present results for the soft-filtered dataset.

The Fine-Scale Distribution of Recombination Varies among Species

As already observed in human, recombination rates were lower in genic than intergenic regions, except in P. tremula and C. sinensis, even when flanking regions containing promoters and other cis-regulatory elements are excluded (Fig. 2a). Here, we used the median instead of the mean to limit the effect of extreme values but using the mean recombination rates provided rather similar results with larger confidence intervals (supplementary fig. S3a, Supplementary Material online). Intergenic regions could be more prone to mapping errors, which could artificially inflate recombination rates. If it is a true signal, it suggests that recombination could also be directed towards intergenic regions, at least in some species. Despite lower recombination rate on average, genic regions were, however, enriched in recombination hotspots in 8 over the 11 species: the number of hotspots overlapping a gene feature was significantly higher than expected (random expectation computed by 1,000 iterations of random shuffling of hotspot ranges) (Table 3, supplementary table S2, Supplementary Material online). Notably, it was also the case in human. In contrast, in Glycine max, S. oleracea, Triticum aestivum the number of genic hotspots was less than expected just by chance.

Fig. 2. Variations in median recombination rates (ρ/kb) around and within genes. a) Recombination rates vary between genic and intergenic regions. In order to remove an effect of 5′ and 3′ regulatory regions at the proximity of genes, buffered intergenic regions were defined by excluding 3 kb flanking regions upstream and downstream of genes. b) Within genic regions, the recombination rate varies between genomic features as well as in 5′ upstream and 3′ downstream flanking regions. The median recombination rate and 95% CIs were estimated by 1,000 bootstraps. Annotations of UTRs were not available for C. lanatus.

Table 3 Number of genic and intergenic hotspots, i.e. hotspots overlapping or not a gene feature (soft-filtered data)

Species	# Intergenic hotspots	# Genic hotspots	# Expected genic (95% CI)	
A. thaliana	215	2,470	655 (614 to 700)	
C. sinensis	8,449	5,649	813 (760 to 867)	
C. lanatus	2,629	1,882	687 (641 to 729)	
G. max	50,076	16,211	21,144 (20,917 to 21,359)	
H. sapiens	43,122	14,830	12,894 (12,703 to 13,091)	
M. sieversii	7,414	6,400	3,101 (3,011 to 3,194)	
O. sativa	2,263	3,871	1,147 (1,090 to 1,206)	
P. vulgaris	4,308	3,265	1,348 (1,287 to 1,413)	
P. tremula	1,356	5,443	1,463 (1,393 to 1,528)	
S. bicolor	3,238	3,086	759 (712 to 807)	
S. oleracea	58,161	15,270	23,883 (23,645 to 24,121)	
T. aestivum	68,074	3,818	11,457 (11,267 to 11,631)	
The expected number of genic hotspots was computed by 1,000 random shuffling of hotspots ranges.

Focusing on genic regions, recombination within UTRs and flanking regions was slightly higher than in coding regions (Fig. 2b, supplementary fig. S3b, Supplementary Material online). In most species, such as G. max and Oryza sativa, recombination rates were lower in UTRs than in flanking regions. However, in P. tremula and C. sinensis recombination was only higher in the 3′ UTR and flanking region, with recombination being similar or higher in exons, introns than in 5′ part. We did not detect any clear pattern for differences in averaged recombination rates between exons and introns.

Based on Fig. 2, we chose three species (A. thaliana, P. tremula, G. max) illustrating the diversity of recombination patterns to characterize more precisely 5′ to 3′ genic recombination gradients (Fig. 3). We found patterns similar to one of these three examples in other species (supplementary fig. S4, Supplementary Material online). In A. thaliana, the 5′ flanking region and UTR exhibited a rather large peak of recombination spanning about 1 kb and recombination rate then steeply decreased from 5′ to 3′ just after the start codon. The 3′ end of the gene showed a small peak within the CDS just before the TTS but was much weaker than the 5′ end plateau. The 5′ peak was associated with a corresponding peak of hotspots density around TSS but we did not detect hotspot enrichment around TTS. This pattern was in agreement with previous results and considered as the standard pattern in plants (Choi et al. 2013). Populus tremula showed marked differences with the pattern observed in A. thaliana. We observed two narrow peaks inside the CDS, centered on the CDS start instead of the TTS and on the 3′ UTR instead of the TTS. Contrary to A. thaliana and G. max we also observed a higher peak in 3′ than in 5′, which both corresponded to an hotspot enrichment. Finally, in G. max we did not observed any peak but only a weak decrease towards the inside of the gene from both 5′ and 3′ ends, without hotpost enrichment (note also the different scales for the recombination rate).

Overall, these results clearly showed that the A. thaliana pattern is not universal in flowering plants and that some species, such as G. max, P. vulgaris, C. lanatus, and possibly T. aestivum may not share the supposed ancestral system of targeting recombination towards promoter features or more generally gene flanking regions.

Fig. 3. Recombination gradients along genes are mostly around TSS and TTS. a,b) The median recombination rate (ρ/kb) was estimated in 200 bp windows as a function of the distance to the ATG start codon (ATG codon, the start of first CDS) or the TTS and averaged among introns and exons. The vertical dashed lines are, respectively, the mean distance to the TSS (i.e. the 5′ UTR, panel a) and the mean 3′ UTR length (panel b). Random controls by randomly resampling intervals in the genome and rescaled to the mean of the gradient. c,d) The density of TSS and TTS as a function of the distance to hotspot center (soft filtered hotspots). For smoothing, TSS and TTS were defined as regions of 1 kb centered on the gene start and end position, respectively. Soft lines are a control by randomly resampling hotspot intervals in a ± 50 kb neighboring region.

Hotspot Location and Gene Structure Shaped Gradients of Recombination within Genes

Previous results suggested that whatever the hotspot location, recombination gradients occurred inside genes. To go further we decomposed genic gradients per exons/intron ranks to remove spatial averaging effects among genes of different lengths and to evaluate the distribution of recombination within gene bodies. We kept only genes with less than 15 exons, as the sample size was too low above 14 exons (supplementary table S3, Supplementary Material online).

We first considered the average gradient by pooling all genes (see the black line in Fig. 4a). In A. thaliana and G. max, recombination decreased over the five to seven first exons and reached a minimal plateau with a slight reincrease at the end. In contrast, in P. tremula, after a decrease in 5′, recombination rates strongly reincreased from the middle of the gene to the 3′ end, yielding a U-shape gradient on average. We observed only slight differences in gradients among exons and introns (Fig. 4, supplementary fig. S6, Supplementary Material online).

Fig. 4. Gradients of recombination (median ρ/kb) along genes as a function of exon/intron rank. a) Comparison of independent gradients as a function of exon rank (CDS part) with genes pooled by their number of exons. The black line is the average gradient (all genes pooled). b) Average gradient of recombination rate as a function of exon of rank i (ei points) or intron of rank i (ii points). c) Gradients in flanking regions representing the differences of recombination rates between genic regions and their flanking noncoding regions. On each side of the gene, the recombination rate of the 5′ and 3′ flanking regions (3 kb windows) is represented by a black dot and the difference between the gene part and its close flanking region is represented by a black line.

We then evaluated how gene structure (rank and number of exons) could shape differences in gradients among species. As such, gradients can also be observed by pooling genes according to their number of exons. We decomposed the averaged gradients across ranks into independent gradients for each exon-number class of genes (colored lines in Fig. 4a). In all three species, shorter genes (genes with fewer exons) recombined more than longer genes on average (Fig. 4a, supplementary fig. S7, Supplementary Material online). However, it was not entirely driven by lower recombination rates in the middle of genes. The first exon(s) of shorter genes experienced higher recombination rates than the first exons of longer genes. Gradients per gene class were more or less parallel to the average gradient in A. thaliana, though a slight increase was detected in 3′ ends of longer genes (longer than seven exons). Gradients per gene class were weaker in G. max and the average gradient was mostly shaped by an excess of recombination in shorter genes (shorter than five exons). The structure of P. tremula gradients were more complex than the two previous patterns. The U-shape average gradient was produced by a combination of per gene class gradients increasing in 3′ and higher recombination rates in the first exons of shorter genes. In P. tremula gradients per gene class were strongly polarized towards the 3′ end. Importantly, the average gradient also depends of the proportion of each class of genes. As short genes are more frequent than long ones, this biases the average gradient towards the 5′ end.

We checked if recombination rates at gene boundaries were correlated with the local recombination rate within flanking regions (Fig. 4c). In A. thaliana we noted a 5′ flanking recombination rate globally similar to the recombination rate of the first exon but elevated recombination rates in 3′ flanking regions compared to the last exon preceding them. In G. max, the differences with flanking regions was even stronger. The gradient within the CDS was weak compared to the increase in recombination rates in flanking regions. In P. tremula, recombination rates at 5′ and 3′ flanking regions were similar among gene sizes and decoupled from variations within the transcript sequence.

Recombination Gradient Inference Is Robust

LD-based recombination rate estimates are population-scaled and ρ is the product of the effective population size Ne and the CO rate r (ρ=4Ner). As such ρ estimates are potentially biased if underlying variations in Ne are not properly taken into account (Dapper and Payseur 2018; Samuk and Noor 2022). Though we properly controlled for the genome-averaged effect of population structure and demography during the estimation of ρ (see material and methods for details), fine-scale variations in Ne due to selection could also potentially impact LD-based estimates (O’Reilly et al. 2008; Barroso and Dutheil 2023). To assess the robustness of our results to these potential biases we compared recombination gradients to fine-scale patterns of polymorphism (SNP density and genetic diversity θπ) pooled in the same manner (Fig. 5a,b). We observed a weak noisy gradient of SNP density in A. thaliana but a clear absence of gradient in G. max and P. tremula (Fig. 5a). The recombination gradients we observed are not likely byproducts of fine-scale variations in polymorphism. The ratio ρ/θπ followed a very similar gradients to ρ/kb. These results were strong evidence that recombination gradients were true gradients and not artifacts generated by underlying variations in selection (Fig. 5b).

Fig. 5. Recombination gradients are robust to variations in polymorphism levels and correlate with CO rates estimated in Rowan et al. (2019). a) Gradients of SNP density (SNP/kb) along genes as a function of CDS rank. SNP density was estimated for each CDS part and the average SNP density was calculated for each CDS part rank and gene size (number of exons). b) Gradients of ρ/θπ along genes as a function of CDS rank. θπ was estimated per site and the average θπ/bp was calculated by the sum of θπ divided by the total sequence length of exons. ρ/θπ is ρ/kb divided by θπ/kb. c) Peaks of LD-based recombination rates (ρ/kb) around CO center (± 5 kb) estimated in Rowan et al. (2019). d) LD-based hotspots are enriched in COs. Mean CO count around the hotspot center (± 20 kb). Mean CO count measured by counting the number of COs overlapping 200 bp windows around the LD hotspot center. e,f) TSS and TTS are enriched in COs (± 20 kb). The blue line is loess regression with 95% CI

To confirm that LD-based recombination rates were effectively correlated to estimates of CO rates, we compared our results to the most resolved pedigree-based genetic map in the plant model A. thaliana (Rowan et al. 2019). The two estimates were highly congruent despite different approaches and data. The LD-based ρ/kb mapped on 17,077 COs (pedigree-based) showed a peak of ρ/kb around the CO center (Fig. 5c). Conversely, the mean CO count around the LD-based hotspot center showed that LD-based hotspots were enriched in COs (Fig. 5d). Despite a relatively low number of COs compared to the number of genes (17,077 COs vs. 25,942 genes), we were able to detect an increase in CO count around the TSS and TTS (Fig. 5e,f).

The Diversity of Gradients across Species Can Be Explained by a Single Model

Applying the same analyses to the other plant species and to human we observed a diversity of patterns more or less similar to the three species we detailed above. Recombination patterns were mostly determined by differences in hotspot location and intensity. Average 5′ to 3′ recombination gradients were common across the 11 plant species and even in human (Fig. 6, black lines). The less defined gradient was in T. aestivum for which the power to detect gradients was limited (see methods for details). The level of LD was high in this species and the LD decay was very slow, even in genic regions, suggesting a lack of power to detect fine-scale patterns (supplementary fig. S8, Supplementary Material online). This diversity of recombination gradients was parallel to hotspot gradients (supplementary figs. S5 and S9, Supplementary Material online). As a consequence, it seemed general among plants that shorter genes recombine more on average (Fig. 6, supplementary fig. S7, Supplementary Material online).

Fig. 6. Diversity of recombination gradients across plant species. The recombination rate (median ρ/kb) was estimated in exons (CDS part) as a function of their rank and genes were grouped by their number of exons. The black line is the average gradient (all genes pooled). The three reference species already shown above are transparent. Only protein-coding genes were used.

To better understand this diversity of patterns, we formalized a simple model based on the Choi and Henderson conceptual model (2015). We assumed that recombination exponentially decreased within gene from both ends, starting with recombination rates in 5′ and 3′ flanking regions, r5 and r3, depending on hotspots location and intensity (see equations (1) to (3) in Material and Methods). Under this mechanistic view, exon and intron lengths, which vary with gene structure (supplementary fig. S10, Supplementary Material online) should play a role and took them explicitly into account. Because we focused on within gene gradients, this model corresponds either to hotspots located in 5′ and 3′ regions with different proportion and/or intensity, such as in A. thaliana and P. tremula, but also if hotspots are targeted away from genic regions, with no recombination peak in 5′ and 3′ but a simple decrease in recombination rate within genes, as in G. max. We fitted this model to the detailed gradients shown in Fig. 6. This simple model fitted the data very well (r2 from 0.64 to 0.96, Table 4, Fig. 7b, c, supplementary fig. S12, Supplementary Material online) and the ratio r5/r3 was congruent with the observed gradients.

Fig. 7. Fit of the model to recombination gradients. a) Conceptual model of recombination gradients adapted from Choi and Henderson (2015). Recombination rates exponentially decay around the TSS and TTS. This decay can be due to hotspots in TSS and/or TTS, or loss of recombination within the coding sequence. Our model only fits recombination rates in exons. b) Recombination rates predicted by the model, as described in equations (1) to (3). Model fitted with nonlinear least squares and the sum of square differences between expected and observed recombination rates as a loss function to minimize. c) Goodness of fit of the model, with observed values as a function of predicted values.

Table 4 Fit of the model and parameter estimates for 11 species

Species	r0	r0 (SE)	r5	r5 (SE)	r3	r3 (SE)	Ratio r5/r3	a	a (SE)	Tract length (bp)	R2	
A. thaliana	− 0.74	0.09	1.43	0.05	1.17	0.05	1.23	0.000505	0.00003	1,982	0.95	
C. sinensis	0.49	0.24	0.69	0.11	1.59	0.16	0.44	0.000128	0.000023	7,817	0.76	
C. lanatus	− 3.58	1.29	2.35	0.7	1.62	0.59	1.45	0.000042	0.000011	23,950	0.72	
G. max	− 0.01	0.04	0.35	0.02	0.28	0.02	1.26	0.000302	0.000031	3,315	0.84	
H. sapiens	0.36	0.02	0.23	0.02	0.15	0.02	1.56	0.000067	0.000012	14,962	0.71	
M. sieversii	0.77	0.07	0.66	0.06	0.99	0.06	0.67	0.000629	0.000076	1,590	0.78	
O. sativa	− 0.07	0.06	0.56	0.04	0.48	0.03	1.15	0.000294	0.000028	3,399	0.86	
P. vulgaris	− 0.08	0.09	0.48	0.05	0.46	0.05	1.04	0.000175	0.000026	5,712	0.84	
P. tremula	− 1.77	0.23	3.4	0.14	5.38	0.15	0.63	0.000442	0.000025	2,263	0.96	
S. bicolor	− 0.02	0.02	0.09	0.01	0.06	0.01	1.44	0.00028	0.000047	3,576	0.69	
S. oleracea	− 0.05	0.03	0.65	0.02	0.59	0.02	1.1	0.000239	0.000013	4,180	0.95	
T. aestivum	− 0.04	0.13	0.44	0.07	0.43	0.07	1.01	0.000186	0.000041	5,365	0.64	
The basal recombination rate r0 (with SE), the 5′ recombination rate r5 (with SE), the 3′ recombination rate r3 (with SE), the ratio r5/r3, the inverse of the tract length a (with SE), the tract length in bp, and the R2 of the nonlinear least squares fitting procedure. The standard errors of the predictions were obtained with the “nls” R function.

Discussion

We characterized recombination patterns in genic regions in 11 plant species, including nine species not described before. We detected recombination hotspots in all species and we generalized the previously observed pattern that recombination hotspots are organized around the boundaries of CDSs, both at the 5′ and 3′ ends of genes. However, we also uncovered more variations than initially envisaged from only a few model species, including some species without clear targeting of recombination hotspots around genic regions. We also leveraged high-resolution recombination maps at the exon scale to characterize detailed genic recombination gradients and we proposed a single model that can explain the diversity of observed patterns through the interaction between hotspot location and gene structure.

Recombination Varies around and within Coding Sequences

At the gene scale, the location of recombination hotspots was often concentrated in gene 5′ and 3′ ends, as previously known for a few plants only (Choi et al. 2013; Hellsten et al. 2013; He et al. 2017; Marand et al. 2017, 2019a). We observed that recombination can be targeted towards diverse domains both in 5′ and 3′ ends (UTRs, TSS/TTS, promoters, flanking) and was not systematically associated with promoters as previously found (Hellsten et al. 2013; Wijnker et al. 2013; He et al. 2017). Despite divergent patterns, it is clear that the genic region is finely structured according to the different genomic domains composing it.

Within 5′ and 3′ ends enriched in COs in A. thaliana, chromatin accessibility redirects meiotic recombination towards gene promoters and terminators (Choi et al. 2013), and the same mechanism also likely occurs in other species (Lloyd 2022). In plants as in yeast, double-strand breaks (DSBs) hotspots, the precursors of COs, are found both in TSS and TTS where nucleosome occupancy is reduced, whereas COs form close to H2A.Z and H3K4me3 histone marks which are epigenetic modifications of chromatin packaging also involved in the regulation of gene expression (Jessop et al. 2005; Choi et al. 2013; Kianian et al. 2018). Low methylation is also usually associated with CO hotspots (Choi et al. 2013). Both TSS and TTS are hypomethylated in A. thaliana, G. max, and P. tremula (Vining et al. 2012; Song et al. 2013; Yelina et al. 2015), which fits well with our results. These findings are globally in agreement with the tethered-loop/axis model, which postulates that DSBs form globally in chromatin loops bound to the chromosome axis, whereas COs occur only in nucleosome-depleted regions that are susceptible to being tethered to the central chromosome axis by recombination-promoting factors (Blat et al. 2002; Tock and Henderson 2018). Indeed, gene promoters and terminators are located in chromatin loops and COs most probably concentrate here due to their highly conserved chromatin state and DNA accessibility (Cooper et al. 2016; Tock and Henderson 2018). However, P. tremula and C. sinensis do not conform exactly to this canonical model as CO hotspots preferentially occur inside rather than outside CDSs, as also seen in M. guttatus (Hellsten et al. 2013). In PRDM9−/− mice, the center of PRDM9-independent hotspots is on the +1 nucleosome positioned downstream of the TSS (Brick et al. 2012). For these species, the detailed description of chromatin and histone patterns during meiosis (especially H2A.Z and H3K4me3), as well as nucleosome maps, could be useful to better understand this noncanonical pattern (Tock and Henderson 2018; Lloyd 2022).

We also identified species, G. max, P. vulgaris, C. lanatus, and T. aestivum (supplementary figs. S4 and S5, Supplementary Material online), where we did not identify a clear signature of hotspot location around 5′ and/or 3′ regions, and some species with hotspots possibly targeting away from genic regions G. max, S. olearacea, and possibly T. aestivum, which showed a significant deficit in hotspot in genic regions (Table 2). Interestingly, P. vulgaris and C. lanatus showed an excess of hotspot in genic regions without clear 5′/3′ pattern, whereas S. olearacea showed a deficit of hotspots in genic regions but clear peaks both around TSS and TTS. If it is a true signal, it would suggest the existence of additional mechanisms, different from the supposed ancestral mechanism, which can direct recombination towards other targets than promoter-like sequences. In addition, also in contrast with previous expectations, we found a signature of 5′ hotspots with a weak 5′ to 3′ recombination gradient in human, despite the PRDM9 mechanism. This reinforces the recent findings in some animals that several competitive mechanisms could be at play in locating recombination hotspots within a genome (Hoge et al. 2024; Joseph et al. 2024).

Genic Recombination Gradients Are Shaped by Hotspot Location and Gene Structure

The specific location of hotspots can generate recombination gradient within genic regions, which can be used as a signature to distinguish promoter-driven and PRDM9-driven mechanisms (Singhal et al. 2015; Hoge et al. 2024). However, in previous studies, only average gradients were considered. We showed that it is key to decompose recombination gradients as a function of exon rank and gene length (in exon number here) to obtain a proper characterization of recombination patterns and better understand the underlying mechanism. Typically, because short genes are more frequent than longer ones, the average gradient is biased towards the 5′ end, such that J-shaped patterns can appear as U-shaped, as in P. tremula or almost flat pattern can appear as a standard 5′ to 3′ gradient as in P. vularis.

The intensity and location of hotspots are globally responsible for the polarity of recombination within genes. Hotspots overlapping 5′ and 3′ ends of genes have a consistent influence on the average recombination rate. More importantly, recombination hotspots seem similar in size but their intensities vary between species (Fig. 1, supplementary fig. S3, Supplementary Material online). Hotspot position also shapes the recombination gradient towards the 5′ or 3′ end. Though the TSS is usually the preferential location of CO hotspots, the increasing gradient towards the 3′ end in P. tremula, M. sieversii and C. sinensis is associated with the location of a significant fraction of recombination hotspots around the TTS (Choi and Henderson 2015). In other species, such as G. max, the intensity of hotspots compared to the genome-wide average recombination rate is rather low and TSS and TTS do not seem to be enriched in hotspots in G. max, which can explain why gradients are less pronounced in such species. Our modeling approach showed that these various patterns can be explained under a single framework taking hotspot location (around or away 5′ and 3′ flanking regions) and gene length and structure into account. This model fits all observed pattern well, even novel patterns undescribed until now, and provides a quantitative measure of the relative asymmetry between 5′ and 3′ hotspots. Note that in human, introns are extremely large compared to plants (supplementary fig. S9b, Supplementary Material online) and could explain why the 5′ gradient is steep and mostly concentrated in the first exon.

Gene length not only affects recombination in the middle of genes but it also lowers recombination for all exon positions. As the proportion of hotspots is invariant with gene length, the same amount of recombination is automatically spread over more exons, making average ρ/kb more or less negatively proportional to gene length. The phenomenological model we propose following Choi and Henderson (2015) can explain this observation simply because shorter genes have a closer influence of hotspots on both sides than longer genes (Fig. 7). A direct consequence is that shorter genes recombine more than longer genes on average, as we observed in 9 out of 12 species.

Finally, our results support the prediction that the similar GC gradients shaped by exon rank and gene class observed in plants can be generated by the underlying recombination gradients we characterized through the effect of gBGC (Ressayre et al. 2015). In A. thaliana and rice, recombination gradients per rank and gene class match GC gradients relatively well (compare our Fig. 6 with Fig. 5 in Ressayre et al. 2015).

LD-Based Estimates Are Reliable to Detect Gradients

We addressed the limitations and reliability of LD-based estimates of recombination rates. We carefully controlled for population structure, mating system and demography during the estimation procedure, following the most recent developments and benchmarks in the literature (Kamm et al. 2016; Raynaud et al. 2023). A recent simulation-based study showed that LDhat was the most accurate method to estimate fine-scale estimates of recombination rates from polymorphism data in presence of hotspots (Dutheil 2024). LDhat was even more accurate when it was combined with LDpop (Kamm et al. 2016). In A. thaliana LD-based estimates correlated well with the Rowan et al. (2019) pedigree-based CO map. The power to detect a gradient was not affected by the mating system or data quality (SNP density). The highly selfing A. thaliana is one of the most resolved gradients. We observed a clear gradient in C. sinensis despite quite low SNP density (0.6 SNP/kb) while the gradient was noisy in T. aestivum with five times more SNPs (3 SNP/kb). We also checked that recombination gradients were not biased by local patterns of polymorphism. Altogether, the extremely fine resolution we achieved in our analyses proves the good accuracy and reliability of our pipeline, which is freely available and easy to use (https://github.com/ThomasBrazier/ldhat-recombination-pipeline.git).

Conclusion

It was known that many organisms, including plants, had recombination hotspots (Mézard 2006; Buard and de Massy 2007; Choi and Henderson 2015), sometimes inducing recombination gradients in genic regions. Our results confirm that it is likely a rather general feature of plant genome but we also discovered an unexpected diversity of patterns. Indirect approaches like ours are important to detect unseen patterns and open new questions for future research, such as the mechanistic bases behind the 3′-driven gradient in P. tremula and C. sinensis.

We also provided a new model taking gene structure into account to unravel the complexity of these recombination gradients that likely emerged from both the location of recombination hotspots and the distribution of genes per class of exon numbers. Our modeling approach allowed us to highlight the relative importance of hotspots in 5′ and 3′ for the shape and intensity of recombination gradients.

We suggest that a detailed characterization of recombination and its known covariates (e.g. methylation, chromatin state) at such a fine scale will be helpful to better understand the evolution of recombination hotspots themselves (Dluzewska et al. 2018; Lloyd 2022), patterns of base composition in genic regions (Glémin et al. 2014), or the evolution of intron structure (Duret 2001). Moreover, the evolutionary significance of these recombination gradients in genic regions may have been underappreciated so far, as they scale with effective population size over thousands of generations and might have a deep impact on CDSs through complex interactions with linked selection or gBGC (Loewe and Charlesworth 2007; Glémin et al. 2014).

Materials and Methods

To produce fine-scale recombination maps we analyzed previously published polymorphism datasets in 11 plant species and human (Table 1) starting from Variant Call Format (VCF) files. We identified datasets by literature search with the keywords “resequencing,” “high-density SNP,” “polymorphism data,” “variant database,” “Whole Genome Sequencing,” or “WGS” in association with “plant,” “angiosperms” and a set of plant species names. We also searched directly for public genomic databases (e.g. 1,001 genomes). We used a previous study of 57 plant species to identify species with interesting and contrasting characteristics, such as monocots/eudicots, small/large genomes, low/high recombination rates, and homogeneous/heterogeneous recombination landscapes (Brazier and Glémin 2022). Before downstream analyses, we kept only biallelic SNPs and filtered them with the same quality criteria for all datasets (minor allele frequency >0.05, <10% of missing data per site, quality score ≥30 when the information was available). When the quality score was not available in the published VCF, we checked in the original study that the dataset we reused was already filtered by a similar quality score.

Genome assemblies and corresponding annotations (Generic Feature Format files, GFF) on which SNPs were originally mapped were retrieved from public databases, mostly NCBI (supplementary table S1, Supplementary Material online). GFFs with gene, mRNA, CDS and exon features were parsed using custom Python scripts to infer introns, 5′ and 3′ UTRs, flanking sequences (1 to 3 kb upstream/downstream, in consecutive windows of 1 kb) and intergenic sequences according to the GFF3 specifications. Scripts are freely available as a Python package (https://github.com/ThomasBrazier/PiSlice.git v0.2). Only protein-coding genes were retained when this information was annotated. For datasets without information about gene biotypes, we considered all genes as protein-coding. To avoid low sample size for longer genes, we removed genes with more than 14 exons as in Ressayre et al. (2015). When more than one splicing variant has been annotated, the splicing variant with the longest CDS sequence was retained, as it is often also the most expressed sequence. We considered only the CDS part of exons.

All statistical analyses were performed with R version 4.3.3 (R Core Team 2022). Genomic coordinates were manipulated with the R package GenomicRanges (Lawrence et al. 2013).

We developed a method for estimating robust LD-based recombination rates and implemented it in a custom pipeline performing data processing, SNP filtering, statistical phasing, controlling for population structure and demography, dealing with selfing and highly homozygous species, estimating recombination rates with LDhat 2.2 (Auton and McVean 2007; Auton et al. 2012), and inferring hotspots with LDhot (Auton et al. 2014). This pipeline was optimized to run large datasets with millions of SNPs. The Snakemake pipeline is freely available (https://github.com/ThomasBrazier/ldhat-recombination-pipeline.git v1.1) and the method is presented in detail below.

Population Sampling

Population-scaled recombination rates and hotspot detection are based on the assumption of a population without genetic structure and absence of gene flow from other populations (Dapper and Payseur 2018; Samuk and Noor 2022). For each dataset, we inferred population structure and sampled one single genetic population as homogeneous as possible and well distinct from other possible genetic clusters. As much as possible, we tried to sample a single large population (>40 diploid individuals or 80 haploid genomes for selfing species) of wild/landrace individuals with the highest possible polymorphism level. Results and metadata of the original publications were used for a presampling. The genetic structure within the presampled population was assessed with FastSTRUCTURE for K=1textto7 (Raj et al. 2014). FastSTRUCTURE was run on a random subset of 100,000 SNPs for computation issues. The best model (best value of K) was inferred by FastSTRUCTURE’s algorithm for multiple choices of K and validated by the visual assessment of admixture plots. If a population substructure was positively assessed, individuals were clustered into different genetic populations according to their mean admixture proportion. Clusters with higher proportions of missing data and low polymorphism (heterozygosity and θ^π) were discarded and we kept individuals in the largest remaining genetic cluster.

Population-Scaled Recombination Rates

After population sampling, the VCF was phased with Shapeit2, with a phasing window size of 2 Mb (O’Connell et al. 2014). Selfing species were treated as haploid individuals, by randomly sampling only one phased haplotype per individual. When the sample size was sufficient, at most 40 diploid individuals (or 80 haploid sequences for selfing species) were randomly subset for LDhat 2.2 for computational issues (see supplementary table S1, Supplementary Material online for sampling sizes).

The population-scaled recombination rates ρ/kb were estimated with LDhat 2.2 (Auton and McVean 2007; Auton et al. 2012). For the estimate of genetic diversity required by “LDhat” we used the genome averaged θ^π estimated with “VCFtools” on the population dataset (Danecek et al. 2011). The per-base mutation rate was set to 1.28×10−8 for all species. Previously to “LDhat,” we estimated the demography of the population per chromosome with “SMC++” (Terhorst et al. 2017) with the default parameters, including eight epochs, and a missing cutoff for runs of homozygosity adjusted per chromosome (see supplementary table S1, Supplementary Material online for parameters). We used “LDpop” to generate a complete demography-aware look-up table with the θ estimate and the fast approximate method (Kamm et al. 2016). The maximum value of ρ was set to 100, as suggested in the “LDhat” manual. Recombination rates between pairwise adjacent SNPs were estimated with the “LDhat interval” program, using different values of block penalty (bpen=5, 15, 25). Since the quality of recombination landscapes was similar between different block penalties, despite a smoothing effect for higher values, we kept estimates with the smaller block penalty (bpen=5). For computational efficiency of the “LDhat interval” part in larger datasets, the dataset was split into intervals of 2,000 SNPs (50 SNPs overlapping at each end), run in parallel and merged at the end of “LDhat interval.” We ran the Markov chain Monte carlo (MCMC) algorithm for 10,000,000 iterations, sampling every 5,000 iterations, with a burn-in of 1,000,000 first iterations. In a few species we adjusted these hyper-parameters to reach convergence (supplementary table S4, Supplementary Material online). The convergence of MCMC sampling chains was checked with “LDhat” diagnostic plots. For species with an already published Marey map in Brazier and Glémin (2022), the visual congruence between the LD map and the Marey map was used as a qualitative validation of our estimates.

Hotspot Detection

Recombination hotspots were inferred with “LDhot” using default settings (Auton et al. 2014). We tested the significance of putative hotspots in 3 kb windows, with a step size of 1 kb along the genome. The background window was ±50 kb around the hotspot center. The hotspots inferred were then summarized using the “LDhot” summary program with a significance cutoff of 0.001 for calling a hotspot and 0.01 for merging adjacent hotspots. In order to control for hotspot quality and reduce false positive hotspots, we compared the raw hotspot dataset to filtered datasets after soft and hard filtering strategies. As hotspots are generally defined as narrow genomic regions, the soft filtering approach removed every hotspot larger than 10 kb. The intensity of a hotspot was estimated by dividing the recombination peak rate by the background recombination rate (±50 kb around the hotspot center). Low-intensity hotspots were considered to be most likely false positive calls. In addition, we considered that extreme values of hotspot intensity were most probably artifactual due to high variance in the LDhat recombination map caused by assembly and/or SNP calling errors. The hard filtering removed hotspots larger than 10 kb or with an intensity lower than 4 or higher than 200.

Recombination Rates and Hotspots in Genomic Features

We calculated the median and average ρ/kb for each genomic feature with the median and the weighted arithmetic mean of the ρ/kb estimates intersecting this genomic interval (R “median” and “weighted.mean” function). The weights were the number of nucleotides spanned by each interval. The number of hotspots overlapping a given genomic feature was counted with “countOverlap” and “findOverlap” from the GenomicRanges R package (Lawrence et al. 2013). The 95% CIs of the mean and the median of ρ/kb within genomic features were estimated by bootstrap (1,000 iterations).

Gradients of Recombination

To assess a putative gradient of recombination at the 5′ end of genes, we estimated LD-based recombination rates (ρ/kb) in consecutive intervals (200 bp) along a region ±5 kb around the TSS, ATG codon (start of the CDS), and TTS. Annotations of start codons are generally more consistent and reliable than TSS among genomes that have been annotated using different methodologies. Intervals overlapping another gene were discarded. We also separated the average gradient from exon/intron-specific gradients. Each interval in a gene was marked as exon-specific (intron-specific, respectively) only if an exon (intron) overlapped at least 70% of the interval. A random control of ρ/kb was calculated by resampling the same intervals at a random position in the genome.

We also produced gradients of recombination by estimating recombination rates as a function of exon/intron ordinal rank within genes as in Ressayre et al. (2015). We computed the mean, weighted mean and median ρ/kb for each exon (CDS part) and intron. The recombination rate of exons/introns was averaged by exon/intron ordinal rank for the average gradient, and by rank and gene size (number of exons) for detailed gradients. As annotation quality and sample size declined for longer genes, we removed all genes longer than 14 exons or 10 kb.

We evaluated the performance of our method to produce gradients by testing whether it had the power to detect a true gradient of recombination or the capacity to reject an artifactual gradient. Firstly, we assessed whether the gradient observed could be an artifact. LDhat recombination rates are estimated within intervals between adjacent SNPs. As a consequence, a lack of SNPs at the beginning of genes where recombination hotspots are supposed to be located could reduce our resolution to accurately position hotspots in this region, thus leading to an artificial leakage of high ρ/kb values at a larger scale than the true biological size of the hotspot. To evaluate this, we recalculated the gradient of ρ/kb as a function of the rank, but this time we kept only ρ/kb estimates beginning in the genomic interval (exon/intron). Thus, we were sure to measure only the actual recombination within the interval. We observed the same patterns using this restricted set of ρ values confirming that the observed gradients were a true biological signal (supplementary fig. S13a,b, Supplementary Material online). Secondly, we assessed the power to detect a true gradient, if there was one, for a given empirical dataset. For each gene, the ranks of exons were sorted by descending ρ/kb values, to simulate the maximal negative 5′ to 3′ gradient that could be observed from empirical observations, if there is one. Ten out of eleven datasets had the necessary power to detect true gradients of recombination, except for T. aestivum lacking resolution (supplementary fig. S13b, Supplementary Material online).

Fitting of the Hotspot Model

We fitted a phenomenological model to recombination gradients explicitly taking into account exon positions (starting at xs and ending at xe), exon length (Lexon=xe−xs), and gene length (Ltot). We considered that recombination within gene was generated by two exponentially decreasing gradients starting from 5′ and 3′ ends, respectively, plus a relative basal rate r0, which can be negative. We assumed a same decreasing rate, a, for each end. If hotspots were centered exactly at the beginning and at the end of genes, 1/a could be interpreted as the length of the recombination tract. However, depending of the species, recombination can start to decrease before of just after gene boundaries. Accordingly, recombination rates in 5′ (r5) and 3′ (r3) do not necessarily correspond to hotspot rates. More precisely, r5=r5,maxe−ad5, where r5,max is the recombination rate at the center of the hotspot and d5 is the distance between the hotspot center and the start of the gene (similar interpretation for r3). Unfortunately, r5 and d5 (respectively, r3 and d3) are not identifiable. For each exon, the predicted recombination rate is thus given by

(1) rpredicted=r0+R1+R2Lexon

with

(2) R1=r5e−axs−e−axea

and

(3) R2=r3e−a(Ltot−xe)−e−a(Ltot−xs)a

The model was fitted with nonlinear least squares implemented in the “nls” R function which returned predictions and standard errors. The goodness of fit was assessed by plotting observed vs. predicted values.

Gradients of Polymorphism

To carefully account for a correlation of recombination rates with levels of polymorphism we estimated both SNP density (SNP/kb) and θ^π in each exon rank and gene classes in a similar manner as for recombination gradients. SNP density was the number of SNPs overlapping a given genomic feature found with “findOverlap” from the GenomicRanges R package (Lawrence et al. 2013). SNP density was averaged by calculating the mean per exon rank and gene classes. We also estimated θ^π per site with “VCFtools” on the final VCF dataset used for “LDhat” (Danecek et al. 2011). The average θ^π per exon rank and gene class was calculated as the sum of all θ^π overlapping a given exon rank divided by the total number of sites covered by these exons.

Supplementary Material

msae183_Supplementary_Data

Acknowledgments

We wish to thank Laurent Duret, Pierre-Alexandre Gagnaire, Marie Raynaud, Julien Joseph, Nicolas Lartillot, and all the other members of the HotRec ANR project for insighful discussions. Elise Rolland helped with the pipeline. We used the GENOUEST computing facility for bio-informatic analyses.

Supplementary Material

Supplementary material is available at Molecular Biology and Evolution online.

Author Contributions

S.G. originated the idea; T.B. analyzed the data; T.B. wrote the first draft; and T.B. and S.G. edited and revised the manuscript.

Funding

This work was supported by the Agence Nationale de la Recherche, Grant ANR-19-CE12472 0019/HotRec.

Data Availability

Analyses scripts and documentation are available at https://github.com/ThomasBrazier/landrec-gradients. Data produced in this study are available at https://osf.io/3aekw/.
==== Refs
References

1001 Genomes Consortium . Electronic address: magnus.nordborg@gmi.oeaw.ac.at and 1001 Genomes Consortium. 1,135 genomes reveal the global pattern of polymorphism in Arabidopsis thaliana. Cell. 2016:166 (2 ):481–491. 10.1016/j.cell.2016.05.063.27293186
Auton  A, Fledel-Alon  A, Pfeifer  S, Venn  O, Ségurel  L, Street  T, Leffler  EM, Bowden  R, Aneas  I, Broxholme  J, et al  A fine-scale chimpanzee genetic map from population sequencing. Science. 2012:336 (6078 ):193–198. 10.1126/science.1216872.22422862
Auton  A, McVean  G. Recombination rate estimation in the presence of hotspots. Genome Res. 2007:17 (8 ):1219–1227. 10.1101/gr.6386707.17623807
Auton  A, Myers  S, McVean  G. 2014. Identifying recombination hotspots using population genetic data, arXiv, arXiv:1403.4264 [q-bio], preprint: not peer reviewed. 10.48550/arXiv.1403.4264
Auton  A, Rui Li  Y, Kidd  J, Oliveira  K, Nadel  J, Holloway  JK, Hayward  JJ, Cohen  PE, Greally  JM, Wang  J, et al  Genetic recombination is targeted towards gene promoter regions in dogs. PLoS Genet. 2013a:9 (12 ):e1003984. 10.1371/journal.pgen.1003984.24348265
Baker  Z, Schumer  M, Haba  Y, Bashkirova  L, Holland  C, Rosenthal  GG, Przeworski  M. Repeated losses of PRDM9-directed recombination despite the conservation of PRDM9 across vertebrates. eLife. 2017:6 :e24133. 10.7554/eLife.24133.28590247
Barroso  GV, Dutheil  JY. The landscape of nucleotide diversity in Drosophila melanogaster is shaped by mutation rate variation. Peer Community J. 2023:3 :e40. 10.24072/pcjournal.267.
Baudat  F, Buard  J, Grey  C, Fledel-Alon  A, Ober  C, Przeworski  M, Coop  G, de Massy  B. PRDM9 is a major determinant of meiotic recombination hotspots in humans and mice. Science. 2010:327 (5967 ):836–840. 10.1126/science.1183439.20044539
Blat  Y, Protacio  RU, Hunter  N, Kleckner  N. Physical and functional interactions among basic chromosome organizational features govern early steps of meiotic chiasma formation. Cell. 2002:111 (6 ):791–802. 10.1016/S0092-8674(02)01167-4.12526806
Boulton  A, Myers  RS, Redfield  RJ. The hotspot conversion paradox and the evolution of meiotic recombination. Proc Natl Acad Sci USA. 1997:94 (15 ):8058–8063. 10.1073/pnas.94.15.8058.9223314
Brazier  T, Glémin  S. Diversity and determinants of recombination landscapes in flowering plants. PLoS Genet. 2022:18 (8 ):e1010141. 10.1371/journal.pgen.1010141.36040927
Brick  K, Smagulova  F, Khil  P, Camerini-Otero  RD, Petukhova  GV. Genetic recombination is directed away from functional genomic elements in mice. Nature. 2012:485 (7400 ):642–645. 10.1038/nature11089.22660327
Buard  J, de Massy  B. Playing hide and seek with mammalian meiotic crossover hotspots. Trends Genet. 2007:23 (6 ):301–309. 10.1016/j.tig.2007.03.014.17434233
Cai  X, Sun  X, Xu  C, Sun  H, Wang  X, Ge  C, Zhang  Z, Wang  Q, Fei  Z, Jiao  C, et al  Genomic analyses provide insights into spinach domestication and the genetic basis of agronomic traits. Nat Commun. 2021:12 (1 ):7246. 10.1038/s41467-021-27432-z.34903739
Choi  K, Henderson  IR. Meiotic recombination hotspots - a comparative view. Plant J. 2015:83 (1 ):52–61. 10.1111/tpj.2015.83.issue-1.25925869
Choi  K, Zhao  X, Kelly  KA, Venn  O, Higgins  JD, Yelina  NE, Hardcastle  TJ, Ziolkowski  PA, Copenhaver  GP, Franklin  FCH, et al  Arabidopsis meiotic crossover hot spots overlap with H2A.Z nucleosomes at gene promoters. Nat Genet. 2013:45 (11 ):1327–1336. 10.1038/ng.2766.24056716
Clément  Y, Sarah  G, Holtz  Y, Homa  F, Pointet  S, Contreras  S, Nabholz  B, Sabot  F, Sauné  L, Ardisson  M, et al  Evolutionary forces affecting synonymous variations in plant genomes. PLoS Genet. 2017:13 (5 ):e1006799. 10.1371/journal.pgen.1006799.28531201
Cooper  TJ, Garcia  V, Neale  MJ. Meiotic DSB patterning: a multifaceted process. Cell Cycle. 2016:15 (1 ):13–21. 10.1080/15384101.2015.1093709.26730703
Danecek  P, Auton  A, Abecasis  G, Albers  CA, Banks  E, DePristo  MA, Handsaker  RE, Lunter  G, Marth  GT, Sherry  ST, et al  The variant call format and VCFtools. Bioinformatics. 2011:27 (15 ):2156–2158. 10.1093/bioinformatics/btr330.21653522
Dapper  AL, Payseur  BA. Effects of demographic history on the detection of recombination hotspots from linkage disequilibrium. Mol Biol Evol. 2018:35 (2 ):335–353. 10.1093/molbev/msx272.29045724
de Massy  B . Initiation of meiotic recombination: how and where? conservation and specificities among eukaryotes. Annu Rev Genet. 2013:47 (1 ):563–599. 10.1146/genet.2013.47.issue-1.24050176
Detloff  P, White  MA, Petes  TD. Analysis of a gene conversion gradient at the HIS4 locus in Saccharomyces cerevisiae. Genetics. 1992:132 (1 ):113–123. 10.1093/genetics/132.1.113.1398048
Dluzewska  J, Szymanska  M, Ziolkowski  PA. Where to cross over? defining crossover sites in plants. Front Genet. 2018:9 :609. 10.3389/fgene.2018.00609.30619450
Dooner  HK, Martínez-Férez  IM. Recombination occurs uniformly within the bronze gene, a meiotic recombination hotspot in the maize genome. Plant Cell. 1997:9 (9 ):1633–1646. 10.1105/tpc.9.9.1633.9338965
Duret  L . Why do genes have introns? Recombination might add a new piece to the puzzle. Trends Genet. 2001:17 (4 ):172–175. 10.1016/S0168-9525(01)02236-3.11275306
Dutheil  JY . On the estimation of genome-average recombination rates. Genetics. 2024:227 (2 ):iyae051. 10.1093/genetics/iyae051.38565705
Fridman  E, Pleban  T, Zamir  D. A recombination hotspot delimits a wild-species quantitative trait locus for tomato sugar content to 484 bp within an invertase gene. Proc Natl Acad Sci USA. 2000:97 (9 ):4718–4723. 10.1073/pnas.97.9.4718.10781077
Glémin  S, Clément  Y, David  J, Ressayre  A. GC content evolution in coding regions of angiosperm genomes: a unifying hypothesis. Trends Genet. 2014:30 (7 ):263–270. 10.1016/j.tig.2014.05.002.24916172
Grey  C, Baudat  F, de Massy  B. PRDM9, a driver of the genetic map. PLoS Genet. 2018:14 (8 ):e1007479. 10.1371/journal.pgen.1007479.30161134
Guo  S, Zhao  S, Sun  H, Wang  X, Wu  S, Lin  T, Ren  Y, Gao  L, Deng  Y, Zhang  J, et al  Resequencing of 414 cultivated and wild watermelon accessions identifies selection for fruit quality traits. Nat Genet. 2019:51 (11 ):1616–1623. 10.1038/s41588-019-0518-4.31676863
Haenel  Q, Laurentino  TG, Roesti  M, Berner  D. Meta-analysis of chromosome-scale crossover rate variation in eukaryotes and its significance to evolutionary genomics. Mol Ecol. 2018:27 (11 ):2477–2497. 10.1111/mec.2018.27.issue-11.29676042
He  Y, Wang  M, Dukowic-Schulze  S, Zhou  A, Tiang  C-L, Shilo  S, Sidhu  GK, Eichten  S, Bradbury  P, Springer  NM, et al  Genomic features shaping the landscape of meiotic double-strand-break hotspots in maize. Proc Natl Acad Sci USA. 2017:114 (46 ):12231–12236. 10.1073/pnas.1713225114.29087335
Hellsten  U, Wright  KM, Jenkins  J, Shu  S, Yuan  Y, Wessler  SR, Schmutz  J, Willis  JH, Rokhsar  DS. Fine-scale variation in meiotic recombination in Mimulus inferred from population shotgun sequencing. Proc Natl Acad Sci USA. 2013:110 (48 ):19478–19482. 10.1073/pnas.1319032110.24225854
Hoge  CR, De Manuel  M, Mahgoub  M, Okami  N, Fuller  ZL, Banerjee  S, Baker  Z, Mcnulty  M, Andolfatto  P, Macfarlan  TS, et al  Patterns of recombination in snakes reveal a tug of war between PRDM9 and promoter-like features. Science. 2024:383 (6685 ):eadj7026. 10.1126/science.adj7026.38386752
Jessop  L, Allers  T, Lichten  M. Infrequent co-conversion of markers flanking a meiotic recombination initiation site in Saccharomyces cerevisiae. Genetics. 2005:169 (3 ):1353–1367. 10.1534/genetics.104.036509.15654098
Joseph  J, Prentout  D, Laverré  A, Tricou  T, Duret  L. High prevalence of Prdm9-independent recombination hotspots in placental mammals. Proc Natl Acad Sci USA. 2024:121 (23 ):e2401973121. 10.1073/pnas.2401973121.38809707
Kamm  JA, Spence  JP, Chan  J, Song  YS. Two-locus likelihoods under variable population size and fine-scale recombination rate estimation. Genetics. 2016:203 (3 ):1381–1399. 10.1534/genetics.115.184820.27182948
Kauppi  L, Jeffreys  AJ, Keeney  S. Where the crossovers are: recombination distributions in mammals. Nat Rev Genet. 2004:5 (6 ):413–424. 10.1038/nrg1346.15153994
Kianian  PMA, Wang  M, Simons  K, Ghavami  F, He  Y, Dukowic-Schulze  S, Sundararajan  A, Sun  Q, Pillardy  J, Mudge  J, et al  High-resolution crossover mapping reveals similarities and differences of male and female recombination in maize. Nat Commun. 2018:9 (1 ):2370. 10.1038/s41467-018-04562-5.29915302
Lam  I, Keeney  S. Nonparadoxical evolutionary stability of the recombination initiation landscape in yeast. Science. 2015:350 (6263 ):932–937. 10.1126/science.aad0814.26586758
Latrille  T, Duret  L, Lartillot  N. The Red Queen model of recombination hot-spot evolution: a theoretical investigation. Philos Trans R Soc B Biol Sci. 2017:372 (1736 ):20160463. 10.1098/rstb.2016.0463.
Lawrence  M, Huber  W, Pagès  H, Aboyoun  P, Carlson  M, Gentleman  R, Morgan  MT, Carey  VJ. Software for computing and annotating genomic ranges. PLoS Comput Biol. 2013:9 (8 ):e1003118. 10.1371/journal.pcbi.1003118.23950696
Li  X, Li  L, Yan  J. Dissecting meiotic recombination based on tetrad analysis by single-microspore sequencing in maize. Nat Commun. 2015:6 (1 ):6648. 10.1038/ncomms7648.25800954
Liu  S, Zhang  L, Sang  Y, Lai  Q, Zhang  X, Jia  C, Long  Z, Wu  J, Ma  T, Mao  K, et al  Demographic history and natural selection shape patterns of deleterious mutation load and barriers to introgression across Populus genome. Mol Biol Evol. 2022:39 (2 ):msac008. 10.1093/molbev/msac008.35022759
Lloyd  A . Crossover patterning in plants. Plant Reprod. 2022:36 (1 ):55–72. 10.1007/s00497-022-00445-4.35834006
Loewe  L, Charlesworth  B. Background selection in single genes may explain patterns of codon bias. Genetics. 2007:175 (3 ):1381–1393. 10.1534/genetics.106.065557.17194784
Lozano  R, Gazave  E, Stetter  M, Valluru  R, Bandillo  N, Fernandes  SB, Brown  PJ, Shakoor  N, Mockler  TC, Ross-Ibarra  J, et al  Comparative evolutionary genetics of deleterious load in sorghum and maize. Nature Plants. 2021:7 (1 ):17–24. 10.1038/s41477-020-00834-5.33452486
Malone  RE, Bullard  S, Lundquist  S, Kim  S, Tarkowski  T. A meiotic gene conversion gradient opposite to the direction of transcription. Nature. 1992:359 (6391 ):154–155. 10.1038/359154a0.1355857
Marand  AP, Jansky  SH, Zhao  H, Leisner  CP, Zhu  X, Zeng  Z, Crisovan  E, Newton  L, Hamernik  AJ, Veilleux  RE, et al  Meiotic crossovers are associated with open chromatin and enriched with Stowaway transposons in potato. Genome Biol. 2017:18 (1 ):203. 10.1186/s13059-017-1326-8.29084572
Marand  AP, Zhao  H, Zhang  W, Zeng  Z, Fang  C, Jiang  J. Historical meiotic crossover hotspots fueled patterns of evolutionary divergence in rice. Plant Cell. 2019a:31 (3 ):645–662. 10.1105/tpc.18.00750.30705136
Mézard  C . Meiotic recombination hotspots in plants. Biochem Soc Trans. 2006:34 (4 ):531–534. 10.1042/BST0340531.16856852
Myers  S . A fine-scale map of recombination rates and hotspots across the human genome. Science. 2005:310 (5746 ):321–324. 10.1126/science.1117196.16224025
O’Connell  J, Gurdasani  D, Delaneau  O, Pirastu  N, Ulivi  S, Cocca  M, Traglia  M, Huang  J, Huffman  JE, Rudan  I, et al  A general approach for haplotype phasing across the full spectrum of relatedness. PLoS Genet. 2014:10 (4 ):e1004234. 10.1371/journal.pgen.1004234.24743097
O’Reilly  PF, Birney  E, Balding  DJ. Confounding between recombination and selection, and the Ped/Pop method for detecting selection. Genome Res. 2008:18 (8 ):1304–1313. 10.1101/gr.067181.107.18617692
Pan  J, Sasaki  M, Kniewel  R, Murakami  H, Blitzblau  HG, Tischfield  SE, Zhu  X, Neale  MJ, Jasin  M, Socci  ND, et al  A hierarchical combination of factors shapes the genome-wide topography of yeast meiotic recombination initiation. Cell. 2011:144 (5 ):719–731. 10.1016/j.cell.2011.02.009.21376234
Raj  A, Stephens  M, Pritchard  JK. fastSTRUCTURE: variational inference of population structure in large SNP data sets. Genetics. 2014:197 (2 ):573–589. 10.1534/genetics.114.164350.24700103
Raynaud  M, Gagnaire  P-A, Galtier  N. Performance and limitations of linkage-disequilibrium-based methods for inferring the genomic landscape of recombination and detecting hotspots: a simulation study. Peer Community J. 2023:3 :e27. 10.24072/pcjournal.254.
Raynaud  M, Sanna  P, Joseph  J, Clément  J, Imai  Y, Lareyre  J-J, Laurent  A, Galtier  N, Baudat  F, Duret  L, et al  2024. PRDM9 drives the location and rapid evolution of recombination hotspots in salmonids, bioRxiv, preprint: not peer reviewed. 10.1101/2024.03.06.583651.
R Core Team . 2022. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing.
Ressayre  A, Glémin  S, Montalent  P, Serre-Giardi  L, Dillmann  C, Joets  J. Introns structure patterns of variation in nucleotide composition in Arabidopsis thaliana and rice protein-coding genes. Genome Biol Evol. 2015:7 (10 ):2913–2928. 10.1093/gbe/evv189.26450849
Rowan  BA, Heavens  D, Feuerborn  TR, Tock  AJ, Henderson  IR, Weigel  D. An ultra high-density Arabidopsis thaliana crossover map that refines the influences of structural variation and epigenetic features. Genetics. 2019:213 (3 ):771–787. 10.1534/genetics.119.302406.31527048
Samuk  K, Noor  MAF. Gene flow biases population genetic inference of recombination rate. G3. 2022:12 (11 ):jkac236. 10.1093/g3journal/jkac236.36103705
Schultes  NP, Szostak  JW. Decreasing gradients of gene conversion on both sides of the initiation site for meiotic recombination at the ARG4 locus in yeast. Genetics. 1990:126 (4 ):813–822. 10.1093/genetics/126.4.813.1981763
Serres-Giardi  L, Belkhir  K, David  J, Glémin  S. Patterns and evolution of nucleotide landscapes in seed plants. Plant Cell. 2012:24 (4 ):1379–1397. 10.1105/tpc.111.093674.22492812
Shilo  S, Melamed-Bessudo  C, Dorone  Y, Barkai  N, Levy  AA. DNA crossover motifs associated with epigenetic modifications delineate open chromatin regions in Arabidopsis. Plant Cell. 2015:27 (9 ):2427–2436. 10.1105/tpc.15.00391.26381163
Singhal  S, Leffler  EM, Sannareddy  K, Turner  I, Venn  O, Hooper  DM, Strand  AI, Li  Q, Raney  B, Balakrishnan  CN, et al  Stable recombination hotspots in birds. Science. 2015:350 (6263 ):928–932. 10.1126/science.aad0843.26586757
Slavov  GT, DiFazio  SP, Martin  J, Schackwitz  W, Muchero  W, Rodgers-Melnick  E, Lipphardt  MF, Pennacchio  CP, Hellsten  U, Pennacchio  LA, et al Genome resequencing reveals multiscale geographic structure and extensive linkage disequilibrium in the forest tree Populus trichocarpa. New Phytol. 2012:196 (3 ):713–725. 10.1111/nph.2012.196.issue-3.22861491
Smagulova  F, Gregoretti  IV, Brick  K, Khil  P, Camerini-Otero  RD, Petukhova  GV. Genome-wide analysis reveals novel molecular features of mouse recombination hotspots. Nature. 2011:472 (7343 ):375–378. 10.1038/nature09869.21460839
Song  Q-X, Lu  X, Li  Q-T, Chen  H, Hu  X-Y, Ma  B, Zhang  W-K, Chen  S-Y, Zhang  J-S. Genome-wide analysis of DNA methylation in soybean. Mol Plant. 2013:6 (6 ):1961–1974. 10.1093/mp/sst123.23966636
Stukenbrock  EH, Dutheil  JY. Fine-scale recombination maps of fungal plant pathogens reveal dynamic recombination landscapes and intragenic hotspots. Genetics. 2018:208 (3 ):1209–1229. 10.1534/genetics.117.300502.29263029
Sudmant  PH, Mallick  S, Nelson  BJ, Hormozdiari  F, Krumm  N, Huddleston  J, Coe  BP, Baker  C, Nordenfelt  S, Bamshad  M,et al Global diversity, population stratification, and selection of human copy-number variation. Science. 2015:349 (6253 ):aab3761. 10.1126/science.aab3761.26249230
Sun  X, Jiao  C, Schwaninger  H, Chao  CT, Ma  Y, Duan  N, Khan  A, Ban  S, Xu  K, Cheng  L, et al  Phased diploid genome assemblies and pan-genomes provide insights into the genetic history of apple domestication. Nat Genet. 2020:52 (12 ):1423–1432. 10.1038/s41588-020-00723-9.33139952
Terhorst  J, Kamm  JA, Song  YS. Robust and scalable inference of population history from hundreds of unphased whole genomes. Nat Genet. 2017:49 (2 ):303–309. 10.1038/ng.3748.28024154
Tock  AJ, Henderson  IR. Hotspots for initiation of meiotic recombination. Front Genet. 2018:9 :521. 10.3389/fgene.2018.00521.30467513
Vining  KJ, Pomraning  KR, Wilhelm  LJ, Priest  HD, Pellegrini  M, Mockler  TC, Freitag  M, Strauss  SH. Dynamic DNA cytosine methylation in the Populus trichocarpa genome: tissue-level variation and relationship to gene expression. BMC Genomics. 2012:13 (1 ):27. 10.1186/1471-2164-13-27.22251412
Wang  W, Mauleon  R, Hu  Z, Chebotarov  D, Tai  S, Wu  Z, Li  M, Zheng  T, Fuentes  RR, Zhang  F, et al  Genomic variation in 3,010 diverse accessions of Asian cultivated rice. Nature. 2018:557 (7703 ):43–49. 10.1038/s41586-018-0063-9.29695866
Webster  MT, Hurst  LD. Direct and indirect consequences of meiotic recombination: implications for genome evolution. Trends Genet. 2012:28 (3 ):101–109. 10.1016/j.tig.2011.11.002.22154475
Wijnker  E, Velikkakam James  G, Ding  J, Becker  F, Klasen  JR, Rawat  V, Rowan  BA, de Jong  DF, de Snoo  CB, Zapata  L, et al  The genomic landscape of meiotic crossovers and gene conversions in Arabidopsis thaliana. eLife. 2013:2 :e01426. 10.7554/eLife.01426.24347547
Wu  J, Wang  L, Fu  J, Chen  J, Wei  S, Zhang  S, Zhang  J, Tang  Y, Chen  M, Zhu  J, et al Resequencing of 683 common bean genotypes identifies yield component trait associations across a north–south cline. Nat Genet. 2020:52 (1 ):118–125. 10.1038/s41588-019-0546-0.31873299
Yang  C, Yan  J, Jiang  S, Li  X, Min  H, Wang  X, Hao  D. Resequencing 250 soybean accessions: new insights into genes associated with agronomic traits and genetic networks. Genom Proteom Bioinform. 2021:20 (1 ):29–41. 10.1016/j.gpb.2021.02.009.
Yelina  NE, Lambing  C, Hardcastle  TJ, Zhao  X, Santos  B, Henderson  IR. DNA methylation epigenetically silences crossover hot spots and controls chromosomal domains of meiotic recombination in Arabidopsis. Genes Dev. 2015:29 (20 ):2183–2202. 10.1101/gad.270876.115.26494791
Zhang  X, Chen  S, Shi  L, Gong  D, Zhang  S, Zhao  Q, Zhan  D, Vasseur  L, Wang  Y, Yu  J, et al  Haplotype-resolved genome assembly provides insights into evolutionary history of the tea plant Camellia sinensis. Nat Genet. 2021:53 (8 ):1250–1259. 10.1038/s41588-021-00895-y.34267370
Zhou  Y, Zhao  X, Li  Y, Xu  J, Bi  A, Kang  L, Xu  D, Chen  H, Wang  Y, Wang  Y-g.  et al  Triticum population sequencing provides insights into wheat adaptation. Nat Genet. 2020:52 (12 ):1412–1422. 10.1038/s41588-020-00722-w.33106631
