Diversity in Recombination Hotspot Characteristics
and Gene Structure Shape Fine-Scale Recombination
Patterns in Plant Genomes

Thomas Brazier, ID <fnref n="1" />,*
Sylvain Glémin ID <fnref n="1" />,<fnref n="2" />

**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.

**Key words:** plant genomics, recombination hotspot, gene structure, recombination, crossovers.

**Introduction**

Meiotic recombination is a general feature of sexually repro-
ducing species. During meiosis, crossovers (COs) ensure
proper segregation of homologous chromosomes and re-
shuffle alleles from the two parental genomes. The distribu-
tion 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 recombin-
ation occurs have many implications because recombin-
ation 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 spe-
cific DNA motifs is evolving fast, thus changing recombin-
ation 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 rearrange-
ments (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

<footnote n="1">¹Unité Mixte de Recherche (UMR) 6553 - ECOBIO (Ecosystems, Biodiversity, Evolution), University of Rennes, CNRS, Rennes, France</footnote>

<footnote n="2">²Department of Ecology and Genetics, Evolutionary Biology Center and Science for Life Laboratory, Uppsala University, Uppsala, Sweden</footnote>

<footnote>*Corresponding author: E-mail: thomas.brazier@univ-rennes.fr.</footnote>

<footnote>Associate editor: Grace Yuh Chwen Lee</footnote>
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 *Caenorhabditis*; (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
<table><caption>Table 1 Summary statistics of 12 datasets (11 plants and 1 human)</caption><thead><tr><th>Species</th><th>Reference</th><th># SNPs</th><th>Mb</th><th>Map masked (%)</th><th>SNP/kb</th><th># Chr.</th><th># Diploid genomes</th><th>Mean ρ/kb</th><th>Median ρ/kb</th><th>θ̂<sub>π</sub></th><th>ρ/θ</th></tr></thead><tbody><tr><td>A. thaliana</td><td>The 1001 Genomes Consortium (2016)</td><td>641,453</td><td>119</td><td>4.7</td><td>5.39</td><td>5</td><td>40</td><td>2.68</td><td>1.21</td><td>0.0014</td><td>1.97</td></tr><tr><td>C. sinensis</td><td>Zhang et al. (2021)</td><td>1,808,673</td><td>2,985</td><td>42.3</td><td>0.6</td><td>15</td><td>40</td><td>4.64</td><td>4.60</td><td>0.0002</td><td>25.75</td></tr><tr><td>C. lanatus</td><td>Guo et al. (2019)</td><td>366,856</td><td>362</td><td>11.7</td><td>1.0</td><td>11</td><td>40</td><td>1.53</td><td>0.07</td><td>0.0004</td><td>4.12</td></tr><tr><td>G. max</td><td>Yang et al. (2021)</td><td>4,337,396</td><td>949</td><td>12.2</td><td>4.57</td><td>20</td><td>40</td><td>12.34</td><td>1.51</td><td>0.0012</td><td>10.38</td></tr><tr><td>H. sapiens</td><td>Sudmant et al. (2015)</td><td>6,584,165</td><td>2,796</td><td>4.7</td><td>2.3</td><td>22</td><td>40</td><td>3.37</td><td>0.57</td><td>0.0008</td><td>3.99</td></tr><tr><td>M. sieversii</td><td>Sun et al. (2020)</td><td>6,295,986</td><td>651</td><td>1</td><td>9.66</td><td>17</td><td>37</td><td>8.06</td><td>3.25</td><td>0.0036</td><td>2.21</td></tr><tr><td>O. sativa indica subsp.</td><td>Wang et al. (2018)</td><td>1,710,624</td><td>372</td><td>8.7</td><td>4.59</td><td>12</td><td>40</td><td>2.77</td><td>2.85</td><td>0.0021</td><td>1.34</td></tr><tr><td>P. vulgaris</td><td>Wu et al. (2020)</td><td>1,840,346</td><td>515</td><td>31.7</td><td>3.6</td><td>11</td><td>40</td><td>7.11</td><td>2.23</td><td>0.0011</td><td>6.69</td></tr><tr><td>P. tremula</td><td>Liu et al. (2022)</td><td>1,663,613</td><td>361</td><td>11.2</td><td>4.6</td><td>19</td><td>36</td><td>4.67</td><td>2.70</td><td>0.0011</td><td>4.24</td></tr><tr><td>S. bicolor</td><td>Lozano et al. (2021)</td><td>422,937</td><td>683</td><td>30.9</td><td>0.6</td><td>10</td><td>32</td><td>0.44</td><td>0.03</td><td>0.0002</td><td>2.77</td></tr><tr><td>S. oleracea</td><td>Cai et al. (2021)</td><td>3,283,574</td><td>879</td><td>0.5</td><td>3.7</td><td>6</td><td>40</td><td>2.68</td><td>1.52</td><td>0.0011</td><td>2.53</td></tr><tr><td>T. aestivum B subgenome</td><td>Zhou et al. (2020)</td><td>15,614,566</td><td>5,177</td><td>27.8</td><td>3.0</td><td>7</td><td>40</td><td>3.12</td><td>2.14</td><td>0.0004</td><td>7.69</td></tr></tbody></table>

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 ρ/θ.

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).

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

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.

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).

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.
<figure><img src="image_1.png" /><figcaption>**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.</figcaption></figure>

