Stepwise expansion of recombination suppression on sex chromosomes and other supergenes through lower load advantage and deleterious mutation sheltering

Stepwise expansion of recombination suppression on sex chromosomes and other supergenes through lower load advantage and deleterious mutation sheltering | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search Confirmatory Results Stepwise expansion of recombination suppression on sex chromosomes and other supergenes through lower load advantage and deleterious mutation sheltering View ORCID Profile Paul Jay , Amandine Véber , Tatiana Giraud doi: https://doi.org/10.1101/2025.06.27.661902 Paul Jay 1 Laboratoire d’Ecologie Alpine, CNRS, Univ. Grenoble Alpes, Univ. Savoie Mont Blanc , 38000 Grenoble, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Paul Jay For correspondence: paul.yann.jay{at}gmail.com Amandine Véber 2 MAP5, CNRS, Université de Paris Cité , 75006 Paris, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Tatiana Giraud 3 Ecologie Société Evolution, IDEEV , Bâtiment 680, 12 route RD128, 91190 Gif-sur-Yvette, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract Many organisms possess sex chromosomes with non-recombining regions that have expanded progressively. Yet, the causes of this stepwise expansion remain poorly understood. Here, using mathematical modeling and stochastic simulations, we show that recombination suppression can expand simply due to the widespread presence of deleterious recessive mutations in genomes. We demonstrate that a significant proportion of new inversions are initially advantageous because they carry fewer mutations than average. However, these less-loaded inversions generally fail to fix on autosomes because, as their frequency increases, the recessive deleterious mutations they carry are more likely to occur in a homozygous state, leading to a selective disadvantage. In contrast, the permanent heterozygosity of Y-like sex chromosomes shelters sex-linked inversions from this disadvantage, facilitating their fixation and thereby the stepwise expansion of non-recombining regions. We show that this sheltering effect leads to fixation probabilities exceeding those expected under drift alone. Once recombination is suppressed, deleterious mutations accumulate on inversions, which could select for recombination restoration. However, we show that the accumulation of overlapping genomic rearrangements following recombination suppression can prevent its restoration. Our theoretical model proposes a simple and testable framework explaining evolutionary strata on sex and mating-type chromosomes, and other supergenes. Introduction Many organisms have sex chromosomes with large non-recombining regions that expand in a stepwise manner, although the underlying mechanisms remain poorly understood ( Jay et al., 2024 ; Ponnikas et al., 2018 ; Saunders & Muyle, 2024 ; Wright et al., 2016 ). It has long been considered that recombination suppression on sex chromosomes gradually expands because selection favors the linkage of sexually antagonistic loci to sex-determining genes ( Rice, 1987 ; Ruzicka et al., 2020 ; Wright et al., 2016 ), generating evolutionary strata of differentiation between sex chromosomes. However, there has been no compelling evidence that sexual antagonism is actually responsible for the expansion of recombination suppression in sex chromosome so far ( Beukeboom & Perrin, 2014 ; Ironside, 2010 ; Jay et al., 2024 ; Ponnikas et al., 2018 but see Wright et al., 2017 ) and recombination suppression has also been reported to gradually expand around many fungal mating-type loci and other super-genes despite the lack of sexual antagonism ( Bazzicalupo et al., 2019 ; Branco et al., 2017 ; Hartmann et al., 2021 ; Jay et al., 2021 , 2024 ; Wang et al., 2013 ; Yan et al., 2020 ). In addition, theoretical issues have been raised about the model of sexual antagonism driving the evolution of sex chromosomes ( Cavoto et al., 2018 ). Altogether, this suggests that other mechanisms can drive the stepwise extension of recombination suppression( Jay et al., 2024 ). In 2021, we developed a general model showing that recombination suppression can extend stepwise around sex-determining or mating-type loci only because of the presence of deleterious mutations in genomes ( Jay et al., 2022 ). Following its publication in PLOS Biology , this model generated significant interest but also faced attacks from several authors ( Charlesworth & Olito, 2024 ; Lenormand & Roze, 2024 ; Olito & Charlesworth, 2023 ), in particular on the controls used to assess the parameter range under which the sheltering mechanism can have a significant effect. After multiple rounds of discussion, the PLOS Biology editorial team decided to retract the paper, arguing that our theory may be correct but that evidence was lacking. The original manuscript, along with the retraction notice, is available in PLOS Biology website ( https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3001698 ), with further discussion on these matters in https://hal.science/hal-04763742 . Here, we present new analyses of our original derivations and simulations, along with new simulations incorporating additional controls, which further support our theory and address previous criticisms. The rationale behind our theory is illustrated in Figure 1 and consists of three steps that can explain the evolution and maintenance of recombination suppression on sex chromosomes. These steps are independent and can each contribute separately or in combination to the evolution of sex chromosomes. This theory is based on the observation that most mutations are deleterious and partially recessive, resulting in a genomic landscape with many deleterious recessive mutations segregating at low frequencies within populations ( Agrawal & Whitlock, 2011 , 2012 ; Eyre-Walker & Keightley, 2007 ). Download figure Open in new tab Figure 1: Schematic diagram of the three steps of the model. Step 1 . Within any population and for any genomic region, diploid individuals (here represented by two homologous chromosomes) and haplotypes (i.e., specific combinations of mutations) carry variable numbers of partially recessive, deleterious mutations. A substantial fraction (about 50%) of haplotypes has fewer deleterious mutations than average in the corresponding DNA fragment and should be favored by selection. Chromosomal inversions therefore have a significant chance of capturing beneficial haplotypes. This is the lower-load advantage. Step 2 . On autosomes, the increase in frequency of a beneficial inversion leads to a higher frequency of homozygotes for the inversion, which have a fitness disadvantage because they are homozygous for all the recessive deleterious mutations carried by the inversion. The mean fitness of the inversion in the population therefore decreases with increasing inversion frequency, preventing the inversion to reach a high frequency. Permanent heterozygosity at a Y-like sex-determining allele protects linked inversions from this homozygosity disadvantage. Hence, beneficial inversions (e.g. carrying fewer deleterious mutations than average) can continue increasing in frequency and become fixed in the population of Y-like chromosomes. This is the sheltering effect. Step 3 . Fixed inversions on Y-like chromosomes suppress recombination within the corresponding segment. As a result of selective interference and Muller’s ratchet-like processes, the non-recombining region accumulates new deleterious mutations. This degeneration should favor the restoration of recombination, for example through inversion reversion. However, the accumulation of additional chromosomal rearrangements during degeneration may prevent the restoration of recombination. The first step in our model corresponds to the selection of inversions carrying by chance a lower load than average in the genomic region ( Figure 1 ). As highlighted by Nei et al., 1967 , the stochastic distribution of deleterious mutations generates variation in fitness at any genomic region due to differences in mutation loads. Chromosomal inversions occurring at random in genomes will capture a specific haplotype of deleterious variants, that can have a higher or lower fitness than average, just by chance. Inversions that capture fewer deleterious mutations than the population average for a given region ( i . e ., have a lower mutation load) benefit from a relative fitness advantage and should therefore increase in frequency ( Figure 1 ). While the inversions themselves are neutral, their association with a favorable haplotype provides a selective benefit. As inversions suppress recombination when heterozygous, inverted advantageous haplotypes cannot recombine with less-advantageous haplotypes and therefore keep their advantage through meiosis, in contrast to recombining haplotypes. Note that we will use the term “inversions” for simplicity, but the process holds with any other mechanism suppressing recombination; in fact, young evolutionary strata can be collinear ( Branco et al., 2018 ; Sun et al., 2017 ; Vittorelli et al., 2023 ). The second step in our model corresponds to the sheltering of the recessive mutations in the inversion due to their linkage with a permanently heterozygous allele, allowing inversions to fix despite their recessive load. Indeed, when a given inversion increases in frequency, homozygotes for the inversion also become more frequent. Provided that the inversion carries at least one recessive deleterious mutation, homozygotes will experience a disadvantage due to the expression of this recessive load. Such selection against homozygotes should prevent this inversion from reaching high frequencies. In contrast, if a loaded inversion captures a permanently heterozygous allele, such as a Y-like sex-determining gene, its deleterious mutations are maintained in a heterozygous state, preventing the expression of its recessive load, and thereby facilitating its fixation on the Y-like chromosome. This protection against recessive load thanks to the linkage to a permanently heterozygous allele is hereafter called the sheltering effect. The successive fixation of additional inversions linked to this Y-fixed inversion by the same process should cause the non-recombining region to expand further, thereby leading to the formation of sex chromosomes with evolutionary strata. It is important to note that the lower-load advantage (step 1) and the sheltering effect (step 2) are two distinct selective effects. The lower load is an intrinsic fitness advantage, driving the increase in frequency of inversions with fewer deleterious mutations than average; once the frequency of the less-loaded inversion has reached an appreciable level, a second phase starts during which the sheltering effect protects inversions linked to permanently heterozygous alleles from exposing their load at a homozygous state, allowing them to reach fixation. The sheltering effect could act in combination with other types of intrinsic advantage of inversions, such as sexual antagonism or local adaptation, provided that the inversions harbour at least one deleterious, partially recessive mutations ( Jay et al., 2024 ). The third step of our model corresponds to the maintenance of recombination suppression despite further mutation accumulation in less-loaded inversions. Once recombination is suppressed, the inversions will indeed accumulate deleterious mutations by selective interference and Muller’s ratchet-like processes ( Bachtrog, 2013 ; Berdan et al., 2021a , 2023 ). The resulting degeneration may eventually select for recombination restoration, for instance via inversion reversion ( Lenormand & Roze, 2022 ). However, the accumulation of overlapping rearrangements in such a region—due to genetic drift or selective mechanisms ( Berdan et al., 2023 )—may hinder recombination restoration. Indeed, in such cases, reversion of the initial inversions would not restore collinearity and therefore recombination. Here, we analyse the different mechanisms acting in the three steps described above, and show that, in combination, they can lead to the evolution and maintenance of recombination suppression on sex chromosomes under a broad range of parameters. We also investigate whether the combination of the lower-load advantage and the sheltering effect can explain the fusions sometimes observed between sex chromosomes and autosomes ( Charlesworth & Wall, 1999 ; Ferchaud et al., 2022 ; Pennell et al., 2015 ; Yashiro et al., 2021 ) and we explore the impact of various parameters on the probability of recombination suppression evolution. The accumulation of deleterious mutations following recombination suppression has been extensively studied ( Bachtrog, 2013 ; Blaser et al., 2014 ; Grossen et al., 2012 ; Nei, 1970 ), but we investigate here the converse, i . e ., that deleterious mutations could be a cause, and not only a consequence, of recombination suppression. Results Method overview We explored this three-step model of stepwise recombination suppression around permanently heterozygous alleles with infinite population deterministic models and individual-based simulations. Very little is known about the dynamics of deleterious mutations in genomic regions with polymorphic inversions ( Berdan et al., 2021b , 2023 ). We therefore had to use simulations to explore realistic scenarios involving various levels of genetic drift. In both the infinite population model and individual-based simulations, we modeled diploid populations with only partially recessive deleterious mutations. Each mutation altered fitness by a factor of 1 + hs in heterozygotes and 1 + s in homozygotes with s<0 and effects combining multiplicatively across loci. We first considered that all mutations occurring in genomes had the same dominance ( h ) and selection ( s ) coefficients and then relaxed this hypothesis. As the effect studied depends on the number of deleterious mutations in the studied region, we present results as a function of U , the region-wide deleterious mutation rate, defined as U = n · μ , where n is the size of the region and μ the per-site mutation rate. For reference, the estimated genome-wide deleterious mutation rates in humans and Drosophila melanogaster are approximately 2.2 and 1.3 mutations/genome.generation, respectively ( Agrawal & Whitlock, 2012 ). Given the genome sizes of these species, a 10 Mb region corresponds to U ≈ 0.007 in humans and U ≈ 0.07 in D. melanogaster . Individuals were considered to have two pairs of chromosomes, one of which harbored a locus with at least one allele permanently or almost permanently heterozygous (see Methods for details). Several situations were considered, mimicking those encountered in XY sex-determination systems, in fungal mating-type systems or in overdominant supergenes. We simulated the evolution of recombination modifiers suppressing recombination across the fragment in which they reside ( i . e ., cis- modifiers), either exclusively in heterozygotes (mimicking for example chromosomal inversions), or in both heterozygotes and homozygotes ( e . g ., histone modifications). Each of these recombination modifiers appeared in a single haplotype (i.e. in a single chromosome of a single individual), and was thus in linkage disequilibrium with a specific set of mutations. We first considered recombination modifiers that were neutral by themselves, so that their fitness was exclusively dependent on the number of deleterious alleles within the captured segment. Therefore, the fate of these inversions was only influenced by their load (step 1) and possibly a sheltering effect (step 2). We then considered other scenarii, under which recombination modifiers were intrinsically beneficial or deleterious, as proposed for inversions for example ( Kirkpatrick, 2010 ; Villoutreix et al., 2021 ). We first compared the dynamics of inversion-mimicking mutations in an autosome to those capturing a male-determining allele in an XY system, males being XY and females XX, the male-determining allele being permanently heterozygous; we then considered other types of recombination modifiers and heterozygosity rules. Less-loaded inversions are frequent in genomes and are advantageous (STEP 1) As noted by several authors ( Connallon & Olito, 2021 ; Nei et al., 1967 ; Olito & Abbott, 2020 ), inversions capturing fewer deleterious mutations than the population average in the focal genomic segment should increase in frequency. In infinite populations, the number of mutations harbored by individuals within a genomic segment with n sites follows a binomial distribution of parameters n and q , with q the mean frequency of mutations at mutation-selection equilibrium (Figure S1). On average, individuals harbor nq mutations on each of their chromosomal segments of size n . Under realistic parameter values, the vast majority of large chromosomal regions therefore carry several deleterious mutations. For instance, considering s =-0.001, h =0.1 and U =0.001 ( e . g. μ =1×10 −09 and n =1Mb), more than 99.999% of chromosomal fragments carry a least one mutation, the mean number of mutations being nq = 10 (Figure S1). If m is the number of recessive deleterious mutations captured by a given inversion, the mean fitness of individuals homozygous for the inversion, heterozygous for the inversion or lacking the inversion, in an infinite population, can easily be expressed as a function of n, q, h, s , and m . [see Methods; ( 28, 37 )]. Once formed, inversions should initially increase in frequency if heterozygotes for the inversion are fitter than homozygotes without the inversion (W NI >W NN ), which is the case if the inversion carries fewer mutations than the population average [ m < nq; ( 28 )]. Repeated sampling from binomial distributions with a wide range of parameters showed that more than half the inversions occurring in genomes captured fewer deleterious mutations than the population average ( figure 2a ; hereafter referred to as “less-loaded inversions”). Indeed, the distribution of mutation number across individuals is almost symmetric when nq is high (as the binomial distribution converges to a normal distribution; Figure S1). However, when nq is low, the distribution is zero-inflated, increasing the probability of inversions being less-loaded. Therefore, a substantial fraction of inversions occurring in genomes are beneficial when they form ( i . e ., when rare enough to only occur as heterozygotes). For instance, considering U =0.001 ( e . g ., with n =1Mb and μ =1e-9), with h values ranging from 0 to 0.5 and s values ranging from -0.001 to -0.25, between 36% and 98% (mean=70%) of the 1 Mb inversions are beneficial, carrying fewer recessive deleterious mutations than average ( Figure 2a ). Download figure Open in new tab Figure 2: Less-loaded inversions are common and confer a fitness advantage. a . Probability that occurring inversions carry fewer mutations than average in an infinite population, as a function of the region-wide mutation rate (U), the selection coefficient (s), and the dominance coefficient (h) of segregating mutations. Mutations are assumed to be at mutation–selection equilibrium (see Methods). An inversion is considered “less-loaded” if the number m of mutations that it has captured is smaller than the average number of mutations, nq, obtained as the product of the size n of the inversion and of the equilibrium frequency of deleterious alleles q. As m is an integer, transitions between different values of m result in discontinuities in the plotted probabilities. b . Average initial selective advantage of less-loaded inversions in an infinite population. For each parameter combination, we averaged the fitness advantage of all inversions satisfying m < nq. The selective advantage of a given inversion was calculated by comparing the fitness of heterozygous carriers to the population mean fitness. On autosomes, the selective advantage of inversions should erode when the inversions increase in frequency because of the increasing frequency of individuals homozygous for the inversions, exposing the recessive load of these inversion. c . Relative number of mutations initially present in the inversions that eventually fixed on sex chromosomes in finite populations (N = 1000), normalized by the average number of mutations across all inversions. Values below 1 indicate that fixed inversions initially had a lower load than average. The scripts used to produce this figure are available on GitHub. The selective coefficient of less-loaded inversions depends on the number of deleterious mutations that they carry relative to the population average ( Figure 2b ). For instance, with U =0.001, h =0.1 and s =-0.001, inversions with five mutations occur with a probability of 0.037 and have a fitness advantage of 0.0005, whereas inversions with two mutations occur with probability 0.0022 and have a fitness advantage of 0.0008. Overall, under the range of parameter values studied, the average fitness advantage of less-loaded inversions upon formation is substantial in the parameter space studied, being for instance 0.0008 when U =0.001 and 0.0027 when U =0.005. Simulations in finite populations of different sizes confirmed that most inversions (mean=66% in the range of parameter values studied) had a fitness advantage upon formation, and therefore tend to increase in frequency. Notably, simulations showed that inversions that went to fixation on the Y sex chromosome tended to show fewer mutations than average ( Figure 2c ), confirming that the lower-load advantage could be a substantial driver of the increase in frequency of inversions. The simulations also showed that inversions could be favored if they captured mutations that were rarer than average (Figures S2-3). Among inversions that went to fixation, 36 % were initially mutation-free, but the remaining 64% carried deleterious mutations ( Figure 2c ). Loaded inversions are much more likely to fix when they capture a Y-like sex-determining locus (Step 1 + Step 2; Deterministic analyses) To determine whether the permanent heterozygosity of Y-like sex chromosome affects the fate of less-loaded inversion due to the sheltering effect (Step 2), we compared the evolutionary trajectory of inversions on autosomes versus those capturing the sex-determining locus on the Y chromosome. First, we considered deterministic trajectories, without drift, and initially assumed that inverted and non-inverted segments no longer accumulated deleterious mutations after their formation (this strong assumption being relaxed later). We found that the frequency of less-loaded inversions tended to remain low in autosomes, whereas many less-loaded inversions became fixed in the population of Y chromosomes ( Figure 3a ). This supports our hypothesis that the permanent heterozygosity of Y chromosome facilitates the fixation of loaded inversions by preventing recessive deleterious mutations to be expressed at the homozygous stage. To determine the conditions allowing autosomal and Y-linked inversions to be maintained polymorphic or fixed in their chromosomal population, we expressed the equilibrium frequency of less-loaded inversions ( Figure 3d ) as a function of m , the number of deleterious mutations carried by the inversion (see Methods and appendix). Assuming that q <<1 and s <<1, we found that inversions on autosomes became fixed when m < qhn/(1-h ), whereas they stabilized at intermediate frequency when qhn/(1-h) < m < nq (see Methods). Inversions on autosomes capturing more than qhn/(1-h) recessive deleterious mutations indeed suffered from a homozygote disadvantage preventing them from reaching high frequencies ( Figure 3a ). This was the case for most autosomal inversions under realistic parameter values ( e . g ., for more than 99.99 % of inversions when U =0.001, h =0.1 and s=-0 .001; Figure 3d ). Download figure Open in new tab Figure 3: Less-loaded inversions are more likely to fix on Y chromosomes than on autosomes. a . Deterministic change in frequency for 2-Mb inversions, either capturing the sex-determining allele on the Y chromosome or occurring on an autosome. For clarity, the frequency shown for Y-linked inversions refers to their frequency within the population of Y chromosomes. Inversions are assumed to carry 5% fewer mutations than the population average ( m = 0.95 × nq ), with mutations at mutation–selection equilibrium (μ = 10− 8 , U = 0.02, h = 0.1). Supplementary Figures S4–S8 illustrate additional cases, including inversions with partial linkage to the Y sex-determining allele, varying mutation loads, and inversions linked to loci with alternative heterozygosity patterns or located on the X chromosome. b . Change in inversion frequency in stochastic simulations of a population of 1,000 individuals, under various parameter-value combinations differing in the mutation rate, as well as in the dominance and selection coefficients of mutations. All mutations are assumed to have the same s and h values. The figure shows the frequency trajectories of 10,000 independent inversions located either on an autosome or a proto-Y chromosome. Each line represents a single inversion (i.e., an independent simulation). Results for inversions of different sizes, X-linked inversions and other parameter sets are shown in Supplementary Figures S10–S13. c . Change in the mean number of mutations carried by the inverted segments across the same 10,000 simulations as in panel b. Each line represents the trajectory of a single inversion. A close-up view of the dynamics of mutation accumulation is provided in Supplementary Figure S9. d . Proportion of inversions that are expected to fix in deterministic simulations without mutation accumulation after inversion formation—defined here as the fraction of inversions that confer a selective advantage over the course of the simulation. For Y-linked inversions, this corresponds to inversions carrying fewer mutations than the population average ( m < nq ). For autosomal inversions, this corresponds to inversion with m < qnh /(1− h ), where q is the equilibrium frequency of deleterious alleles. e . Proportion of inversions that reached fixation after 10,000 generations in stochastic simulations with N =1000 and N =10,000, across different combinations of parameter values. This panel extends the results of the panel b to additional values of the dominance coefficient ( h ), selection coefficient (s), and region-wide mutation rate ( U ). For each parameter combination, 10,000 inversions capturing a random genomic fragment were simulated. Results for other value of U I and s are available in supplementary Figure S14. In panels b, c and e, the dataset used is the same as the one used in Jay et al. (2022) . In contrast to Figure 3c in Jay et al., 2022 all inversions are shown here, and not only loaded inversions that survived the first 20 generations By contrast, permanently heterozygous alleles protect inversions from this homozygote disadvantage, allowing these inversions to fix. We found that, without drift and without mutation accumulation after inversion formation, all inversions capturing the Y-like sex-determining allele became fixed if they carried fewer mutations than average ( i . e. m < nq ). We therefore found that, within a realistic range of parameters, inversions were much more likely to become fixed if they captured the male-determining allele on the Y chromosome than if they were unlinked to this allele ( e . g ., on an autosome; Figure 3d and Figures S4-9). For instance, with U =0.001 ( e . g ., with n =1Mb and μ =1e-9), h =0.1, s=-0 .001, 47% of inversions occurring on the Y chromosome are expected to become fixed in the absence of genetic drift, versus only 0.0045% of inversions on autosomes. Genetic drift and mutation accumulation do not necessarily prevent Y-linked inversion fixation (Step 1 + Step 2 ; Simulations with drift) The dynamics of inversion frequency become more difficult to predict deterministically when accounting for the accumulation of deleterious mutations after the inversion arises. Indeed, while an approximation for infinite populations can be made ( Connallon & Olito, 2021 ; Nei et al., 1967 ; Olito & Abbott, 2025 ), it assumes that inverted segments evolve under mutation-selection dynamics with recombination. However, such assumptions do not hold in many situations. Indeed, low-frequency inversions in finite populations should almost never occur as homozygotes and should therefore evolve with almost no recombination. This is also true for inversions capturing a permanently heterozygous allele, which never undergoes recombination. We confirmed, by simulations, that the above approximations for an infinite population strongly departed from the situation observed in finite populations, in which the dynamics of mutation accumulation in inversions involve a mixture of a Muller’s ratchet-like regime (in sex chromosomes or when rare in autosomes, at the start of their spread) and a mutation-selection-drift regime with recombination ( e . g ., Figure S9). We therefore used individual-based simulations to study the fate of inversions in finite populations and accumulating deleterious mutations. Simulations of the inversion frequency trajectories in finite populations with mutations occurring confirmed the tendency identified in the deterministic model without mutation accumulation: over most of the parameter space explored, inversions were much more likely to fix if they captured the sex-determining allele on the Y chromosome than if they were located on autosomes ( Figure 3b,e and Figures S10-14). The probabilities of inversion fixation were higher on the Y chromosome than on autosomes across 94.9% of the parameter space explored with N=1000. On average, inversions were seven times more likely to fix on Y chromosomes than on autosomes with N =1000 and 341 times when N =10,000 ( Figures 3e and S14). Difference on inversion fixation rate between Y chromosome and autosomes were much pronounced when mutations segregating were relatively highly recessive (h<0.2), being 15 times higher with N=1000 and 2295 times higher when N =10,000. A non-negligible number of autosomal inversions carrying a mutation load segregated for hundreds of generations. For example, with N =1000, U =0.01, s= -0.01, h =0.1, 97 out of 10,000 inversions were still segregating after 500 generations. However, in this example, all these autosomal inversions were lost at the end of simulations ( i . e ., after 10,000 generations, Figure 3b ). These inversions initially increased in frequency because, by chance, they had a lower-than-average mutation load, but their homozygote disadvantage prevented them from reaching fixation. They were eventually lost because they accumulated further mutations relative to non-inverted segments ( Figures 3c and S9). Under certain conditions, in particular with low s and h values and high U values, autosomal inversions tended to remain at an intermediate frequency for a long time, reflecting their overdominant behavior ( Figure 3b and S10-13). Under the parameter range studied, some autosomal inversions nevertheless went to fixation, 72% of them being initially mutation-free ( Figure 3b,e and Figure S10-14). By contrast, substantial fractions of less-loaded inversions capturing the permanently heterozygous sex-determining allele on the Y chromosome reached fixation in the Y chromosome population; this was the case even for inversions that were not mutation-free ( Figures 3b,e and S10-14). For example, for N =1000, U =0.01, s=- 0.01, h =0.1, 55 out of the 10,000 Y-linked inversions were still segregating after 500 generations, and 30 of them became fixed in the Y chromosome population; in contrast, no inversion reached fixation on autosomes under the same conditions ( Figure 3b ). New mutations occurred on Y-linked inversions, but they did not accumulate rapidly enough to prevent these 30 inversions from increasing in frequency and reaching fixation ( Figure 3c and Figure S13). Only 36 % of the Y-linked inversions that fixed under the parameter range studied were initially mutation-free. Permanent heterozygosity effectively resulted in directional selection for the less-loaded inversions, and not only for those free of mutations ( Figure 3 and S10-14), leading to rapid fixation before the accumulation of too many new deleterious mutations. To check that the higher proportion of inversions fixing on the Y chromosome was not only due to the smaller effective population size of these chromosomes relative to autosomes, we compared the fixation rates of 100,000 inversions on Y chromosomes and autosomes with identical effective population sizes ( i . e ., autosomes with ¼ of the standard population size, Figure 4 ). This comparison was performed on a reduced parameter space to limit computation time. Our results showed that inversions were more likely to fix on the Y sex chromosome than on autosomes with identical population sizes across 75% of the parameter space explored ( Figure 4 ). On average, inversions were 18 times more likely to fix on Y chromosomes than on autosomes with the same effective population size ( Figure 4 ; note that, due to the limited parameter space explored, these numbers should not be considered as providing estimates of what occurs in nature). In contrast, the rate of inversion fixation on autosomes in populations with a ¾ reduction in size (N = 3,125) was only slightly higher than in autosomes in populations with the standard size (N = 12,500). These findings indicate that, while genetic drift plays a role in the fixation of inversions on Y-like sex chromosomes, the sheltering effect is the primary factor underlying the higher fixation rate of inversions on the Y sex chromosome compared to autosomes in our model. Furthermore, this comparison shows that the magnitude of the sheltering effect is large, allowing inversion fixation events on Y-like chromosomes to occur nearly 20 times more frequently than genetic drift alone. Download figure Open in new tab Figure 4: The sheltering effect is a major factor contributing to inversion fixation rate on Y-like sex chromosomes. Proportion of inversions that reached fixation after 10,000 generations in stochastic simulations, across different combinations of parameter values. This is the same analysis as the one performed in Figure 3e, but with different conditions. This figure reports new simulations, absent from Jay et al., (2022) . In this second set of simulations, mutations had their fitness coefficients sampled from a gamma distribution of mean -0.001 or -0.01, and their dominance coefficient randomly sampled with uniform probabilities among either [0.1, 0.2, 0.3, 0.4, 0.5] (mean=0.3, no fully recessive mutations), or [0.0, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5] (mean=0.22, with fully recessive mutations). In addition to the simulations of the inversion dynamics on Y chromosomes and on autosomes with identical population sizes (N=12,500), we also performed simulations with an autosomal population of size equivalent to the effective population size of the Y chromosome when the total population size is 12,500, i.e. 3,125 (¼ of 12,500). Finally, to show that the sheltering effect can act even in the absence of the lower-load advantage, we also performed simulations of inversions benefiting from an intrinsic selective advantage (sInv=0.01), in addition to the advantage or disadvantage conferred by their relative mutation load. A total of 100,000 inversions of 2Mb and 5Mb were studied for each combination of parameter values. Our simulations thus demonstrate that the sheltering effect significantly contributes to the fixation of loaded inversions on Y-like sex chromosomes by preventing the expression of captured partly recessive deleterious mutations in a homozygous state. Parameters impacting the probability of inversion spread and fixation (Step 1 + Step 2) Across all combinations of parameters analyzed, the difference in fixation rates between Y-linked and autosomal inversions was most pronounced when segregating mutations were substantially recessive ( i . e ., h < 0.2), assuming that all mutations shared the same dominance coefficient. In natural populations, however, dominance coefficients vary depending on the functional effects of mutations. While several studies have estimated that the average dominance coefficient of deleterious mutations is around h ≈ 0.25 ( Agrawal & Whitlock, 2012 ; Manna et al., 2011 , 2012), the full distribution of dominance effects around this mean remains poorly characterized. To study how the dominance coefficient distribution affects the probability of inversion spread and fixation, we conducted simulations using two distributions, both with an average close to h = 0.25, but one distribution allowing for fully recessive mutations, and the other excluding them ( Figure 4 ). In both scenarios, we observed that inversions were consistently more likely to fix on the Y chromosome than on autosomes. This suggests that the presence of a substantial proportion of relatively recessive mutations ( h < 0.2) is sufficient to drive the observed difference in fixation dynamics between Y-linked and autosomal inversions. As for any beneficial variant, even inversions benefiting from a selective advantage could be lost by drift, the probability of such loss depending on the number and type of mutations initially captured relative to the average for the population ( Figures 3 , 4 and S10-14). Compared to inversions occurring in mutation-poor regions, inversions occurring in regions in which large numbers of mutations segregate are more likely to capture many fewer mutations than average, and therefore to have a stronger relative advantage. In other words, inversions in mutation-dense regions ( i . e . with higher U ) have a wider fitness distribution. The probability of Y-linked inversion fixation thus increases with increasing inversion size and mutation rate ( Figures 3 , 4 , S15, S16 and S19). As expected ( Kimura & Ohta, 1969 ), we found that beneficial inversions took more generations to fix when population sizes were larger. This longer time favored the additional accumulation of deleterious mutations in inversions, decreasing their fitness and, in some cases, preventing their fixation. The probability of inversion fixation population therefore deceased with increasing population size both on autosomes and Y chromosomes ( Figures 3e and S14). Yet, this decrease was more pronounced on autosomes, resulting in a much greater difference in inversion fixation rates between Y chromosomes and autosomes when N = 10,000 (ratio 341:1 overall, 2295:1 when h<0.2) compared to N = 1,000 (ratio 7:1 overall, 15:1 when h<0.2; Figures 3e and S14). In addition, it is important to note that, depending on parameter values, the lower fixation rate of inversions in larger populations may be offset by the higher number of inversions expected to arise in such populations. The sheltering effect also contributes to fixation of inversions with other advantage than a lower load or with a small disadvantage To further highlight that the sheltering effect (step 2) and the lower-load advantage (step 1) correspond to distinct selective effects, we ran simulations with inversions that had an intrinsic selective advantage, e . g ., due to gene disruption or change in expression caused by the inversion breakpoints ( Dobzhansky, 1972 ; Kirkpatrick, 2010 ; Villoutreix et al., 2021 ) ( Figure 4 ). The initial fitness of these inversions depended on both their intrinsic advantage and the random load of deleterious mutations captured during their formation. Without the sheltering effect, we would expect a higher fixation rate of intrinsically beneficial inversions on autosomes than on sex chromosomes due to the more efficient selection allowed by their higher effective population size. In contrast, we found that intrinsically beneficial inversions were more likely to fix on the Y-like sex chromosome in 100% of the parameter space due to the sheltering of their recessive load. On average, such inversions were 18 times more likely to fix on the Y-like sex chromosome than on autosomes. Less-loaded inversions linked to a permanently heterozygous allele could also fix when assuming that inversion heterozygotes suffered from a fitness cost (Figures S18-19). However, with such a heterozygosity cost, only inversions with much fewer mutations than average were beneficial and could reach fixation, and not all less-loaded inversions (see Methods and appendix). Other systems with permanently heterozygous alleles, other recombination modifiers and sex chromosome-autosome fusion (Step 1 + Step 2) Further simulations showed that inversions were also more frequently fixed around other types of supergenes with permanently heterozygous alleles than elsewhere on autosomes, and in particular when: i) two or more permanently heterozygous alleles segregated at a mating incompatibility locus, modeling plant self-incompatibility or fungal mating-type systems [Figures S5 and S15; ( Hartmann et al., 2021 ; Takayama & Isogai, 2005 )]; ii) the alleles were not permanently heterozygous, but were strongly overdominant, thus occurring mostly in a heterozygous state, as for several supergenes controlling color polymorphism [Figure S8 ; ( Jay et al., 2021 ; Küpper et al., 2016 ; Wang et al., 2013 )]; iii) the fitness effects of mutations occurring across the genome were drawn from a gamma distribution ( i . e ., the mutations segregating in populations had different fitness effects from one another; Figure S16); iv) recombination modifiers suppressed recombination even when homozygous, as it is the case for histone modifications or methylation ( Boideau et al., 2021 ), rather than solely when heterozygous, as in the case of inversions (Figure S15); and v) the inversion was in strong but incomplete linkage with the permanently heterozygous allele ( e . g . 0.1 cM away from the allele, Figure S4-8). Our model can thus account for the existence of inversions very close, but not fully linked to a mating-type locus, as reported in the chestnut blast fungus ( Hartmann et al., 2025 ; Kubisiak & Milgroom, 2006 ; Stauber et al., 2021 ) In addition, we found that the sheltering of deleterious recessive mutations could also lead to the fusion of the permanently heterozygous sex chromosome with an autosome, when simulating chromosome fusions associated with an extension of the non-recombining region to the newly fused autosome (Figure S17). Less-loaded inversions linked to a permanently heterozygous allele could also increase in frequency and fix in simulations assuming a haplodiplontic life cycle ( i . e ., an alternation of haploid and diploid phases; Figures S20-21). However, the occurrence of a haploid phase enhanced the purge of recessive deleterious mutations, leading to a lower population mutation load and thereby to a narrower fitness distribution of inversions upon formation [Figure S21 ; ( Kondrashov & Crow, 1991 ; Scott & Rescan, 2017 )]. A lower proportion of inversions were therefore fixed in haplodiplontic populations than in purely diploid populations (Figure S20). Evolution of non-recombining sex chromosomes with evolutionary strata despite possible reversions (step 1 + step 2 + step 3) We have shown that, thanks to a combination of lower load (step 1) and sheltering (step 2), inversions capturing a permanently heterozygous allele can fix on Y-like chromosomes across a broad parameter range, with fixation probabilities exceeding those expected under drift alone. This process should thus occur repeatedly around permanently heterozygous alleles, leading to the formation of non-recombining sex chromosomes with a typical pattern of evolutionary strata. We investigated the formation of evolutionary strata through the successive fixation of multiple, potentially overlapping inversions, by simulating the evolution of large chromosomes experiencing recurrent inversion events, using parameter values representative of those observed in mammals ( Figures 5 , 6 and S22). A key focus was the potential for recombination restoration via inversion reversion, which may be favored when a fixed inversion on the Y chromosome accumulates a sufficiently high mutational load to exceed the population average ( Lenormand & Roze, 2022 ). Download figure Open in new tab Figure 5: Accumulation of successive inversions around a Y-like sex-determining allele in an XY system, leading to the formation of non-recombining sex chromosomes with evolutionary strata. Results of a simulation of N=1000 individuals, each with two pairs of 100 Mb chromosomes, over 100,000 generations. Chromosome 1 harbors an X/Y sex-determining locus at 50 Mb (individuals are XX or XY). In each generation, one inversion appears, on average, in the whole population, in an individual sampled uniformly at random, with the two recombination breakpoints sampled uniformly at random from k=100 potential breakpoints. a . Overview of chromosomal inversion frequency and position for 10 different generations. The width of the box represents the position of the inversion and the height of the box indicates inversion frequency. Inversions appearing on the Y chromosome are depicted in yellow, those appearing on the X chromosomes are depicted in gray. The colors are not entirely opaque, so that regions with overlapping inversions appear darker. Previously fixed inversions may be lost due to the occurrence of beneficial reversions and selection. b . Changes in the relative rate of recombination over the entire course of the simulation. The numbers of recombination events occurring at each position (binned in 1 Mb windows) are recorded at the formation of each offspring, across all homologous chromosomes in the population. Only recombination events between the X and the Y chromosomes are shown for chromosome 1 (i.e., recombination events between the two X chromosomes in females are not shown). Unlike chromosome 1, chromosome 2 harbors no permanently heterozygous alleles. All inversions on this chromosome suffer from homozygote disadvantage and very few inversions therefore become fixed on chromosome 2. See Figure S22 for a simulation with N=10,000 individuals. The datasets and scripts used to produce this figure are available on Figshare (doi:10.6084/m9.figshare.19704457) and GitHub. Download figure Open in new tab Figure 6: Effect of the number of potential inversion breakpoints on the evolution of recombination suppression in sex chromosomes. Each dot represents the result of a simulation with N=1000 individuals. For each number of breakpoints, 10 simulations were conducted. See Figure 5 for an example of such a simulation. a . Fraction of the length of the Y sex chromosome not recombining after 100,000 generations. b . Number of reversions occurring over the course of the 100,000 generations. Boxplot elements: central line: median, box limits: 25th and 75th percentiles, whiskers: 1.5x interquartile range. The dataset and script used to produce this figure are available on Figshare (doi:10.6084/m9.figshare.19704457) and GitHub. We simulated, over 100,000 generations, populations of N =10,000 or N =1,000 individuals, carrying two 100 Mb chromosomes, one of which harbored a mammalian-type sex-determining locus (XY males and XX females). Individuals experienced only deleterious or weakly deleterious mutations, with mutation rates and fitness effects similar to those observed in humans (fitness effect being drawn from a gamma distribution). The dominance coefficient of each mutation was chosen uniformly at random from a wide set of values (see Methods for details). At the start of each simulation, we randomly sampled k genomic positions that could be used as inversion breakpoints, with k being 10, 100, 1000 or 10,000. These genomic positions represent inversion hotspots, such as those that can be generated by repeats in genomes ( Porubsky et al., 2022 ). In each generation, we introduced j inversions, j being sampled from a Poisson distribution. To limit simulation times, we used high inversion rates, with one inversion occurring on average each generation in the whole population. The two breakpoints of each inversion were chosen at random from the k positions. It was, therefore, possible for two independent inversions to appear at the same position, allowing, in particular, the reversion of an inversion to its ancestral orientation, thereby restoring recombination. We assumed that inversions that partially overlapped another inversion or that captured a smaller inversion could not be reversed, i . e ., that recombination could not be restored in such situations even if subsequent inversions re-used the same breakpoints. Indeed, reversions of partially overlapped inversions do not restore ancestral arrangements but instead result in complex reshufflings of gene order and orientation and are therefore unlikely to restore recombination. We ran 10 simulations for each set of parameters. In all simulations assuming relatively large numbers of potential inversion breakpoints ( k =100, 1000, 10,000), the Y chromosome progressively stopped recombining with the X chromosome as it accumulated successive inversions fully linked to the male sex-determining allele ( Figures 5 , 6 and S22). After the occurrence of an initial inversion capturing the Y sex-determining allele, multiple inversions partially overlapping this inversion or other Y-fixed inversions were selected for, thereby generating a growing chaos of overlapping chromosomal rearrangements. The non-recombining region, thus, extended around the sex-determining locus in a stepwise manner, perfectly reflecting the evolution of sex chromosomes and other supergenes with evolutionary strata ( Figure 5 ). Some events in the gradual extension of recombination suppression were reversed, due to the occasional occurrence of beneficial reversions ( Figures 5 and 6 ). The accumulation of overlapping inversions was, however, more rapid than the occurrence of beneficial reversions, leading to a progressive extension of the non-recombining region ( Figures 5 and 6 ). By contrast, when we assumed a smaller number of potential breakpoints ( k = 10), recombination suppression between sex chromosomes evolved in only two of the 10 simulations, and over only a small genomic region ( Figure 6 ), owing to the more frequent occurrence of beneficial reversions. Discussion Evolution of recombination suppression around permanently heterozygous alleles through the combination of a lower-load advantage and the sheltering effect Our results show that recombination suppression on sex chromosomes and other supergenes can evolve simply because genomes harbor many partially recessive, deleterious variants. Our model for the evolution of sex chromosomes, and supergenes in general, is based on simple and widespread phenomena: i) inversions (or any recombination suppressor) can be favored solely because they contain fewer deleterious mutations than the population average, a situation applying to a substantial fraction of the inversions formed; ii) such inversions tend to display overdominance: they are beneficial in the heterozygous state but suffer from a homozygote disadvantage, which prevents them from reaching high frequencies and becoming fixed on autosomes; iii) when, by chance, such less-loaded inversions capture a permanently heterozygous allele, they do not suffer from this homozygote disadvantage and are therefore able to increase in frequency until they are fully associated with the permanently heterozygous allele ( e . g . they become fixed in the Y chromosome population). These three phenomena have been reported independently in several studies, but, to our knowledge, never in interaction [see ( Lenormand & Roze, 2022 ; Nei et al., 1967 ; Olito & Abbott, 2025 ) for i), ( Kirkpatrick, 2010 ; Ohta, 1971 ) for ii), and ( Antonovics & Abrams, 2004 ; Hartmann et al., 2021 ) for iii)]. The combined influence of the mechanisms related to ii) and iii) has been shown to promote sex chromosome-autosome fusion in highly inbred populations ( Charlesworth & Wall, 1999 ). We show here that such mechanisms can readily lead to the stepwise extension of the non-recombining region on sex chromosomes themselves, without the need for inbreeding. Moreover, unlike previous studies [ e . g . ( Charlesworth et al., 1987 ; Connallon et al., 2018 ; Olito & Abbott, 2025 )], we show that the higher probability of inversion fixation on Y chromosomes is not restricted to mutation-free inversions, but applies to any inversion loaded with deleterious recessive mutations. The proposed combined effects of lower load and deleterious mutation sheltering can also explain the fusions sometimes observed between sex chromosomes and autosomes, with the recombination suppression extending to a part of the fused autosome ( Charlesworth & Wall, 1999 ; Ferchaud et al., 2022 ; Pennell et al., 2015 ; Yashiro et al., 2021 ). The theory proposed here to explain the stepwise evolution of recombination suppression applies to any locus with at least one permanently or nearly permanently heterozygous allele. It only requires that, within a population, individuals carry different numbers of partially recessive deleterious mutations in a genomic region that can be subjected to recombination suppression. As discussed below, this situation is probably frequent in diploid, dikaryotic or heterokaryotic organisms. ( Berdan et al., 2023 ; Giner-Delgado et al., 2019 ; Porubsky et al., 2022 ; Wellenreuther & Bernatchez, 2018 ) It is important to emphasize that the lower-load advantage and the sheltering effect are two distinct and independent mechanisms that can influence the evolutionary fate of inversions. The lower load corresponds to a selective advantage, arising when an inversion captures by chance a genomic segment carrying fewer deleterious mutations than the population average in this genomic region. The fitness benefit conferred by such an inversion increases with the absolute difference in mutation load between the captured segment and the population mean. Consequently, larger inversions—by capturing more sites—have a wider fitness distribution, and can thus experience a stronger selective advantage. The sheltering effect corresponds to the protection from the homozygote disadvantage of inversion loaded with partially recessive deleterious mutations when they become frequent. Indeed, autosomal inversions carrying partially recessive deleterious mutations experience a reduction in fitness as their frequency increases, as they increasingly occur at the homozygous state. In contrast, inversions linked to the permanently heterozygous Y chromosome are protected from this homozygous disadvantage, allowing even highly loaded inversions to spread and potentially fix. As shown in our analyses, this difference in fate between autosomal and Y-linked inversions due to the sheltering effect does not only apply to less-loaded inversions. It applies to any inversion—whether its spread is driven by genetic drift or a selective advantage—as long as it carries recessive deleterious mutations. Consequently, this effect is expected to facilitate the spread of virtually any inversions of substantial size on Y-like chromosomes. The right controls to study the lower-load and sheltering effects In this manuscript, as well as in Jay et al. (2022) , we did not conduct any comparison with a fully neutral situation, i . e ., without any deleterious mutations in genomes (corresponding to s = 0), a situation under which the fixation of inversions would be influenced solely by genetic drift. Indeed, such a comparison, which was central in the debate surrounding the Jay et al. 2022 study ( Charlesworth & Olito, 2024 ; Olito & Charlesworth, 2023 ) is not relevant. The sheltering effect and the lower-load advantage only arise because the presence of deleterious mutations changes the fitness landscape. The observation that inversions might not be more likely to fix on Y chromosomes for certain parameter values in the presence of deleterious mutations compared to a scenario without any deleterious mutation does not negate the role of the sheltering effect and the lower-load advantage in promoting the fixation of inversions when deleterious mutations are segregating , in contrast to previous claims ( Charlesworth & Olito, 2024 ; Olito & Charlesworth, 2023 ). It is not relevant to ask whether the sheltering effect and the lower-load advantage lead to higher fixation rates of inversions on sex chromosomes in a scenario with deleterious mutations compared to a scenario without any deleterious mutations in genomes, unless we also take into account how frequent is this situation with only “purely neutral” mutations genomes in nature compared to genomes loaded with recessive deleterious mutations. Indeed, the probability of inversion fixation are conditional to the assumed scenario. One could argue that the scenario s =0 is not a realistic scenario but a control for the genetic drift effect, removing any selective effect. However, the sheltering effect is not an intrinsic benefit that can be directly compared to a case without any selective effect, as classically done in population genetics. The sheltering effect only mitigates the selective disadvantage caused by deleterious mutations; removing those mutations eliminates both the disadvantage and the sheltering effect itself, making the “ s =0” scenario an inappropriate control for assessing the impact of sheltering as compared to genetic drift (see also Box 3 in Jay et al., 2024 , and discussion in Saunders & Muyle, 2024 ). Indeed, deleterious mutations prevent the fixation of most inversions on autosomes (an effect that is absent when s=0), while the sheltering effect protects inversions from this fitness loss on Y-like chromosomes. The sheltering effect thus saves inversions specifically on Y-like chromosomes when deleterious mutations are segregating. Because deleterious mutations are known to segregate in natural populations ( Agrawal & Whitlock, 2012 ; Eyre-Walker & Keightley, 2007 ), these effects need to be considered if we are to understand the evolution of recombination suppression on sex chromosome in nature. As an analogy, consider a scenario where one aims to evaluate the effect of a drug (representing the sheltering effect) on the likelihood of developing a disease (analogous to the fitness cost of homozygosity), which is caused by pathogenic bacteria (representing recessive deleterious mutations). Comparing the survival probability of individuals who take the drug (Y-like chromosomes) in the presence of bacteria (s<0) and in the absence of bacteria (s=0) provides no information about the drug efficacy against the bacteria. Observing that, in the absence of bacteria , individuals who take the drug have a lower survival rate than those who do not tells us nothing about how effective the drug is when the bacteria are present . To assess the drug benefit, one must compare survival between individuals taking the drug (Y-like chromosomes) and those who do not take the drug (autosomes), bacteria being present (s<0). The protective effect of the drug (the sheltering effect) can indeed only be evaluated in the context where the underlying cause of the disease (here, the deleterious mutations) is present. It is however true that Y-like chromosomes have a reduced effective population size, which can also increase inversion fixation probabilities by inducing higher genetic drift. The appropriate control for genetic drift in our model is however not the case s =0, as explained above, but an autosome with the same effective population size as the Y sex chromosome, as it retains the effect of deleterious mutations while incurring similar genetic drift. This allows us to assess whether the observed higher fixation probability on the Y chromosome is only due to genetic drift, all other conditions being equal. Our new simulations using this control show that genetic drift only plays a minor role in explaining the higher inversion fixation rate on Y-like sex chromosomes than on autosomes, the major role being played by the sheltering effect. Conditions for the occurrence of the lower-load advantage and the sheltering effect In this study, we explored a broad parameter space, with N ranging from 1,000 to 12,500; s from -0.001 to -0.5; h from 0 to 0.5; μ from 1e-09 to 1e-08; and inversion sizes from 500 kb to 5 Mb (or region-wide mutation rate U from 0.0005 to 0.05). Our results demonstrate that the number of segregating deleterious mutations and their dominance coefficients are crucial factors for the occurrence of the lower-load advantage and the sheltering effect, with more mutations and more recessive mutations favoring inversion fixation. These quantitative effects are much more informative than an all-or-nothing comparison with no deleterious mutations at all and further show the potential crucial role of deleterious mutations in the evolution of stepwise recombination suppression. Numerous studies in nature have shown that a substantial proportion of new and segregating mutations are deleterious ( Agrawal & Whitlock, 2012 ; Eyre-Walker & Keightley, 2007 ). While the precise numbers and effects of these mutations are debated, it is widely accepted that genomes carry tens of thousands of harmful mutations. Therefore, any large inversion ( e . g ., 1 Mb) is expected to harbor multiple deleterious mutations. For instance, the average human genome contains approximately 4.1 to 5.0 million polymorphic sites, with an estimated 25% of these mutations being deleterious ( Racimo & Schraiber, 2014 ). Even if this estimate was overestimated by a factor of 100, megabase-scale inversions would still contain multiple deleterious variants, setting the stage for the lower-load advantage and the sheltering effect. Importantly, megabase-scale inversions, as we modelled, are commonly observed in natural populations. For instance, the two most recent evolutionary strata on the human Y chromosome span approximately 1 and 4 Mb, respectively ( Zhou et al., 2023 ), which corresponds to U =0.0007 and U =0.0028 when considering a genome-wide deleterious mutation rate of 2.2. Although empirical estimates of the dominance coefficient for mutations in natural populations are limited, studies in Drosophila , yeasts and nematodes have estimated the average dominance coefficient for deleterious mutations at approximately h = 0.25 ( Agrawal & Whitlock, 2012 ; Manna et al., 2011 ). This suggests that many deleterious mutations have dominance coefficients below 0.25. For example, yeast gene knockout data led to estimates of average dominance coefficient of 0.046 for mutations affecting catalytic functions ( Agrawal & Whitlock, 2011 ). In addition, there is substantial evidence that the distribution of dominance coefficients for deleterious mutations is right-skewed, with a notable overrepresentation of strongly recessive mutations ( Mrnjavac et al., 2025 ). These estimates indicate that the conditions for the sheltering effect are commonly met in natural populations, as we have shown that the effect is more pronounced with mutations having dominance coefficients below 0.2. Finally, it is important to note that evolutionary strata rarely evolve in natural populations. For example, the mammalian Y sex chromosome only experienced five successive events of recombination suppression across 250 million years ( Cortez et al., 2014 ; Zhou et al., 2023 ). Therefore, only a few lucky inversions with the right combinations of mutation number, selective and dominance coefficients are necessary to explain natural patterns of stepwise evolution of recombination suppression. Mutation accumulation does not necessarily prevent recombination suppression extensions On autosomes, inversions maintained at low frequencies due to their homozygous disadvantage are typically rapidly lost, as they continue to accumulate additional deleterious mutations over time. On the Y chromosome, however, we show that the rate of mutation accumulation following the formation of an inversion can be slow enough that it does not fully offset the initial selective advantage of less-loaded inversions. This can allow Y-linked inversions to reach fixation, despite ongoing degeneration. However, due to computational limitations, our analyses were restricted to relatively small population sizes. In larger populations, the time to fixation may be considerably longer, potentially allowing more deleterious mutations to accumulate before fixation occurs—thereby reducing the likelihood that such inversions will ultimately fix. To understand the evolutionary dynamics of inversion fixation in large populations more precisely, in particular under the influence of mutation accumulation, more complex models are needed that would capture the dynamics of deleterious mutations in spreading inversion, incorporating both Muller’s rachet and mutation-selection-drift dynamics. Long-term persistence of recombination suppression on sex chromosomes Our result show that deleterious mutations accumulate following Y-linked inversion fixation, as observed on many non-recombining sex chromosomes ( Bachtrog, 2013 ). This may lead to selection for inversion reversion, thereby restoring recombination ( Lenormand & Roze, 2022 ). The question as to whether this actually occurs is common to all mechanisms explaining evolutionary strata, including sexual antagonism. In this paper we showed that overlapping genomic rearrangements accumulating following recombination suppression could prevent recombination to be restored, depending on the number of putative inversion breakpoints. Indeed, partially overlapping inversions result in a complex reshuffling of gene order and orientation, preventing the restoration of recombination even in situations in which reversion could be selected for. When the number of potential breakpoints is relatively high, the probability of partially overlapping inversions occurring is higher than the probability of a reversion occurring, regardless of the relative rates of inversions and deleterious mutations. The reversion of inversions has in fact been reported only in rare cases in which inversions occur at specific breakpoints rich in repeated elements ( Cui et al., 2012 ; Hanson et al., 2014 ), as in our simulations when assuming only a few possible breakpoints in the genome. The number of genomic positions at which inversions can occur in natural conditions is unknown, but this number is likely to be high given the chaos of rearrangements observed in some sex and mating-type chromosomes with hundreds of different breakpoints ( Badouin et al., 2015 ; Branco et al., 2017 ; Carey et al., 2021 ; Porubsky et al., 2022 ). A recent study reported the recurrent appearance of inversions at the same positions, but also the existence of numerous potential breakpoints of inversions in human genomes ( Porubsky et al., 2022 ). Moreover, several studies have shown that chromosomal rearrangements rapidly accumulate in recently established regions of non-recombination ( Badouin et al., 2015 ; Carpentier et al., 2022 ) and overlapping inversions are observed on many sex chromosomes, which should prevent the restoration of recombination by reversions ( Bellott et al., 2014 ; Branco et al., 2017 ; Carey et al., 2021 ; Lemaitre et al., 2009 ; Skinner et al., 2021 ). Therefore, over a wide range of realistic parameter values, the reversion of inversions should not prevent the stepwise formation of non-recombining sex chromosomes. Lenormand & Roze, 2024 argued that, on a longer term and without dosage compensation, the fitness of individuals carrying the inversions having accumulated a further load should decrease to a point where species could go extinct, and therefore that the combined effect of the lower-load advantage and sheltering effect could not explain the evolution of recombination suppression on sex chromosomes. However, that recombination suppression on sex chromosomes can lead to species extinction or sex chromosome turn-over unless dosage compensation or other mitigation mechanisms eventually evolve does not in itself negates the potential role of the sheltering effect and overlapping rearrangement in allowing the evolution and short-term maintenance of recombination suppression (see also discussion in Saunders & Muyle, 2024 ). Conversely, our theory does not negate a possible role of dosage compensation for the long-term maintenance of sex chromosomes, but it shows that early dosage compensation or sexual antagonism may not be required to explain sex chromosome evolution, in contrast to the conclusions reached in previous studies ( Lenormand & Roze, 2022 , 2024 ). As in the case of asexuality for example, a selection on the short term does not mean that it will not lead to species extinction on the long term if mitigation mechanism does not evolve ( de Vienne et al., 2013 ). Furthermore, the early dosage compensation hypothesis requires a form of sexual antagonism, as it assumes the existence of numerous sex-specific gene regulators for genes on proto-Y sex chromosome, even before the evolution of recombination suppression ( Jay et al., 2024 ; Lenormand & Roze, 2022 ); this early dosage compensation mechanism ( Lenormand & Roze, 2022 ) cannot, therefore, explain evolutionary strata on fungal mating-type chromosomes, with which no form of antagonistic selection is associated ( Bazzicalupo et al., 2019 ; Hartmann et al., 2021 ). Parameters restricting the extension of recombination suppression Our model shows that stepwise recombination suppression can also evolve around supergenes carrying more than two permanently heterozygous alleles, such as plant self-incompatibility or mushroom (Agaricomycotina) mating-type loci. Such multi-allelic loci are however not known to often display extensions of recombination suppression beyond incompatibility loci ( Le Veve et al., 2022 ). This may be because the existence of multiple alleles allows the loss of degenerated alleles over the longer term: if a permanently heterozygous allele evolves a large non-recombining region, for instance because of an inversion fixation, and then degenerates through Muller’s ratchet-like processes, it can be lost by selection, in contrast to alleles in bi-allelic compatibility systems. In fungi, only biallelic mating-type loci have been reported so far to experience recombination suppression around their mating-type locus so far ( Jay et al., 2024 ). In addition, recessive self-incompatibility alleles can be homozygous in sporophytic systems ( Llaurens et al., 2009 ). In multi-allelic systems, we therefore expect to observe extensions of recombination cessation only around permanently heterozygous alleles, i . e ., dominant self-incompatibility alleles, and only rare, recent and non-degenerated inversions; long-read polymorphism sequencing data would allow testing this hypothesis. Assuming a reduced fertility in heterozygotes for inversions due to segregation issues possibly arising during meiosis ( Kirkpatrick, 2010 ) decreased the probability of inversion fixation, as expected. We found that large inversions were more likely to fix than small inversions, as they were more likely to capture haplotypes with much fewer mutations than average, thereby compensating the heterozygous cost; larger inversions may however be also more likely to induce segregation issues, although little data is available to date. Given the lack of knowledge about inversion rates in natural conditions and the computation challenge represented by the simulation of complex patterns of recombination, we used reasonably high inversion rates in our sex-chromosome evolution simulations, making it possible to observe the stepwise extension of a non-recombining region within 100,000 generations. Recent studies suggested that inversions may not be too rare and that inversion breakpoints may be widely distributed throughout the genome ( Hämälä et al., 2021 ; Porubsky et al., 2022 ; Todesco et al., 2020 ; Zhou et al., 2019 ). The use of different inversion rates might result in much shorter or much longer times for the stepwise extension of the non-recombining region, but should not change the final outcome. Of course, lower inversion rates and higher costs of inversions in terms of fertility in heterozygotes should decrease the rate of inversion fixation. In nature, the stepwise extension of the non-recombining region between sex or mating-type chromosomes often occurs over time scales of the order of tens of millions of generations ( e . g ., about 250 M years in human; ( Zhou et al., 2023 )), suggesting that the fixation of inversions are relatively rare events. Predictions regarding the variability of sex-chromosome structure across species We show that inversions are more likely to spread and fix in regions harboring many segregating deleterious mutations, in which they have a greater chance of capturing highly advantageous haplotypes. Our model therefore predicts that species harboring a large number of deleterious recessive variants, due to their small population size, short haploid phase, outcrossing mating system, high mutation rate and high levels of mutation recessiveness, for example, will be more prone to the evolution of large non-recombining regions with evolutionary strata on sex chromosomes than species with a low mutation load. In addition, in species with large population sizes, the time required for an inversion to become fixed may exceed the time during which the number of deleterious mutations in the inversion remains below average. This should prevent some inversions from becoming fixed, potentially decreasing the expansion rate of non-recombining regions on sex chromosomes in species with large population sizes. Depending on mutation effects and on the relative rates of inversions and mutations, this could however be compensated by the occurrence of a higher number of inversions each generation in large populations. Consequently, variations in population size, mutation rate, length of haploid phase and mating system (outcrossing versus selfing, or inbreeding) across lineages may account for the large variation in sex-chromosome structures in nature, with some organisms maintaining homomorphic sex chromosomes and others evolving highly differentiated sex chromosomes with multiple evolutionary strata ( Abbott et al., 2017 ). The more efficient purging of recessive deleterious mutations in species with an extended haploid phase ( Kondrashov & Crow, 1991 ; Scott & Rescan, 2017 ) could therefore potentially account for the smaller non-recombining regions observed on the sex chromosomes of plants and algae than in animals ( Coelho et al., 2018 ; Filatov, 2015 ). In fungi for example, multiple species with a diploid-like life cycle have repeatedly and independently evolved stepwise recombination suppression around mating-type loci, but not their closely related species with haploid-like life cycles ( De Filippo et al., 2025 ; Hartmann et al., 2021 ; Jay et al., 2024 ). A test of our model could thus look for association between the size of non-recombining regions or the number of evolutionary strata on sex chromosomes and estimates of the number of deleterious mutations segregating in genomes, or parameters predicting the accumulation of deleterious mutations, such as population size, mating systems, length of the haploid phase and mutation rates ( Jay et al., 2024 ). Conclusion Given its simplicity and wide scope of application, our model of sex chromosome evolution is a powerful alternative to other explanations ( Ponnikas et al., 2018 ), although the various theories are not mutually exclusive. The strength of our model lies in the absence of strong assumptions, such as sexually antagonistic selection or small population sizes. Our model, based on the often-overlooked observation that recessive deleterious mutations are widespread in genomes within natural populations, can explain the evolution of stepwise recombination suppression over a wide range of realistic parameter values. Furthermore, it can explain why some supergenes, such as fungal mating-type chromosomes and autosomal supergenes in butterflies and ants, display evolutionary strata ( Carpentier et al., 2022 ; De-Kayne et al., 2025 ; Jay et al., 2021 ; Yan et al., 2020 ), and why many mating-type loci in fungi display a stepwise extension of non-recombining regions despite the absence of antagonistic selection ( Branco et al., 2017 ; Hartmann et al., 2021 ). Our model can also explain why meiotic drivers, which are often permanently heterozygous, are often associated with large non-recombining region involving polymorphic chromosomal inversions ( Dyer et al., 2007 ; Reinhardt et al., 2014 ). Our model therefore provides a general and simple framework for understanding the evolution of non-recombining regions around many kinds of loci carrying permanently heterozygous alleles. Materials and Methods Infinite population deterministic model We consider the discrete-time evolution of an infinite size, randomly mating population experiencing only deleterious recessive mutations, with heterozygotes and homozygotes suffering from a 1-hs and 1-s reduction in fitness, respectively. At all n sites, mutations are at the same mutation-selection equilibrium frequency, denoted q . We used non-approximated values of q as derived by ( 78 ): We assume that the sites are independent. The number of mutations carried by a chromosomal segment of length n then follows a binomial distribution of parameters ( n, q ). We follow the frequency of an inversion I of size n , considering that it captures m mutations and that it appears in a population carrying only non-inverted segments. The mean fitness of a non-inverted homozygote can be computed as follows (see Appendix, section 2): Note that in the parameter regimes where we can make the approximation q ≈ u / hs with q < < 1, we have , in accordance with mutation load theory ( 37, 79 ). Similarly, an individual who is heterozygous for an inversion with m mutations has a mean fitness of: Assuming that q < < 1, Nei et al., 1967 considered that individuals heterozygous for an inversion had no mutations homozygous, such that W NI ≈ (1 − hs ) m (1 − hs ) nq ≈ (1 − hs ) nq+ m Note that we use non-approximated values in computations. An individual who is homozygous for a segment I with m mutations is homozygous for all these mutations. Its fitness can therefore be expressed as: W II = (1 − s ) m . The inversion frequency trajectory can be determined with a simple two-locus two-allele model. We considered four different situations, depending on the possible heterozygosity at the locus with permanently heterozygous alleles. The Appendix, sections 4-8, presents the evolution in time of the frequency of the inversion in detail in these four situations. Here, we briefly describe the results for inversions more or less linked to a permanently heterozygous allele in a XY system. The change in frequency of inversions on the Y chromosome, on the X chromosome or on autosomes between generations t and t + 1( Figure 2c ) is described by where is the frequency of inversions in the population of Y chromosomes at time t , is the frequency of the inversion on the X chromosome in females (respectively in males), r is the rate of recombination between the inversion and the sex-determining locus, D is their linkage disequilibrium such that and is the mean male fitness ( Figure 2c ; see Appendix, section 8 for details). When r =0.5, this system of equations describes the evolution of inversions on autosomes. Unless stated otherwise, the deterministic simulations presented here ( Figures 2c and S4-8) were performed with an initial D =-0.01 or D =0.01, depending on the allele which the inversion appeared linked to. Results for various initial linkage disequilibrium values are presented in Figure S6. When the inversion captures the male-determining allele, r =0, F XIm = F XIf = 0, and F XNf = F XNm = 1 at any time t . The equations for , the frequency of inversions capturing the male-determining allele then reduces to with . Substituting F YN by 1 − F YN and withdrawing F YN , we have 1 − F YN ,we have The search F YN for such that Δ F YN = 0 readily gives two equilibria: 0 and 1. Since and 0 ≤ F YN ≤ 1, we can conclude that: In the case of an inversion appearing on an autosome ( r =0.5), F XIm = F XIf = F YI and , so that the change in inversion frequency in the population is: and The equilibria are 0 and the roots of the polynomial P, which are 1 and . We thus have: Inversion equilibrium frequencies therefore depend, as expected, on the relative fitness of homozygotes and heterozygotes for the inversion and of the non-inverted homozygotes. We thus derive the conditions for the inversion to be favoured or disfavoured as a function of the number of mutations captured by the inversion ( m) . A straightforward computation gives (see the Appendix, section 9): with Assuming q <<1 and s <<1, these quantities can be approximated by and β 2 ≈ q . When qhn/ (1 − h ) < m < nq , inversions should thus go to fixation on Y chromosomes and stabilize at intermediate frequency on autosomes. When m < qhn/ (1 − h ), inversions should go to fixation on autosomes and on the Y chromosome. Observe that the closer h is to ½ ( i . e ., the scenario without dominance), the smaller the difference between the thresholds β 1 and β 2 is. When h=0.5, β 1 = β 2 , inversions should therefore go to fixation on autosome when they have fewer mutations than average, as they have then no homozygous disadvantage preventing their fixation. In contrast, when h is small, β 1 is significantly smaller than β 2 , showing that the condition for heterozygotes to be favoured over non-inverted homozygotes ( W NI > W NN ) is much easier to meet than for inversion homozygotes to be favoured over inversion heterozygotes ( W II > W NI ): inversions are therefore much likely to be maintained at intermediate frequencies on autosomes than to fix. See Figure S1 for a graphical representation of these results and the Appendix, section 9, for derivation details. We used these equilibrium frequencies as functions of m to compute the expected equilibrium frequency of inversions that can occur in the genome ( F equ , Figure 2d ). To do so, we sum, for all m ∈ {0, …, n } the equilibrium frequency of inversions ( F equ,m ) weighted by their occurrence probability ( P m ). For inversions on the Y chromosome, we therefore have: For inversions on autosomes, we obtain, using the approximate values for β 1 and β 2 derived earlier in the case q<<1 and s<<1: The expected equilibrium frequency of less-loaded inversions ( Figure 2d ) is therefore the expected equilibrium frequency of all inversions divided by the probability of occurrence of less-loaded inversions: The fate of inversions associated with a heterozygote fitness cost is described in Section 10 of the Appendix. Individual-based simulations Initial simulations from Jay et al. 2022 In the initial simulations from Jay et al. 2022 , we used SLiM V3.2 ( Haller & Messer, 2019 ) to simulate the evolution of a single panmictic population of N =1000 or N =10,000 individuals in a Wright-Fisher model. To assess the fate of inversions under various conditions ( Figure 3 ), we simulated individuals with two pairs of 10Mb chromosomes on which mutations occurred at a rate u , with u ranging from 10 −8 to 10 −9 per bp, their dominance coefficient h ranged from 0 to 0.5 (0, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5) and their selection coefficient s from -0.5 to -0.001 (−0.001, -0.01, -0.1, -0.25, -0.5). We also simulated populations where each mutation had its selection coefficient s drawn from a gamma distribution with a shape of 0.2, and its dominance coefficient h randomly sampled among 0, 0.001, 0.01, 0.1, 0.25, 0.5 with uniform probabilities (Figure S17). We considered recombination rates of 10 −6 and 10 −5 per bp, which gave similar results. Results are presented for analyses in which the recombination rate was 10 −6 . Among the 5,000,000 sites on chromosome 1, a single segregating locus was subject to balancing selection, for which several situations were considered: (i) the locus had two alleles, only one of which was permanently heterozygous, mimicking a classical XY (or ZW) determining system ( Figure 3 ), (ii) the locus had two permanently heterozygous alleles, mimicking, for instance, the situation encountered at most fungal mating-type loci, iii) the locus had three (or more) permanently heterozygous alleles, mimicking, for instance, the situation encountered in plant self-incompatibility systems and mushroom ( Agaromycotina ) mating-type loci. For each parameter combination ( u, h, s, N , heterozygosity rule at the locus with a permanently heterozygous allele), a simulation was run for 15,000 generations, to allow the population to reach an equilibrium for the number of segregating mutations. Figure S23 shows that each population had reached equilibrium by the end of the burn-in period. The population state was saved at the end of this initialization phase. These saved states (one for each parameter combination) were repeatedly used as initial states for studying the dynamics of recombination modifiers. Recombination modifiers mimicking inversions of 500kb, 1000kb, 2000kb and 5000kb were then introduced on chromosome 1 around the locus under balancing selection (X-linked or Y-linked inversions) or on chromosome 2 (autosomal inversions). For each parameter combination ( h, s, u, N , heterozygosity rule, size of the region affected by the recombination modification and position on the genome), we ran 10,000 independent simulations starting with the introduction of a single recombination modifier in the same saved initial population. These inversion-mimicking, recombination modifier mutations were introduced on a single, randomly selected chromosome and, when heterozygous, they suppressed recombination across the region in which they reside ( i . e ., as a cis -recombination modifier). We monitored the frequency of these inversion-mimicking mutations during 10,000 generations, during which all evolutionary processes (such as point mutation, recombination and mating) remained unchanged, e . g ., mutations were still appearing on inversions following their formation. Under the same assumptions and parameters, we also studied the dynamics of recombination modifiers suppressing recombination also when homozygous and not only when heterozygous, again across a fragment in which they reside (Figure S15). To study the effect of the existence of a haploid phase on the accumulation of deleterious mutations and the spread of inversions on autosomes and sex chromosomes, we performed additional simulations, involving 15,000 generations of burn-in and the introduction of 10,000 inversions under each combination of parameters in these initial populations, as above. The populations were considered to harbour a single locus with two permanently heterozygous alleles (similarly to the previous situation ii). Every x generation, populations experimented a haploid generation, with x taking the values 2, 3, 5, 10 or 100. For simulating haploidy, all mutations from one chromosome of each pair (“genome2” in SLiM) were removed, and the dominance coefficient of mutations on the other chromosome (“genome1” of each pair) was set to 1. The recombination rate was set to 0, and mating could only occur between gametes that were derived from the first chromosome of each pair. Therefore, during haploid generations, selection acted only on the first chromosome of each pair, and the second chromosome had no contribution to the following generation. These modifications allowed simulating the occurrence of haploid phases without changing most parameters or model behavior. Note that, during the haploid phase, the number of individuals remained unchanged, but the number of haploid genomes was divided by two, because the second haploid genome of each pair had no contribution to the following generation. To study the evolution of chromosomal fusion (Figure S17), we simulated populations with no recombination between X and Y chromosomes (chromosome 1), between the position 1Mb and 9Mb during 15,000 generations (burn-in), mimicking the evolution of a population with old sex chromosomes and small pseudo-autosomal regions (1Mb on each chromosome edge). Then, we introduced a fusion-mimicking mutation resulting in the linkage of one sex chromosome (chromosome 1, X or Y) and one autosome (chromosome 2), and suppressing the recombination when heterozygous over 1Mb of the fused side of each chromosome (see figure S18 for a graphical representation). Therefore, these mutations behave like 2Mb inversions that would also lead to chromosome fusion and result in the extension of the size of the non-recombining region from 8Mb to 10Mb. We tracked the frequency of 10,000 X-autosome and Y-autosome fusion-mimicking mutations for each parameter combination (as done before for inversion-mimicking mutations). To study more specifically the formation of evolutionary strata on sex chromosomes ( Figures 4 and 5 ), we also simulated the evolution of two 100 Mb chromosomes, one of which carried an XY sex-determining locus at the 50 Mb position, over 115,000 generations (including an initial burn-in of 15,000 generations): individuals could be either XX or XY and could only mate with individuals of a different genotype at this locus. We simulated randomly mating populations of N =1000 and N =10,000 individuals. Point mutations appeared at a rate of μ =10 −9 per bp, and their individual selection coefficients were determined by sampling a gamma distribution with a mean of -0.03 and with a shape of 0.2; these parameter values were set according to observations in humans ( Eyre-Walker & Keightley, 2007 ; Kim et al., 2017 ). For each new mutation, a dominance coefficient was chosen from the following values, considered to have uniform probabilities: 0, 0.001, 0.01, 0.1, 0.25, 0.5. At the beginning of each simulation, we randomly sampled k genomic positions over the two chromosomes that could be used as inversion breakpoints, with k being 10, 100, 1000 or 10,000. After the 15,000 generations of the burn-in period allowing populations to reach an equilibrium in terms of the number of segregating mutations, we introduced each generation j inversions in the population, j being sampled from a Poisson distribution of parameter λ , with λ = N * k * u i , u i being the inversion rate. In order to keep the simulation time tractable, we used inversion rates allowing one inversion to occur on average each generation in the population ( i . e ., N * k * u i = 1). For each inversion, the first breakpoint was randomly chosen among the k positions, and the second among the potential breakpoint positions less than 20Mb apart on the same chromosome (considering therefore only a subset of the k positions for the second breakpoints). Two independent inversions could use the same breakpoints, allowing in particular inversion reversion restoring recombination. If two independent inversions occurred with the same breakpoints on different haplotypes (for instance on the X and on the Y chromosome), we assumed that recombination was restored between these haplotypes, as a reversion would do. The occurrence of partially overlapping inversions ( i . e ., with different breakpoints) on different haplotypes did not restore recombination between these haplotypes. We assumed that inversions subsequently partially overlapped by another inversion or that captured a smaller inversion could not reverse, i . e . the restoration of recombination by further inversions, even if using the same breakpoints, was prevented. (For N =1000, for each value of k ( i . e ., 10, 100, 1000, 10000), we ran 10 simulations ( Figure 5 ). For N =10,000, because of computing limitation (each simulation taking about three weeks to run), we only ran one simulation per number of breakpoints. New simulations in this study compared to Jay et al (2022) for controlling for the effect of genetic drift To further support the conclusions from Jay et al. 2022 , we performed new simulations, in particular to compare the rate of inversion fixation on Y-like sex chromosomes and autosomes with similar population sizes ( Figure 4 ). We used SLiM V4.1 ( Haller & Messer, 2023 ) to simulate the evolution of a single panmictic population of N =3,125 or N =12,500 individuals under a Wright-Fisher model. To assess the fate of inversions under various conditions, we simulated diploid individuals with a pair of 5Mb chromosomes on which mutations occurred at a rate μ= 1.2*10 −9 per bp and recombination occurred at a rate r=1.2*10 −9 . Each new mutation had its selection coefficient s drawn from a gamma distribution with a shape of 0.2 and a mean of -0.01 or -0.001. The dominance coefficient of mutation ( h) was randomly sampled with uniform probabilities among either [0.1, 0.2, 0.3, 0.4, 0.5] (mean=0.3, no fully recessive mutations), or [0.0, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5] (mean=0.22, includes full recessive mutations). We performed simulations of autosomes and sex chromosomes. For the latter, a single segregating locus with two alleles at the chromosome center was subject to balancing selection: one of the two alleles was permanently heterozygous, mimicking a classical XY system. For each parameter combination ( h, s, N , chromosome type), a simulation was run for 200,000 generations, to allow the population to reach an equilibrium for the number of segregating mutations. At the end of this initialization phase, the nucleotide diversity of populations (π) ranged from 0.000115 (with s=-0.01 and no fully recessive mutations) and 0.000312 (with s=-0.001 and fully recessive mutations). The non-neutral (with s<-1/N) genetic diversity ranged from 0.000040 and 0.000098, indicating that diploid individuals carried on average one deleterious mutation every 10000-15000 base pair. These levels of genetic diversity in terms of deleterious mutations are in the order of magnitude of those estimated in natural populations: for instance, in humans, there is one heterozygous site every 1000bp and about 25% of segregating point mutations have been estimated to be deleterious ( Racimo & Schraiber, 2014 ). The population state was saved at the end of the initialization phase. These saved states (one for each parameter combination) were repeatedly used as initial states for studying the dynamics of chromosomal inversions. Recombination modifiers mimicking inversions of 2Mb and 5Mb were then introduced at the center of the autosome or Y chromosome. These inversions, when heterozygous, suppressed recombination across the region in which they resided. We considered inversions either with no intrinsic advantage (s inv = 0.0, meaning their fitness depended solely on the set of mutations they captured) or with an intrinsic fitness advantage (s inv = 0.01 or 0.05). For each parameter combination ( h, s, N , chromosome type, size of the chromosomal inversion, fitness effect of the inversion), we ran 100,000 independent simulations starting with the introduction of a single inversion in a single randomly-sampled individual, and we used the same saved initial population for all simulations. We monitored the frequency of these inversion-mimicking mutations during 25,000 generations, during which all evolutionary processes (such as point mutation, recombination and mating) remained unchanged, e . g ., mutations were still appearing on inversions following their formation. In order to reduce simulation time, simulations were stopped when the inversions reached fixation, i . e ., when the inversion reached a frequency of 1.0 on autosomes or 0.25 on sex chromosomes. Simulations were parallelized with GNU Parallel ( 82 ).and plot were made with ggplot2. All scripts are available at git{at}github.com :PaulYannJay/MutationShelteringTheoryV2.git. List of supplementary materials Figures S1-23 Appendix (supplementary methods) Acknowledgment All authors acknowledge the contribution of Emilie Tezenas, co-author of ( Jay et al. 2022 ), who has now left academia and is therefore not an author on this new version of the work. They also thank Ricardo Rodriguez de la Vega, Fanny Hartmann, Jacqui Shykoff, Sylvain Billiard, Janis Antonovics, Olivier Tenaillon and Diala Abu Awad for interesting discussions and comments on previous draft versions of the manuscript. AV thanks Denis Roze for insightful discussions. AV acknowledges support from the chaire program « Mathematical modelling and biodiversity » (Ecole Polytechnique, Museum National d’Histoire Naturelle, Veolia Environnement, Fondation X). Funder Information Declared International Human Frontier Science Program Organization, https://ror.org/02ebx7v45 European Research Council, https://ror.org/0472cxd90 Footnotes Funding: This work was supported by the European Research Council (ERC) EvolSexChrom (832352) grant to TG, a Louis D. Foundation (Institut de France) prize to TG, and a Human Frontier Science Program (HFSP) fellowship to PJ. Competing interests: The authors have no competing interests to declare. Data and materials availability: The SLiM and R scripts used for produce all main and supplementary figures are available from GitHub (git{at}github.com:PaulYannJay/MutationShelteringTheoryV2.git). The numerical outputs of the original simulations from Jay et al. 2022 used to produce this manuscript figures are available from Figshare ( doi:10.6084/m9.figshare.19704457 ). The output from the new simulations will be uploaded on Figshare upon publication. Details concerning the mathematical modeling are available in the appendix. References ↵ Abbott , J. K. , Nordén , A. K. , & Hansson , B. ( 2017 ). Sex chromosome evolution : Historical insights and future perspectives . Proceedings. Biological Sciences , 284 ( 1854 ), Article 1854 . doi: 10.1098/rspb.2016.2806 OpenUrl CrossRef PubMed ↵ Agrawal , A. F. , & Whitlock , M. C. ( 2011 ). Inferences about the distribution of dominance drawn from yeast gene knockout data . Genetics , 187 ( 2 ), 553 – 566 . doi: 10.1534/genetics.110.124560 OpenUrl Abstract / FREE Full Text ↵ Agrawal , A. F. , & Whitlock , M. C. ( 2012 ). Mutation load : The fitness of individuals in populations where deleterious alleles are abundant . Annual Review of Ecology Evolution and Systematics , 43 , 115 – 135 . doi: 10.1146/annurev-ecolsys-110411-160257 OpenUrl CrossRef ↵ Antonovics , J. , & Abrams , J. Y. ( 2004 ). Intratetrad mating and the evolution of linkage relationships . Evolution; International Journal of Organic Evolution , 58 ( 4 ), Article 4 . doi: 10.1111/j.0014-3820.2004.tb00403.x OpenUrl CrossRef ↵ Bachtrog , D. ( 2013 ). Y-chromosome evolution : Emerging insights into processes of Y-chromosome degeneration . Nature Reviews. Genetics , 14 ( 2 ), Article 2 . doi: 10.1038/nrg3366 OpenUrl CrossRef ↵ Badouin , H. , Hood , M. E. , Gouzy , J. , Aguileta , G. , Siguenza , S. , Perlin , M. H. , Cuomo , C. A. , Fairhead , C. , Branca , A. , & Giraud , T. ( 2015 ). Chaos of rearrangements in the mating-type chromosomes of the anther-smut fungus Microbotryum lychnidis-dioicae . Genetics , 200 ( 4 ), Article 4 . doi: 10.1534/genetics.115.177709 OpenUrl Abstract / FREE Full Text ↵ Bazzicalupo , A. L. , Carpentier , F. , Otto , S. P. , & Giraud , T. ( 2019 ). Little Evidence of Antagonistic Selection in the Evolutionary Strata of Fungal Mating-Type Chromosomes (Microbotryum lychnidis-dioicae) . G3: Genes|Genomes|Genetics , 9 ( 6 ), 1987 – 1998 . doi: 10.1534/g3.119.400242 OpenUrl Abstract / FREE Full Text ↵ Bellott , D. W. , Hughes , J. F. , Skaletsky , H. , Brown , L. G. , Pyntikova , T. , Cho , T.-J. , Koutseva , N. , Zaghlul , S. , Graves , T. , Rock , S. , Kremitzki , C. , Fulton , R. S. , Dugan , S. , Ding , Y. , Morton , D. , Khan , Z. , Lewis , L. , Buhay , C. , Wang , Q. , … Page , D. C. ( 2014 ). Mammalian Y chromosomes retain widely expressed dosage-sensitive regulators . Nature , 508 ( 7497 ), Article 7497 . doi: 10.1038/nature13206 OpenUrl CrossRef ↵ Berdan , E. L. , Barton , N. H. , Butlin , R. , Charlesworth , B. , Faria , R. , Fragata , I. , Gilbert , K. J. , Jay , P. , Kapun , M. , Lotterhos , K. E. , Mérot , C. , Durmaz Mitchell , E. , Pascual , M. , Peichel , C. L. , Rafajlovic , M. , Westram , A. M. , Schaeffer , S. W. , Johannesson , K. , & Flatt , T. ( 2023 ). How chromosomal inversions reorient the evolutionary process . Journal of Evolutionary Biology , 36 ( 12 ), Article 12 . doi: 10.1111/jeb.14242 OpenUrl CrossRef ↵ Berdan , E. L. , Blanckaert , A. , Butlin , R. K. , & Bank , C. ( 2021a ). Deleterious mutation accumulation and the long-term fate of chromosomal inversions . PLOS Genetics , 17 ( 3 ), Article 3 . doi: 10.1371/journal.pgen.1009411 OpenUrl CrossRef ↵ Berdan , E. L. , Blanckaert , A. , Butlin , R. K. , & Bank , C. ( 2021b ). Deleterious mutation accumulation and the long-term fate of chromosomal inversions . PLOS Genetics , 17 ( 3 ), e1009411 . doi: 10.1371/journal.pgen.1009411 OpenUrl CrossRef PubMed ↵ Beukeboom , L. , & Perrin , N. ( 2014 ). The Evolution of Sex Determination . Oxford University Press . ↵ Blaser , O. , Neuenschwander , S. , Perrin , N. , Gardner , A. E. A. , & Day , E. T. ( 2014 ). Sex-chromosome turnovers : The hot-potato model . The American Naturalist , 183 ( 1 ), Article 1 . doi: 10.1086/674026 OpenUrl CrossRef PubMed ↵ Boideau , F. , Richard , G. , Coriton , O. , Huteau , V. , Belser , C. , Deniot , G. , Eber , F. , Falentin , C. , Ferreira de Carvalho , J. , Gilet , M. , Lodé-Taburel , M. , Maillet , L. , Morice , J. , Trotoux , G. , Aury , J.-M. , Chèvre , A.-M. , & Rousseau-Gueutin , M. ( 2021 ). Epigenomic and structural events preclude recombination in Brassica napus . New Phytologist . doi: 10.1111/nph.18004 OpenUrl CrossRef ↵ Branco , S. , Badouin , H. , Vega , R.C.R.dela , Gouzy , J. , Carpentier , F. , Aguileta , G. , Siguenza , S. , Brandenburg , J.-T. , Coelho , M. A. , Hood , M. E. , & Giraud , T. ( 2017 ). Evolutionary strata on young mating-type chromosomes despite the lack of sexual antagonism . Proceedings of the National Academy of Sciences , 114 ( 27 ), Article 27 . doi: 10.1073/pnas.1701658114 OpenUrl Abstract / FREE Full Text ↵ Branco , S. , Carpentier , F. , Rodríguez de la Vega , R. C., Badouin , H. , Snirc , A. , Le Prieur , S. , Coelho , M. A. , de Vienne , D. M. , Hartmann , F. E. , Begerow , D. , Hood , M. E. , & Giraud , T. ( 2018 ). Multiple convergent supergene evolution events in mating-type chromosomes . Nature Communications , 9 ( 1 ), Article 1 . doi: 10.1038/s41467-018-04380-9 OpenUrl CrossRef PubMed ↵ Carey , S. B. , Jenkins , J. , Lovell , J. T. , Maumus , F. , Sreedasyam , A. , Payton , A. C. , Shu , S. , Tiley , G. P. , Fernandez-Pozo , N. , Healey , A. , Barry , K. , Chen , C. , Wang , M. , Lipzen , A. , Daum , C. , Saski , C. A. , McBreen , J. C. , Conrad , R. E. , Kollar , L. M. , … McDaniel , S. F. ( 2021 ). Gene-rich UV sex chromosomes harbor conserved regulators of sexual development . Science Advances , 7 ( 27 ), Article 27 . doi: 10.1126/sciadv.abh2488 OpenUrl CrossRef ↵ Carpentier , F. , Rodríguez de la Vega , R.C. , Jay , P. , Duhamel , M. , Shykoff , J. A. , Perlin , M. H. , Wallen , R. M. , Hood , M. E. , & Giraud , T. ( 2022 ). Tempo of Degeneration Across Independently Evolved Nonrecombining Regions . Molecular Biology and Evolution , 39 ( 4 ), Article 4 . doi: 10.1093/molbev/msac060 OpenUrl CrossRef ↵ Cavoto , E. , Neuenschwander , S. , Goudet , J. , & Perrin , N. ( 2018 ). Sex-antagonistic genes, XY recombination and feminized Y chromosomes . Journal of Evolutionary Biology , 31 ( 3 ), Article 3 . doi: 10.1111/jeb.13235 OpenUrl CrossRef ↵ Charlesworth , B. , Coyne , J. A. , & Barton , N. H. ( 1987 ). The relative rates of evolution of sex chromosomes and autosomes . The American Naturalist , 130 ( 1 ), Article 1 . doi: 10.1086/284701 OpenUrl CrossRef Web of Science ↵ Charlesworth , B. , & Olito , C. ( 2024 ). Making sense of recent models of the “sheltering” hypothesis for recombination arrest between sex chromosomes . Evolution , 78 ( 12 ), 1891 – 1899 . doi: 10.1093/evolut/qpae119 OpenUrl CrossRef PubMed ↵ Charlesworth , B. , & Wall , J. D. ( 1999 ). Inbreeding, heterozygote advantage and the evolution of neo–X and neo–Y sex chromosomes . Proceedings of the Royal Society of London. Series B: Biological Sciences , 266 ( 1414 ), Article 1414 . doi: 10.1098/rspb.1999.0603 OpenUrl CrossRef ↵ Coelho , S. , Gueno , J. , Lipinska , A. , Cock , J. , & Umen , J. ( 2018 ). UV chromosomes and haploid sexual systems . Trends in plant science , 23 ( 9 ), Article 9 . doi: 10.1016/j.tplants.2018.06.005 OpenUrl CrossRef ↵ Connallon , T. , & Olito , C. ( 2021 ). Natural selection and the distribution of chromosomal inversion lengths . Molecular Ecology . doi: 10.1111/mec.16091 OpenUrl CrossRef ↵ Connallon , T. , Olito , C. , Dutoit , L. , Papoli , H. , Ruzicka , F. , & Yong , L. ( 2018 ). Local adaptation and the evolution of inversions on sex chromosomes and autosomes . Philosophical Transactions of the Royal Society B: Biological Sciences , 373 ( 1757 ), Article 1757 . doi: 10.1098/rstb.2017.0423 OpenUrl CrossRef ↵ Cortez , D. , Marin , R. , Toledo-Flores , D. , Froidevaux , L. , Liechti , A. , Waters , P. D. , Grützner , F. , & Kaessmann , H. ( 2014 ). Origins and functional evolution of Y chromosomes across mammals . Nature , 508 ( 7497 ), 488 – 493 . doi: 10.1038/nature13151 OpenUrl CrossRef PubMed Web of Science ↵ Cui , L. , Neoh , H. , Iwamoto , A. , & Hiramatsu , K. ( 2012 ). Coordinated phenotype switching with large-scale chromosome flip-flop inversion observed in bacteria . Proceedings of the National Academy of Sciences , 109 ( 25 ), Article 25 . doi: 10.1073/pnas.1204307109 OpenUrl Abstract / FREE Full Text ↵ De Filippo , E. , Chahine , E. , Legendre--Despas , J. , Snirc , A. , Labat , A. , Michel , P. , Grognet , P. , Gautier , V. , Levert , E. , Lalanne , C. , Silar , P. , Giraud , T. , & Hartmann , F. E. ( 2025 ). Stepwise recombination suppression around the mating-type locus associated with a diploid-like life cycle in Schizothecium fungi (p. 2025.04.18.645578). bioRxiv . doi: 10.1101/2025.04.18.645578 OpenUrl Abstract / FREE Full Text ↵ De-Kayne , R. , Gordon , I. J. , Terblanche , R. F. , Collins , S. , Omufwoko , K. S. , Martins , D. J. , & Martin , S. H. ( 2025 ). Incomplete recombination suppression fuels extensive haplotype diversity in a butterfly colour pattern supergene . PLOS Biology , 23 ( 2 ), e3003043 . doi: 10.1371/journal.pbio.3003043 OpenUrl CrossRef PubMed ↵ de Vienne , D. M. , Giraud , T. , & Gouyon , P.-H. ( 2013 ). Lineage Selection and the Maintenance of Sex . PLOS ONE , 8 ( 6 ), e66906 . doi: 10.1371/journal.pone.0066906 OpenUrl CrossRef PubMed ↵ Dobzhansky , T. ( 1972 ). Genetics of the Evolutionary Process : Columbia Classics edition . Columbia University Press . ↵ Dyer , K. A. , Charlesworth , B. , & Jaenike , J. ( 2007 ). Chromosome-wide linkage disequilibrium as a consequence of meiotic drive . Proceedings of the National Academy of Sciences , 104 ( 5 ), Article 5 . doi: 10.1073/pnas.0605578104 OpenUrl Abstract / FREE Full Text ↵ Eyre-Walker , A. , & Keightley , P. D. ( 2007 ). The distribution of fitness effects of new mutations . Nature Reviews Genetics , 8 ( 8 ), Article 8 . doi: 10.1038/nrg2146 OpenUrl CrossRef Web of Science ↵ Ferchaud , A.-L. , Mérot , C. , Normandeau , E. , Ragoussis , J. , Babin , C. , Djambazian , H. , Bérubé , P. , Audet , C. , Treble , M. , Walkusz , W. , & Bernatchez , L. ( 2022 ). Chromosome-level assembly reveals a putative Y-autosomal fusion in the sex determination system of the Greenland Halibut (Reinhardtius hippoglossoides) . G3 Genes\textbarGenomes\textbarGenetics , 12 ( 1 ), Article 1 . doi: 10.1093/g3journal/jkab376 OpenUrl CrossRef ↵ Filatov , D. A. ( 2015 ). Homomorphic plant sex chromosomes are coming of age . Molecular Ecology , 24 ( 13 ), Article 13 . doi: 10.1111/mec.13268 OpenUrl CrossRef ↵ Giner-Delgado , C. , Villatoro , S. , Lerga-Jaso , J. , Gayà-Vidal , M. , Oliva , M. , Castellano , D. , Pantano , L. , Bitarello , B. D. , Izquierdo , D. , Noguera , I. , Olalde , I. , Delprat , A. , Blancher , A. , Lalueza-Fox , C. , Esko , T. , O’Reilly , P. F. , Andrés , A. M. , Ferretti , L. , Puig , M. , & Cáceres , M. ( 2019 ). Evolutionary and functional impact of common polymorphic inversions in the human genome . Nature Communications , 10 ( 1 ), Article 1 . doi: 10.1038/s41467-019-12173-x OpenUrl CrossRef PubMed ↵ Grossen , C. , Neuenschwander , S. , & Perrin , N. ( 2012 ). The evolution of xy recombination : Sexually antagonistic selection versus deleterious mutation load . Evolution , 66 ( 10 ), Article 10 . doi: 10.1111/j.1558-5646.2012.01661.x OpenUrl CrossRef ↵ Haller , B. C. , & Messer , P. W. ( 2019 ). SLiM 3 : Forward genetic simulations beyond the wright–fisher model . Molecular Biology and Evolution , 36 ( 3 ), Article 3 . doi: 10.1093/molbev/msy228 OpenUrl CrossRef PubMed ↵ Haller , B. C. , & Messer , P. W. ( 2023 ). SLiM 4 : Multispecies Eco-Evolutionary Modeling . The American Naturalist , 201 ( 5 ), E127 – E139 . doi: 10.1086/723601 OpenUrl CrossRef PubMed ↵ Hämälä , T. , Wafula , E. K. , Guiltinan , M. J. , Ralph , P. E. , dePamphilis , C. W. , & Tiffin , P. ( 2021 ). Genomic structural variants constrain and facilitate adaptation in natural populations of Theobroma cacao, the chocolate tree . Proceedings of the National Academy of Sciences , 118 ( 35 ), Article 35 . doi: 10.1073/pnas.2102914118 OpenUrl CrossRef ↵ Hanson , S. J. , Byrne , K. P. , & Wolfe , K. H. ( 2014 ). Mating-type switching by chromosomal inversion in methylotrophic yeasts suggests an origin for the three-locus Saccharomyces cerevisiae system . Proceedings of the National Academy of Sciences of the United States of America , 111 ( 45 ), Article 45 . doi: 10.1073/pnas.1416014111 OpenUrl Abstract / FREE Full Text ↵ Hartmann , F. E. , Duhamel , M. , Carpentier , F. , Hood , M. E. , Foulongne-Oriol , M. , Silar , P. , Malagnac , F. , Grognet , P. , & Giraud , T. ( 2021 ). Recombination suppression and evolutionary strata around mating-type loci in fungi : Documenting patterns and understanding evolutionary and mechanistic causes . New Phytologist , 229 ( 5 ), Article 5 . doi: 10.1111/nph.17039 OpenUrl CrossRef PubMed ↵ Hartmann , F. E. , Rodríguez de la Vega , R. C., Demené , A. , Badet , T. , Vernadet , J.-P. , Rougemont , Q. , Labat , A. , Snirc , A. , Stauber , L. , Croll , D. , Prospero , S. , Dutech , C. , & Giraud , T. ( 2025 ). An Inversion Polymorphism Under Balancing Selection, Involving Giant Mobile Elements, in an Invasive Fungal Pathogen . Molecular Biology and Evolution , 42 ( 2 ), msaf026 . doi: 10.1093/molbev/msaf026 OpenUrl CrossRef PubMed ↵ Ironside , J. E. ( 2010 ). No amicable divorce? Challenging the notion that sexual antagonism drives sex chromosome evolution . BioEssays , 32 ( 8 ), Article 8 . doi: 10.1002/bies.200900124 OpenUrl CrossRef ↵ Jay , P. , Chouteau , M. , Whibley , A. , Bastide , H. , Parrinello , H. , Llaurens , V. , & Joron , M. ( 2021 ). Mutation load at a mimicry supergene sheds new light on the evolution of inversion polymorphisms . Nature Genetics , 53 ( 3 ), Article 3 . doi: 10.1038/s41588-020-00771-1 OpenUrl CrossRef ↵ Jay , P. , Jeffries , D. , Hartmann , F. E. , Véber , A. , & Giraud , T. ( 2024 ). Why do sex chromosomes progressively lose recombination? Trends in Genetics: TIG , 40 ( 7 ), 564 – 579 . doi: 10.1016/j.tig.2024.03.005 OpenUrl CrossRef PubMed ↵ Jay , P. , Tezenas , E. , Véber , A. , & Giraud , T. ( 2022 ). Sheltering of deleterious mutations explains the stepwise extension of recombination suppression on sex chromosomes and other supergenes . PLOS Biology , 20 ( 7 ), Article 7 . doi: 10.1371/journal.pbio.3001698 OpenUrl CrossRef ↵ Kim , B. Y. , Huber , C. D. , & Lohmueller , K. E. ( 2017 ). Inference of the distribution of selection coefficients for new nonsynonymous mutations using large samples . Genetics , 206 ( 1 ), Article 1 . doi: 10.1534/genetics.116.197145 OpenUrl Abstract / FREE Full Text ↵ Kimura , M. , & Ohta , T. ( 1969 ). The average number of generations until fixation of a mutant gene in a finite population . Genetics , 61 ( 3 ), Article 3 . OpenUrl ↵ Kirkpatrick , M. ( 2010 ). How and why chromosome inversions evolve . PLoS Biology , 8 ( 9 ), Article 9 . doi: 10.1371/journal.pbio.1000501 OpenUrl CrossRef ↵ Kondrashov , A. S. , & Crow , J. F. ( 1991 ). Haploidy or diploidy : Which is better? Nature , 351 ( 6324 ), Article 6324 . doi: 10.1038/351314a0 OpenUrl CrossRef ↵ Kubisiak , T. L. , & Milgroom , M. G. ( 2006 ). Markers linked to vegetative incompatibility (vic) genes and a region of high heterogeneity and reduced recombination near the mating type locus (MAT) in Cryphonectria parasitica . Fungal Genetics and Biology: FG & B , 43 ( 6 ), Article 6 . doi: 10.1016/j.fgb.2006.02.002 OpenUrl CrossRef ↵ Küpper , C. , Stocks , M. , Risse , J. E. , dos Remedios , N. , Farrell , L. L. , McRae , S. B. , Morgan , T. C. , Karlionova , N. , Pinchuk , P. , Verkuil , Y. I. , Kitaysky , A. S. , Wingfield , J. C. , Piersma , T. , Zeng , K. , Slate , J. , Blaxter , M. , Lank , D. B. , & Burke , T. ( 2016 ). A supergene determines highly divergent male reproductive morphs in the ruff . Nature Genetics , 48 ( 1 ), Article 1 . doi: 10.1038/ng.3443 OpenUrl CrossRef ↵ Le Veve , A. , Burghgraeve , N. , Genete , M. , Lepers-Blassiau , C. , Takou , M. , Meaux , J. D. , Mable , B. K. , Durand , E. , Vekemans , X. , & Castric , V. ( 2022 ). Long-term balancing selection and the genetic load linked to the self-incompatibility locus in Arabidopsis halleri and A. lyrata (p. 2022.04.12.487987). bioRxiv . doi: 10.1101/2022.04.12.487987 OpenUrl Abstract / FREE Full Text ↵ Lemaitre , C. , Braga , M. D. V. , Gautier , C. , Sagot , M.-F. , Tannier , E. , & Marais , G. A. B. ( 2009 ). Footprints of inversions at present and past pseudoautosomal boundaries in human sex chromosomes . Genome Biology and Evolution , 1 , 56 – 66 . doi: 10.1093/gbe/evp006 OpenUrl CrossRef PubMed ↵ Lenormand , T. , & Roze , D. ( 2022 ). Y recombination arrest and degeneration in the absence of sexual dimorphism . Science , 375 ( 6581 ), Article 6581 . doi: 10.1126/science.abj1813 OpenUrl CrossRef ↵ Lenormand , T. , & Roze , D. ( 2024 ). Can mechanistic constraints on recombination reestablishment explain the long-term maintenance of degenerate sex chromosomes? Peer Community Journal , 4 . doi: 10.24072/pcjournal.373 OpenUrl CrossRef ↵ Llaurens , V. , Billiard , S. , Castric , V. , & Vekemans , X. ( 2009 ). Evolution of Dominance in Sporophytic Self-Incompatibility Systems : I. Genetic Load and Coevolution of Levels of Dominance in Pollen and Pistil . Evolution , 63 ( 9 ), Article 9 . doi: 10.1111/j.1558-5646.2009.00709.x OpenUrl CrossRef Manna , F. , Gallet , R. , Martin , G. , & Lenormand , T. ( 2012 ). The high-throughput yeast deletion fitness data and the theories of dominance . Journal of Evolutionary Biology , 25 ( 5 ), 892 – 903 . doi: 10.1111/j.1420-9101.2012.02483.x OpenUrl CrossRef PubMed ↵ Manna , F. , Martin , G. , & Lenormand , T. ( 2011 ). Fitness landscapes : An alternative theory for the dominance of mutation . Genetics , 189 ( 3 ), 923 – 937 . doi: 10.1534/genetics.111.132944 OpenUrl Abstract / FREE Full Text ↵ Mrnjavac , A. , Vicoso , B. , & Connallon , T. ( 2025 ). A re-evaluation of Muller’s sheltering hypothesis for the evolution of sex chromosome gene content (p. 2025.04.24.650552). bioRxiv . doi: 10.1101/2025.04.24.650552 OpenUrl Abstract / FREE Full Text ↵ Nei , M. ( 1970 ). Accumulation of nonfunctional genes on sheltered chromosomes . The American Naturalist , 104 ( 938 ), Article 938 . OpenUrl ↵ Nei , M. , Kojima , K.-I. , & Schaffer , H. E. ( 1967 ). Frequency changes of new inversions in populations under mutation-selection equilibria . Genetics , 57 ( 4 ), Article 4 . OpenUrl ↵ Ohta , T. ( 1971 ). Associative overdominance caused by linked detrimental mutations . Genetics Research , 18 ( 3 ), Article 3 . doi: 10.1017/S0016672300012684 OpenUrl CrossRef ↵ Olito , C. , & Abbott , J. K. ( 2020 ). The evolution of suppressed recombination between sex chromosomes by chromosomal inversions . bioRxiv . doi: 10.1101/2020.03.23.003558 OpenUrl Abstract / FREE Full Text ↵ Olito , C. , & Abbott , J. K. ( 2025 ). The evolution of suppressed recombination between sex chromosomes and the lengths of evolutionary strata . Evolution, qpaf045 . doi: 10.1093/evolut/qpaf045 OpenUrl CrossRef ↵ Olito , C. , & Charlesworth , B. ( 2023 ). Do deleterious mutations promote the evolution of recombination suppression between X and Y chromosomes? (p. 2023.11.27.568803) . bioRxiv . doi: 10.1101/2023.11.27.568803 OpenUrl Abstract / FREE Full Text ↵ Pennell , M. W. , Kirkpatrick , M. , Otto , S. P. , Vamosi , J. C. , Peichel , C. L. , Valenzuela , N. , & Kitano , J. ( 2015 ). Y Fuse? Sex Chromosome Fusions in Fishes and Reptiles . PLOS Genetics , 11 ( 5 ), Article 5 . doi: 10.1371/journal.pgen.1005237 OpenUrl CrossRef ↵ Ponnikas , S. , Sigeman , H. , Abbott , J. K. , & Hansson , B. ( 2018 ). Why do sex chromosomes stop recombining? Trends in Genetics , 34 ( 7 ), Article 7 . doi: 10.1016/j.tig.2018.04.001 OpenUrl CrossRef ↵ Porubsky , D. , Höps , W. , Ashraf , H. , Hsieh , P. , Rodriguez-Martin , B. , Yilmaz , F. , Ebler , J. , Hallast , P. , Maria Maggiolini , F. A. , Harvey , W. T. , Henning , B. , Audano , P. A. , Gordon , D. S. , Ebert , P. , Hasenfeld , P. , Benito , E. , Zhu , Q. , Human Genome Structural Variation Consortium (HGSVC), Lee , C. , … Korbel , J. O. ( 2022 ). Recurrent inversion polymorphisms in humans associate with genetic instability and genomic disorders . Cell , 185 ( 11 ), Article 11 . doi: 10.1016/j.cell.2022.04.017 OpenUrl CrossRef ↵ Racimo , F. , & Schraiber , J. G. ( 2014 ). Approximation to the Distribution of Fitness Effects across Functional Categories in Human Segregating Polymorphisms . PLOS Genetics , 10 ( 11 ), Article 11 . doi: 10.1371/journal.pgen.1004697 OpenUrl CrossRef ↵ Reinhardt , J. A. , Brand , C. L. , Paczolt , K. A. , Johns , P. M. , Baker , R. H. , & Wilkinson , G. S. ( 2014 ). Meiotic Drive Impacts Expression and Evolution of X-Linked Genes in Stalk-Eyed Flies . PLOS Genetics , 10 ( 5 ), Article 5 . doi: 10.1371/journal.pgen.1004362 OpenUrl CrossRef ↵ Rice , W. R. ( 1987 ). The accumulation of sexually antagonistic genes as a selective agent promoting the evolution of reduced recombination between primitive sex chromosomes . Evolution , 41 ( 4 ), Article 4 . doi: 10.2307/2408899 OpenUrl CrossRef ↵ Ruzicka , F. , Dutoit , L. , Czuppon , P. , Jordan , C. Y. , Li , X.-Y. , Olito , C. , Runemark , A. , Svensson , E. I. , Yazdi , H. P. , & Connallon , T. ( 2020 ). The search for sexually antagonistic genes : Practical insights from studies of local adaptation and statistical genomics . Evolution Letters , 4 ( 5 ), Article 5 . doi: 10.1002/evl3.192 OpenUrl CrossRef ↵ Saunders , P. A. , & Muyle , A. ( 2024 ). Sex Chromosome Evolution : Hallmarks and Question Marks . Molecular Biology and Evolution , 41 ( 11 ), msae218 . doi: 10.1093/molbev/msae218 OpenUrl CrossRef PubMed ↵ Scott , M. F. , & Rescan , M. ( 2017 ). Evolution of haploid–diploid life cycles when haploid and diploid fitnesses are not equal . Evolution , 71 ( 2 ), Article 2 . OpenUrl ↵ Skinner , B. M. , Sargent , C. A. , Churcher , C. , Hunt , T. , Herrero , J. , Loveland , J. E. , Dunn , M. , Louzada , S. , Fu , B. , Chow , W. , Gilbert , J. , Austin-Guest , S. , Beal , K. , Carvalho-Silva , D. , Cheng , W. , Gordon , D. , Grafham , D. , Hardy , M. , Harley , J. , … Tyler-Smith , C. ( 2021 , juin 8). The pig X and Y Chromosomes : Structure, sequence, and evolution . doi: 10.1101/gr.188839.114 OpenUrl Abstract / FREE Full Text ↵ Stauber , L. , Badet , T. , Feurtey , A. , Prospero , S. , & Croll , D. ( 2021 ). Emergence and diversification of a highly invasive chestnut pathogen lineage across southeastern Europe . eLife , 10 , e56279 . doi: 10.7554/eLife.56279 OpenUrl CrossRef PubMed ↵ Sun , Y. , Svedberg , J. , Hiltunen , M. , Corcoran , P. , & Johannesson , H. ( 2017 ). Large-scale suppression of recombination predates genomic rearrangements in Neurospora tetrasperma . Nature Communications , 8 ( 1 ), 1140 . doi: 10.1038/s41467-017-01317-6 OpenUrl CrossRef PubMed ↵ Takayama , S. , & Isogai , A. ( 2005 ). Self-incompatibility in plants . Annual Review of Plant Biology , 56 , 467 – 489 . doi: 10.1146/annurev.arplant.56.032604.144249 OpenUrl CrossRef PubMed ↵ Todesco , M. , Owens , G. L. , Bercovich , N. , Légaré , J.-S. , Soudi , S. , Burge , D. O. , Huang , K. , Ostevik , K. L. , Drummond , E. B. M. , Imerovski , I. , Lande , K. , Pascual-Robles , M. A. , Nanavati , M. , Jahani , M. , Cheung , W. , Staton , S. E. , Muños , S. , Nielsen , R. , Donovan , L. A. , … Rieseberg , L. H. ( 2020 ). Massive haplotypes underlie ecotypic differentiation in sunflowers . Nature , 584 ( 7822 ), Article 7822 . doi: 10.1038/s41586-020-2467-6 OpenUrl CrossRef PubMed ↵ Villoutreix , R. , Ayala , D. , Joron , M. , Gompert , Z. , Feder , J. L. , & Nosil , P. ( 2021 ). Inversion breakpoints and the evolution of supergenes . Molecular Ecology , 30 ( 12 ), Article 12 . doi: 10.1111/mec.15907 OpenUrl CrossRef ↵ Vittorelli , N. , Rodríguez de la Vega , R. C., Snirc , A. , Levert , E. , Gautier , V. , Lalanne , C. , De Filippo , E. , Gladieux , P. , Guillou , S. , Zhang , Y. , Tejomurthula , S. , Grigoriev , I. V. , Debuchy , R. , Silar , P. , Giraud , T. , & Hartmann , F. E. ( 2023 ). Stepwise recombination suppression around the mating-type locus in an ascomycete fungus with self-fertile spores . PLoS Genetics , 19 ( 2 ), e1010347 . doi: 10.1371/journal.pgen.1010347 OpenUrl CrossRef PubMed ↵ Wang , J. , Wurm , Y. , Nipitwattanaphon , M. , Riba-Grognuz , O. , Huang , Y.-C. , Shoemaker , D. , & Keller , L. ( 2013 ). A Y-like social chromosome causes alternative colony organization in fire ants . Nature , 493 ( 7434 ), Article 7434 . doi: 10.1038/nature11832 OpenUrl CrossRef ↵ Wellenreuther , M. , & Bernatchez , L. ( 2018 ). Eco-evolutionary genomics of chromosomal inversions . Trends in Ecology & Evolution , 33 ( 6 ), Article 6 . doi: 10.1016/j.tree.2018.04.002 OpenUrl CrossRef ↵ Wright , A. E. , Darolti , I. , Bloch , N. I. , Oostra , V. , Sandkam , B. , Buechel , S. D. , Kolm , N. , Breden , F. , Vicoso , B. , & Mank , J. E. ( 2017 ). Convergent recombination suppression suggests role of sexual selection in guppy sex chromosome formation . Nature Communications , 8 ( 1 ), Article 1 . doi: 10.1038/ncomms14251 OpenUrl CrossRef PubMed ↵ Wright , A. E. , Dean , R. , Zimmer , F. , & Mank , J. E. ( 2016 ). How to make a sex chromosome . Nature Communications , 7 ( 1 ), Article 1 . doi: 10.1038/ncomms12087 OpenUrl CrossRef ↵ Yan , Z. , Martin , S. H. , Gotzek , D. , Arsenault , S. V. , Duchen , P. , Helleu , Q. , Riba-Grognuz , O. , Hunt , B. G. , Salamin , N. , Shoemaker , D. , Ross , K. G. , & Keller , L. ( 2020 ). Evolution of a supergene that regulates a trans-species social polymorphism . Nature Ecology & Evolution , 4 ( 2 ), Article 2 . doi: 10.1038/s41559-019-1081-1 OpenUrl CrossRef PubMed ↵ Yashiro , T. , Tea , Y.-K. , Van Der Wal , C. , Nozaki , T. , Mizumoto , N. , Hellemans , S. , Matsuura , K. , & Lo , N. ( 2021 ). Enhanced heterozygosity from male meiotic chromosome chains is superseded by hybrid female asexuality in termites . Proceedings of the National Academy of Sciences , 118 ( 51 ), Article 51 . doi: 10.1073/pnas.2009533118 OpenUrl CrossRef ↵ Zhou , Y. , Minio , A. , Massonnet , M. , Solares , E. , Lv , Y. , Beridze , T. , Cantu , D. , & Gaut , B. S. ( 2019 ). The population genetics of structural variants in grapevine domestication . Nature Plants , 5 ( 9 ), Article 9 . doi: 10.1038/s41477-019-0507-8 OpenUrl CrossRef PubMed ↵ Zhou , Y. , Zhan , X. , Jin , J. , Zhou , L. , Bergman , J. , Li , X. , Rousselle , M. M. C. , Belles , M. R. , Zhao , L. , Fang , M. , Chen , J. , Fang , Q. , Kuderna , L. , Marques-Bonet , T. , Kitayama , H. , Hayakawa , T. , Yao , Y.-G. , Yang , H. , Cooper , D. N. , … Zhang , G. ( 2023 ). Eighty million years of rapid evolution of the primate Y chromosome . Nature Ecology & Evolution , 7 ( 7 ), Article 7 . doi: 10.1038/s41559-022-01974-x OpenUrl CrossRef View the discussion thread. Back to top Previous Next Posted July 01, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Stepwise expansion of recombination suppression on sex chromosomes and other supergenes through lower load advantage and deleterious mutation sheltering Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Stepwise expansion of recombination suppression on sex chromosomes and other supergenes through lower load advantage and deleterious mutation sheltering Paul Jay , Amandine Véber , Tatiana Giraud bioRxiv 2025.06.27.661902; doi: https://doi.org/10.1101/2025.06.27.661902 Share This Article: Copy Citation Tools Stepwise expansion of recombination suppression on sex chromosomes and other supergenes through lower load advantage and deleterious mutation sheltering Paul Jay , Amandine Véber , Tatiana Giraud bioRxiv 2025.06.27.661902; doi: https://doi.org/10.1101/2025.06.27.661902 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Evolutionary Biology Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17691) Bioengineering (13892) Bioinformatics (41937) Biophysics (21452) Cancer Biology (18589) Cell Biology (25504) Clinical Trials (138) Developmental Biology (13378) Ecology (19899) Epidemiology (2067) Evolutionary Biology (24320) Genetics (15609) Genomics (22506) Immunology (17736) Microbiology (40394) Molecular Biology (17181) Neuroscience (88605) Paleontology (666) Pathology (2832) Pharmacology and Toxicology (4824) Physiology (7641) Plant Biology (15156) Scientific Communication and Education (2045) Synthetic Biology (4294) Systems Biology (9825) Zoology (2271)