## 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.

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
<figure><img src="image_2.png" /><figcaption>**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*.</figcaption></figure>

<table><caption>**Table 3** Number of genic and intergenic hotspots, i.e. hotspots overlapping or not a gene feature (soft-filtered data)</caption><thead><tr><th>Species</th><th># Intergenic hotspots</th><th># Genic hotspots</th><th># Expected genic (95% CI)</th></tr></thead><tbody><tr><td><i>A. thaliana</i></td><td>215</td><td>2,470</td><td>655 (614 to 700)</td></tr><tr><td><i>C. sinensis</i></td><td>8,449</td><td>5,649</td><td>813 (760 to 867)</td></tr><tr><td><i>C. lanatus</i></td><td>2,629</td><td>1,882</td><td>687 (641 to 729)</td></tr><tr><td><i>G. max</i></td><td>50,076</td><td>16,211</td><td>21,144 (20,917 to 21,359)</td></tr><tr><td><i>H. sapiens</i></td><td>43,122</td><td>14,830</td><td>12,894 (12,703 to 13,091)</td></tr><tr><td><i>M. sieversii</i></td><td>7,414</td><td>6,400</td><td>3,101 (3,011 to 3,194)</td></tr><tr><td><i>O. sativa</i></td><td>2,263</td><td>3,871</td><td>1,147 (1,090 to 1,206)</td></tr><tr><td><i>P. vulgaris</i></td><td>4,308</td><td>3,265</td><td>1,348 (1,287 to 1,413)</td></tr><tr><td><i>P. tremula</i></td><td>1,356</td><td>5,443</td><td>1,463 (1,393 to 1,528)</td></tr><tr><td><i>S. bicolor</i></td><td>3,238</td><td>3,086</td><td>759 (712 to 807)</td></tr><tr><td><i>S. oleracea</i></td><td>58,161</td><td>15,270</td><td>23,883 (23,645 to 24,121)</td></tr><tr><td><i>T. aestivum</i></td><td>68,074</td><td>3,818</td><td>11,457 (11,267 to 11,631)</td></tr></tbody></table>

The expected number of genic hotspots was computed by 1,000 random shuffling of hotspots ranges.

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).
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

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.

## 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 increased 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).

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
<figure><img src="image_4.png" /><figcaption>**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* (*e*<sub>*i*</sub> points) or intron of rank *i* (*i*<sub>*i*</sub> 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.</figcaption></figure>

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 $\rho$ is the product of the effective population size $N_e$ and the CO rate $r$ ($\rho = 4N_e r$). As such $\rho$ estimates are potentially biased if underlying variations in $N_e$ 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 $\rho$ (see material and methods for details), fine-scale variations in $N_e$ 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 $\theta_\pi$) pooled in the same manner (Fig. 5a,b). We observed a weak noisy gradient of SNP density in *A. thaliana* but a clear
<figure><img src="image_5.png" /><figcaption>**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</figcaption></figure>

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).

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
<figure><img src="image_6.png" /><figcaption>**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.</figcaption></figure>

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).

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, $r_5$ and $r_3$, 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 ($r^2$ from 0.64 to 0.96, Table 4, Fig. 7b, c, supplementary fig. S12, Supplementary Material online) and the ratio $r_5/r_3$ was congruent with the observed gradients.

## 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.
<table><caption>**Table 4** Fit of the model and parameter estimates for 11 species</caption><thead><tr><th>Species</th><th>r<sub>0</sub></th><th>r<sub>0</sub> (SE)</th><th>r<sub>5</sub></th><th>r<sub>5</sub> (SE)</th><th>r<sub>3</sub></th><th>r<sub>3</sub> (SE)</th><th>Ratio r<sub>5</sub>/r<sub>3</sub></th><th>a</th><th>a (SE)</th><th>Tract length (bp)</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>A. thaliana</td><td>-0.74</td><td>0.09</td><td>1.43</td><td>0.05</td><td>1.17</td><td>0.05</td><td>1.23</td><td>0.000505</td><td>0.00003</td><td>1,982</td><td>0.95</td></tr><tr><td>C. sinensis</td><td>0.49</td><td>0.24</td><td>0.69</td><td>0.11</td><td>1.59</td><td>0.16</td><td>0.44</td><td>0.000128</td><td>0.000023</td><td>7,817</td><td>0.76</td></tr><tr><td>C. lanatus</td><td>-3.58</td><td>1.29</td><td>2.35</td><td>0.7</td><td>1.62</td><td>0.59</td><td>1.45</td><td>0.000042</td><td>0.000011</td><td>23,950</td><td>0.72</td></tr><tr><td>G. max</td><td>-0.01</td><td>0.04</td><td>0.35</td><td>0.02</td><td>0.28</td><td>0.02</td><td>1.26</td><td>0.000302</td><td>0.000031</td><td>3,315</td><td>0.84</td></tr><tr><td>H. sapiens</td><td>0.36</td><td>0.02</td><td>0.23</td><td>0.02</td><td>0.15</td><td>0.02</td><td>1.56</td><td>0.000067</td><td>0.000012</td><td>14,962</td><td>0.71</td></tr><tr><td>M. sieversii</td><td>0.77</td><td>0.07</td><td>0.66</td><td>0.06</td><td>0.99</td><td>0.06</td><td>0.67</td><td>0.000629</td><td>0.000076</td><td>1,590</td><td>0.78</td></tr><tr><td>O. sativa</td><td>-0.07</td><td>0.06</td><td>0.56</td><td>0.04</td><td>0.48</td><td>0.03</td><td>1.15</td><td>0.000294</td><td>0.000028</td><td>3,399</td><td>0.86</td></tr><tr><td>P. vulgaris</td><td>-0.08</td><td>0.09</td><td>0.48</td><td>0.05</td><td>0.46</td><td>0.05</td><td>1.04</td><td>0.000175</td><td>0.000026</td><td>5,712</td><td>0.84</td></tr><tr><td>P. tremula</td><td>-1.77</td><td>0.23</td><td>3.4</td><td>0.14</td><td>5.38</td><td>0.15</td><td>0.63</td><td>0.000442</td><td>0.000025</td><td>2,263</td><td>0.96</td></tr><tr><td>S. bicolor</td><td>-0.02</td><td>0.02</td><td>0.09</td><td>0.01</td><td>0.06</td><td>0.01</td><td>1.44</td><td>0.00028</td><td>0.000047</td><td>3,576</td><td>0.69</td></tr><tr><td>S. oleracea</td><td>-0.05</td><td>0.03</td><td>0.65</td><td>0.02</td><td>0.59</td><td>0.02</td><td>1.1</td><td>0.000239</td><td>0.000013</td><td>4,180</td><td>0.95</td></tr><tr><td>T. aestivum</td><td>-0.04</td><td>0.13</td><td>0.44</td><td>0.07</td><td>0.43</td><td>0.07</td><td>1.01</td><td>0.000186</td><td>0.000041</td><td>5,365</td><td>0.64</td></tr></tbody></table>

The basal recombination rate *r*<sub>0</sub> (with SE), the 5' recombination rate *r*<sub>5</sub> (with SE), the 3' recombination rate *r*<sub>3</sub> (with SE), the ratio *r*<sub>5</sub>/*r*<sub>3</sub>, the inverse of the tract length *a* (with SE), the tract length in bp, and the *R*<sup>2</sup> of the nonlinear least squares fitting procedure. The standard errors of the predictions were obtained with the "nls" R function.

<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>
## 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*<sup>-/-</sup> 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. vulgaris.

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; Guard 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* = 1 *textto* 7 (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 $\hat{\theta}_{\pi}$) 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 $\rho$/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 $\hat{\theta}_{\pi}$ estimated with “VCFtools” on the population dataset (Danecek et al. 2011). The per-base mutation rate was set to $1.28 \times 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 $\theta$ estimate and the fast approximate method (Kamm et al. 2016). The maximum value of $\rho$ 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 $\pm 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 $x_s$ and ending at $x_e$), exon length ($L_{exon} = x_e - x_s$), and gene length ($L_{tot}$). 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 $r_0$, 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 on the species, recombination can start to decrease before of just after gene boundaries. Accordingly, recombination rates in 5' ($r_5$) and 3' ($r_3$) do not necessarily correspond to hotspot rates. More precisely, $r_5 = r_{5,max}e^{-a d_5}$, where $r_{5,max}$ is the recombination rate at the center of the hotspot and $d_5$ is the distance between the hotspot center and the start of the gene (similar interpretation for $r_3$). Unfortunately, $r_5$ and $d_5$ (respectively, $r_3$ and $d_3$) are not identifiable. For each exon, the predicted recombination rate is thus given by

$$r_{\text{predicted}} = r_0 + \frac{R_1 + R_2}{L_{\text{exon}}} \quad (1)$$

with

$$R_1 = r_5 \frac{e^{-a x_s} - e^{-a x_e}}{a} \quad (2)$$

and

$$R_2 = r_3 \frac{e^{-a(L_{\text{tot}} - x_e)} - e^{-a(L_{\text{tot}} - x_s)}}{a} \quad (3)$$
The model was fitted with nonlinear least squares imple-
mented in the “nls” R function which returned predictions
and standard errors. The goodness of fit was assessed by plot-
ting 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 θ̂<sub>π</sub> in each exon rank and
gene classes in a similar manner as for recombination gra-
dients. 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 θ̂<sub>π</sub> per site with
"VCFtools" on the final VCF dataset used for "LDhat"
(Danecek et al. 2011). The average θ̂<sub>π</sub> per exon rank and
gene class was calculated as the sum of all θ̂<sub>π</sub> overlapping
a given exon rank divided by the total number of sites cov-
ered by these exons.

**Supplementary Material**

Supplementary material is available at *Molecular Biology and Evolution* online.

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 insightful discussions. Elise Rolland helped with the
pipeline. We used the GENOUEST computing facility for
bio-informatic analyses.

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.

Conflict of Interest

None declared.

Data Availability

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

**References**

1001 Genomes Consortium. Electronic address: magnus.nordborg@gmi.oegaw.ac.at and 1001 Genomes Consortium. 1,135 genomes reveal the global pattern of polymorphism in *Arabidopsis thaliana*. Cell. 2016;166(2):481–491. https://doi.org/10.1016/j.cell.2016.05.063.

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. https://doi.org/10.1126/science.1216872.

Auton A, McVean G. Recombination rate estimation in the presence of hotspots. Genome Res. 2007;17(8):1219–1227. https://doi.org/10.1101/gr.6386707.

Auton A, Myers S, McVean G. 2014. Identifying recombination hot-spots using population genetic data, arXiv, arXiv:1403.4264 [q-bio], preprint: not peer reviewed. https://doi.org/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. https://doi.org/10.1371/journal.pgen.1003984.

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. https://doi.org/10.7554/eLife.24133.

Barroso GV, Dutheil JY. The landscape of nucleotide diversity in *Drosophila melanogaster* is shaped by mutation rate variation. Peer Community J. 2023;3:e40. https://doi.org/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. https://doi.org/10.1126/science.1183439.

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. https://doi.org/10.1016/S0092-8674(02)01167-4.

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. https://doi.org/10.1073/pnas.94.15.8058.

Brazier T, Glémin S. Diversity and determinants of recombination landscapes in flowering plants. PLoS Genet. 2022;18(8):e1010141. https://doi.org/10.1371/journal.pgen.1010141.

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. https://doi.org/10.1038/nature11089.

Buard J, de Massy B. Playing hide and seek with mammalian meiotic crossover hotspots. Trends Genet. 2007;23(6):301–309. https://doi.org/10.1016/j.tig.2007.03.014.

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. https://doi.org/10.1038/s41467-021-27432-z.

Choi K, Henderson IR. Meiotic recombination hotspots - a comparative view. Plant J. 2015;83(1):52–61. https://doi.org/10.1111/tpj.15152.

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. https://doi.org/10.1038/ng.2766.

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. https://doi.org/10.1371/journal.pgen.1006799.
Cooper TJ, Garcia V, Neale MJ. Meiotic DSB patterning: a multifaceted process. Cell Cycle. 2016:15(1):13–21. https://doi.org/10.1080/15384101.2015.1093709.

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. https://doi.org/10.1093/bioinformatics/btr330.

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. https://doi.org/10.1093/molbev/msx272.

de Massy B. Initiation of meiotic recombination: how and where? conservation and specificities among eukaryotes. Annu Rev Genet. 2013:47(1):563–599. https://doi.org/10.1146/genet.2013.47 issue 1.

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. https://doi.org/10.1093/genetics/132.1.113.

Dluzewska J, Szymanska M, Ziolkowski PA. Where to cross over? defining crossover sites in plants. Front Genet. 2018:9:609. https://doi.org/10.3389/fgene.2018.00609.

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. https://doi.org/10.1105/tpc.9.9.1633.

Duret L. Why do genes have introns? Recombination might add a new piece to the puzzle. Trends Genet. 2001:17(4):172–175. https://doi.org/10.1016/S0168-9525(01)02236-3.

Dutheil JY. On the estimation of genome-average recombination rates. Genetics. 2024:227(2):iyae051. https://doi.org/10.1093/genetics/iyae051.

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. https://doi.org/10.1073/pnas.97.9.4718.

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. https://doi.org/10.1016/j.tig.2014.05.002.

Grey C, Baudat F, de Massy B. PRDM9, a driver of the genetic map. PLoS Genet. 2018:14(8):e1007479. https://doi.org/10.1371/journal.pgen.1007479.

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. https://doi.org/10.1038/s41588-019-0518-4.

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. https://doi.org/10.1111/mec.2018.27 issue 11.

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. https://doi.org/10.1073/pnas.1713225114.

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. https://doi.org/10.1073/pnas.1319032110.

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. https://doi.org/10.1126/science.adj7026.

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. https://doi.org/10.1534/genetics.104.036509.

Joseph J, Prentout D, Lavéré A, Tricou T, Duret L. High prevalence of Prdm9-independent recombination hotspots in placental mammals. Proc Natl Acad Sci USA. 2024:121(23):e2401973121. https://doi.org/10.1073/pnas.2401973121.

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. https://doi.org/10.1534/genetics.115.184820.

Kauppi L, Jeffreys AJ, Keeney S. Where the crossovers are: recombination distributions in mammals. Nat Rev Genet. 2004:5(6):413–424. https://doi.org/10.1038/nrg1346.

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. https://doi.org/10.1038/s41467-018-04562-5.

Lam I, Keeney S. Nonparadoxical evolutionary stability of the recombination initiation landscape in yeast. Science. 2015:350(6263):932–937. https://doi.org/10.1126/science.aad0814.

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. https://doi.org/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. https://doi.org/10.1371/journal.pcbi.1003118.

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. https://doi.org/10.1038/ncomms7648.

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. https://doi.org/10.1093/molbev/msac008.

Lloyd A. Crossover patterning in plants. Plant Reprod. 2022:36(1):55–72. https://doi.org/10.1007/s00497-022-00445-4.

Loewe L, Charlesworth B. Background selection in single genes may explain patterns of codon bias. Genetics. 2007:175(3):1381–1393. https://doi.org/10.1534/genetics.106.065557.

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. https://doi.org/10.1038/s41477-020-00834-5.

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. https://doi.org/10.1038/359154a0.

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. https://doi.org/10.1186/s13059-017-1326-8.

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. https://doi.org/10.1105/tpc.18.00750.

Mézard C. Meiotic recombination hotspots in plants. Biochem Soc Trans. 2006:34(4):531–534. https://doi.org/10.1042/BST0340531.

Myers S. A fine-scale map of recombination rates and hotspots across the human genome. Science. 2005:310(5746):321–324. https://doi.org/10.1126/science.1117196.

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. https://doi.org/10.1371/journal.pgen.1004234.

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. https://doi.org/10.1101/gr.067181.107.

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. https://doi.org/10.1016/j.cell.2011.02.009.

Raj A, Stephens M, Pritchard JK. fastSTRUCTURE: variational inference of population structure in large SNP data sets. Genetics. 2014;197(2):573–589. https://doi.org/10.1534/genetics.114.164350.

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. https://doi.org/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 salmornids, bioRxiv, preprint: not peer reviewed. https://doi.org/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. https://doi.org/10.1093/gbe/evv189.

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. https://doi.org/10.1534/genetics.119.302406.

Samuk K, Noor MAF. Gene flow biases population genetic inference of recombination rate. G3. 2022;12(11):jkac236. https://doi.org/10.1093/g3journal/jkac236.

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. https://doi.org/10.1093/genetics/126.4.813.

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. https://doi.org/10.1105/tpc.111.093674.

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. https://doi.org/10.1105/tpc.15.00391.

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. https://doi.org/10.1126/science.aad0843.

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. https://doi.org/10.1111/nph.1202.196 issue-3.

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. https://doi.org/10.1038/nature09869.

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. https://doi.org/10.1093/mp/sst123.

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. https://doi.org/10.1534/genetics.117.300502.

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. https://doi.org/10.1126/science.aab3761.

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. https://doi.org/10.1038/s41588-020-00723-9.

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. https://doi.org/10.1038/ng.3748.

Tock AJ, Henderson IR. Hotspots for initiation of meiotic recombination. Front Genet. 2018;9:521. https://doi.org/10.3389/fgene.2018.00521.

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. https://doi.org/10.1186/1471-2164-13-27.

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. https://doi.org/10.1038/s41586-018-0063-9.

Webster MT, Hurst LD. Direct and indirect consequences of meiotic recombination: implications for genome evolution. Trends Genet. 2012;28(3):101–109. https://doi.org/10.1016/j.tig.2011.11.002.

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. https://doi.org/10.7554/eLife.01426.

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. https://doi.org/10.1038/s41588-019-0546-0.

Yang C, Yan J, Jiang S, Li X, Min H, Wang X, Hao D. Resequencing 250 soybean accessions: new insights into genes associated with agro-nomic traits and genetic networks. Genom Proteom Bioinform. 2021;20(1):29–41. https://doi.org/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. https://doi.org/10.1101/gad.270876.115.

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. https://doi.org/10.1038/s41588-021-00895-y.

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. https://doi.org/10.1038/s41588-020-00722-w.