Consequences of partially recessive deleterious genetic variation for the evolution of inversions suppressing recombination between sex chromosomes

Abstract The evolution of suppressed recombination between sex chromosomes is widely hypothesized to be driven by sexually antagonistic selection (SA), where tighter linkage between the sex‐determining gene(s) and nearby SA loci is favored when it couples male‐beneficial alleles to the proto‐Y chromosome, and female‐beneficial alleles to the proto‐X. Despite limited empirical evidence, the SA selection hypothesis overshadows several alternatives, including an incomplete but often‐repeated “sheltering hypothesis” that suggests that expansion of the sex‐linked region (SLR) reduces homozygous expression of partially recessive deleterious mutations at selected loci. Here, we use population genetic models to evaluate the consequences of deleterious mutational variation for the evolution of neutral chromosomal inversions expanding the SLR on proto‐Y chromosomes. We find that SLR‐expanding inversions face a race against time: lightly loaded inversions are initially beneficial, but eventually become deleterious as they accumulate new mutations, and must fix before this window of opportunity closes. The outcome of this race is strongly influenced by inversion size, the mutation rate, and the dominance coefficient of deleterious mutations. Yet, small inversions have elevated fixation probabilities relative to neutral expectations for biologically plausible parameter values. Our results demonstrate that deleterious genetic variation can plausibly drive recombination suppression in small steps and would be most consistent with empirical patterns of small evolutionary strata or gradual recombination arrest.

Sex chromosomes have evolved from homologous pairs of autosomes repeatedly within many eukaryotic lineages across the tree of life (reviewed by Beukeboom & Perrin 2014, Bachtrog et al. 2014. A striking feature of many sex chromosome systems is the evolution of recombination suppression, which profoundly influences the long-term fate of the chromosomes. Once recombination stops, the subsequent evolution of sequence divergence and functional degeneration of the non-recombining region of the sex-limited chromosome, and possibly dosage compensation to retain adequate gene expression levels in both sexes, can all contribute to the gradual evolution of sex chromosome heteromorphy (Bull 1983, Bachtrog 2006, Beukeboom & Perrin 2014, Charlesworth et al. 2005a, Lenormand & Roze 2022, Rice 1987, 1996.
The initial loss of recombination between sex chromosomes is widely hypothesized to be caused by selection favoring linkage disequilibrium between the sex-determining gene(s) and nearby loci experiencing sex-differences in selection (e.g., sexually antagonistic loci with alleles that have opposite fitness effect between sexes) (Bull 1983, Charlesworth & Charlesworth 1980, Lenormand 2003, Otto 2019, Rice 1987, 1996. Indeed, even though strong empirical support for the sexual antagonism hypothesis remains elusive, it continues to overshadow a variety of alternatives that have received less theoretical or empirical attention (Jefferies et al. 2021, Lenormand & Roze 2022, Olito & Abbott 2020, Ponnikas et al. 2018. Several of these alternative hypotheses revolve around the idea that a chromosomal rearrangement -typically an inversion -expanding the male-limited region of a Y chromosome (or female-limited region of a W) may be selectively favored due to a form of heterozygote advantage, or "sheltering," arising from the combination of wild-type and partially recessive deleterious alleles that it captures (Branco et al. 2017, Ironside 2010, Jay et al. 2021, Ponnikas et al. 2018. In fact, a tangle of at least three distinct sheltering hypotheses have been described in varying detail, ranging from loose verbal models to mathematical and simulation models. First, Charlesworth & Wall (1999) cautiously suggested that selection might favor linkage between the sex-determining locus and multiple selected loci exhibiting associative overdominance due to partially recessive deleterious variation, an idea that was echoed in subsequent reviews (Ironside 2010, Ponnikas et al. 2018. A later verbal model by Branco et al. (2017) proposed that partially recessive deleterious mutations in partial linkage with the sex-determining region could cause selection for recombination arrest to avoid homozygous expression in recombinant genotypes. Most recently, Jay et al. (2021) modeled the evolution of inversions expanding the sex-linked region on Y chromosomes under deleterious mutation pressure using constant fitness effects for inversion genotypes depending on the alleles they initially capture, with accompanying individual-based simulations (we review the various hypotheses in detail in Appendix A of the Supporting Information). Each of these models proposes, in some way, that an inversion linking alleles at selected loci to the heterozygous male-determining allele on a Y chromosome will reduce the homozygous expression of partially recessive deleterious alleles at those loci. This sheltering effect is hypothesized to cause higher fitness for inverted relative to non-inverted Y chromosomes.
The red thread running through each of these sheltering hypotheses is that deleterious mutational variation is pervasive throughout the genome (Crow 1970, Charlesworth et al. 1993, Muller 1950) and the fate of new inversions expanding the sexlinked region may therefore be strongly influenced by the random sample of that variation that they happen to capture. Moreover, deleterious genetic variation is known to have important implications for the evolution of inversions on autosomes: autosomal inversions undergo a complex time-dependent selection process that can result in diverse evolutionary outcomes, including fixation, extinction, and balanced polymorphism, depending on the set of alleles they initially capture (Charlesworth & Charlesworth 1973, Connallon & Olito 2020, Nei et al. 1967. However, there is a crucial difference between autosomal inversions and those expanding the sex-linked region on a Y chromosome. By capturing the dominant sex-determining factor (or expanding the chromosomal region already linked to it), the latter are prevented from occurring in both X and Y chromosomes. At first glance, this might appear to enforce heterozygosity of partially recessive deleterious alleles captured by the inversion, thereby paving the way for inversion fixation as suggested by earlier studies (Ironside 2010, Jay et al. 2021). As we demonstrate below, this is an oversimplification. In fact, when partially recessive deleterious genetic variation is present, inversions expanding the sex-linked region (SLR hereafter) on Y chromosomes experience fundamentally different time-dependent selection processes compared to autosomal ones. Careful consideration of these time-inhomogeneous selection processes is necessary to fully understand the evolutionary dynamics of inversions contributing to recombination suppression between sex chromosomes.
In this paper, we develop a population genetic model to describe the evolutionary dynamics of new inversion mutations that capture the dominant sex-determining factor in randomly mating populations, while explicitly considering the consequences of standing deleterious mutational variation. We briefly describe the deterministic frequency dynamics predicted by the model before turning our attention to the fixation probabilities for inversions of different lengths, which we calculate using stochastic Wright-Fisher simulations. Our results illuminate two important features of inversions on Y chromosomes expanding the SLR: (i) because partially recessive deleterious mutations segregate outside the SLR on both X and Y chromosomes, mutations initially captured by an inversion are expressed as homozygotes at a rate equal to their frequency in X chromosomes, so that inversions capturing even a single deleterious allele will carry a permanent, but temporally dynamic, deleterious mutation load; (ii) inversions initially capturing fewer than the average number of deleterious mutations over the chromosomal segment they span will initially be beneficial, but this selective advantage erodes over time as new mutations accumulate on descendent copies of the inversion until the benefit becomes smaller than the permanent load carried by the inversion; at this point the overall fitness effect of the inversion becomes irreversibly deleterious and its chances of fixation negligible. Hence, the evolutionary fate of SLR-expanding inversions is a race against time; initially beneficial inversions must fix before their selective advantage decays and their window of opportunity closes permanently. The outcome of this race can be non-intuitive and is determined jointly by the size of new inversions as well as several key population genetic parameters, including the deleterious mutation rate, dominance and selection parameters, and population size. We close by discussing the implications of our findings for existing theories of recombination arrest between sex-chromosomes.

OVERVIEW OF THE MODEL
Consider a population of diploid, randomly mating individuals with discrete generations, in which sex is determined genetically by a dominant male-determining factor (as in a maleheterogametic X-Y system). The model is equally applicable to female heterogametic Z-W systems if male/female labels are reversed. The order of life history events proceeds as follows: fertilization, mutation, selection, then meiosis. We model idealized sex chromosomes that can be divided into two regions: (i) a non-recombining SLR that could be limited to just the sexdetermining (SD) gene(s), or a more extensive region harboring the SD gene(s); and (ii) a pseudoautosomal region (PAR) in which recombination can still occur, and functional homologs of any genes are still present on both X and Y chromosomes. Hence, our models are most applicable to genetic systems in which the evolution of recombination suppression between sex chromosomes is incomplete, and the PAR accounts for a sizeable fraction of the proto sex chromosomes (e.g., Charlesworth et al. 2005b).
We model the evolution of new chromosomal inversion mutations arising on a Y chromosome that would expand the SLR were they to fix among the Y chromosomes in the population. Our goal is to predict how new inversions will respond to indirect selection against deleterious mutations segregating within the population at the loci they span. To isolate these indirect selection effects, we assume that the inversion itself is neutral (i.e., inversions cause no breakpoint effects or meiotic dysfunction; Corbett-Detig 2016, Krimbas & Powell 1992, Olito & Abbott 2020, Villoutreix et al. 2021. For simplicity, we also assume that loci within the SLR do not contribute to indirect selection on inversions. This second assumption can be justified in different ways: (i) there has been sufficient differentiation and functional degeneration within the SLR on Y chromosomes that few functional genes remain in this region; (ii) the SLR is small relative to the length of inversions, such that there are few loci other than the sex-determining genes within the SLR; and (iii) any loci within the SLR that are captured by an inversion will contribute minimally to indirect selection favoring suppressed recombination because they are already fully sex-linked.
Our model relies on several other important simplifying assumptions. First, we assume that inversions completely suppress recombination between inverted-Y and X chromosomes over the chromosomal region they span (in fact, genetic exchange could occur via double crossovers or gene conversion; Krimbas & Pow-ell 1992, Korunes & Noor 2019). Second, we assume that new inversions occur rarely enough that the evolutionary fate of a given inversion is independent of others (i.e., we assume "strong selection, weak mutation" with respect to inversions; Gillespie 1991). Our results therefore preclude the possibility that multiple inversions segregate simultaneously within the population. Third, we assume that deleterious alleles segregate at mutation-selection balance outside of the SLR, with no epistasis, and no linkage disequilibrium among loci or with the SLR prior to a new inversion mutation. This requires strong purifying selection against deleterious variants relative to mutation or genetic drift. Finally, we assume that fitness is multiplicative over the loci spanned by the inversion.
Below, we develop a deterministic model describing the frequency dynamics of new inversions expanding the SLR on Y chromosomes in the presence of partially recessive deleterious mutational variation and illustrate important features of the model predictions. We emphasize the role of joint changes in the inversion frequency as well as deleterious allele frequencies on X chromosomes. For simplicity, we present deterministic results for the idealized case where mutation and selection coefficients at all loci spanned by the inversion are equal. We then present Wright-Fisher (W-F) simulations that incorporate stochastic fluctuations in inversion frequencies due to random gamete sampling in a finite population. We present most of the mathematical details in Appendix B. Computer code needed to reproduce the simulations and main figures is available on GitHub (https://github. com/colin-olito/shelteringOnSexChrom) and a version of record is archived on Zenodo (Olito et al. 2022).

DETERMINISTIC MODEL
We first define two useful terms: the total number of selected loci located outside of the SLR on the chromosome arm containing it (i.e., within the PAR), n tot , and the length of a new inversion, x, expressed as the proportion of the PAR that the new inversion links to the ancestral SLR. Assuming that the selected loci are distributed uniformly along the chromosome arm, the number of loci spanned by a new inversion of length x will be n = n tot x. Each of the n loci are assumed to be diallelic, with a wild-type allele (A) at the i th locus that mutates to a deleterious variant (a) at a rate µ i per meiosis (we ignore backmutation from a → A), with locus-specific genotypic relative fitnesses of w i,AA = 1, w i,Aa = 1 − h i s i , w i,aa = 1 − s i . Deleterious alleles segregate at each locus at their mutation-selection balance equilibrium frequency ofq i = µ i /(h i s i ). A new inversion mutation will capture a random sample of the standing deleterious variation at these n loci, which can be conveniently divided into two classes: loci where the inversion initially captures a deleterious allele (denoted by a superscript D) and those where it captures a wild-type allele (denoted by a superscript W ).

B R I E F C O M M U N I C AT I O N
To describe the deterministic frequency dynamics of an SLR-expanding inversion, it is necessary to track allele frequency changes for these two classes of loci, D and W , within four different chromosome classes: X's in ovules/eggs, X's in pollen/sperm, non-inverted Y's, and inverted Y's, which we denote X f , X m , Y , and Y I respectively (see Otto 2014, Olito & Abbott 2020. The frequencies at time t for each of the n loci can be described using the following notation: where q refers to the deleterious allele frequency at the i th locus. Note that, by assumption, q D,i Y I = 1 for all t because all descendent copies of the inversion will carry deleterious alleles at D loci. Likewise, q W,i Y I = 0 for t = 1 but increases due to new mutations on descendent copies of the inversion. We present the full development of the recursions in Appendix B. Under the simplifying assumption that mutation and selection parameters are constant across all loci captured by the inversion (µ i = µ, s i = s, h i = h), the number of deleterious alleles initially captured by the inversion, r, will be Poisson distributed: r ∼ Poisson(U x/hs). In this idealized case, the deleterious allele frequencies at all loci within a given chromosomal class will follow the same trajectory (e.g., each inversion-captured locus that is initially free of a deleterious allele will evolve the same as other loci in the inversion with the same, mutation-free initial state). That is, q W,i class,t = q W class,t , and q D,i This allows us to define the following simplified recursion for the frequency of an inversion that initially captures r deleterious alleles in terms of the allele frequencies: and we have used the convention p · ·,t = 1 − q · ·,t to simplify notation. The deterministic frequency dynamics of the inversion can now be fully described by a system of eight recursions corresponding to the deleterious allele frequencies in each of the seven ) and the frequency of the inversion, Equations (1) and (2).
Equation (1) offers immediate insight into the different contributions of D vs. W loci to relative fitnesses of inversion genotypes. The fitness effects of D loci depend solely on the frequency of deleterious alleles in X chromosomes in ovules/eggs (selection terms in bracketed expressions with an exponent of r involve only p D Xf ,t and q D Xf ,t ) because all descendent copies of the inversion already carry a deleterious allele at these loci. Under recurrent mutation, q D Xf ,t will always be nonzero, and so it is immediately clear that D loci will impart a permanent fitness cost to the inversion. Meanwhile, W loci depend jointly on the frequency of deleterious alleles in ovule/egg-derived X chromosomes and the accumulation of new deleterious mutations on descendent copies of the inversion at these loci (selection terms involve p W Xf ,t , q W Xf ,t , and q W Y I ,t ). Compared to the fitness of standard Y chromosomes (the second expression in square brackets in Equation 2) it is clear that any temporary fitness advantage of new inversions must come from the initially low frequency of deleterious alleles at W loci among inverted Y chromosomes (q W Y I ,t ). The key questions become: when, and for how long, does the temporary fitness benefit from W loci outweigh the permanent load associated with D loci? and does it result in elevated fixation probabilities for new SLR-expanding inversions?

DETERMINISTIC FREQUENCY DYNAMICS
When deleterious alleles are approximately codominant (i.e., h i ≈ 1/2), most purifying selection occurs in heterozygotes and an inversion initially loaded with even a single deleterious allele (i.e., when r > 0) will not invade (see Connallon & Olito 2020, Olito & Abbott 2020. However, when deleterious mutations are partially recessive (0 < h i < 1/2), as expected by theory and supported by empirical data (Agrawal & Whitlock 2011, Huber et al. 2018, Manna et al. 2011reviewed in Billiard et al. 2021), lightly loaded inversions expanding the SLR on Y chromosomes can deterministically rise to high frequency under restrictive conditions. that are initially beneficial because they capture fewer deleterious mutations than the population average over the chromosomal segment they span. The time course of inversion fitness relative to non-inverted Y chromosomes (Fig. 1a-c) illustrates both the decay of the initial fitness benefit as new mutations accumulate on descendent copies of the inversion-bearing chromosome at W loci (inversions that initially have relative fitness greater than one), and the permanent deleterious load due to D loci (all inversions with r > 0 eventually become deleterious). The tipping-point where the relative fitness of loaded inversions drops below one occurs when the transient benefit to the inversion of initially capturing fewer-than-average deleterious mutations no longer compensates for the cost of being fixed for those few mutations. The deterministic inversion frequency dynamics (Fig. 1d-f) reflect these changes in relative fitness and highlight that the deterministic outcome for initially unloaded inversions is to fix among Y chromosomes in the population. When deleterious mutations are strongly recessive (h ≈ 0.01), inversions capturing several deleterious mutations can deterministically rise to intermediate to high frequencies before the decline in their initial fitness advantage results in them becoming deleterious and crashing to extinction (Fig. 1c,f,i). The restrictive conditions under which an initially loaded inversion is expected to deterministically 'fix' occur when all deleterious mutations are strongly recessive and a large inversion captures very few of them (see Figs. S1-S3 in Appendix C).
Importantly, the load carried by inversions due to D loci is temporally dynamic because the frequencies of partially recessive deleterious alleles on X (and non-inverted Y) chromosomes change over time in response to the inversion frequency ( Fig. 1G-I). We illustrate this for representative cases of inversions initially loaded with relatively few deleterious mutations (r = 1 in panels G,H; r = 10 in panel I) by showing the deterministic frequency dynamics for three important loci × chromosome classes: W loci on the inversion (q W YI ), and both W and D loci on X chromosomes in ovules/eggs (q W Xf and q D Xf ). The red line shows the accumulation of deleterious mutations on the inversion at W loci, which causes the decline in initial fitness benefit due to these loci. The black solid and dotted lines show the corresponding changes in deleterious allele frequency among ovule/egg-derived X chromosomes at W loci (q W Xf ) and D loci (q D Xf ), respectively. Selection against deleterious alleles at D loci on ovule/egg-derived X chromosomes intensifies as the initially beneficial inversion increases in frequency because all inverted-Y-bearing sons who inherit a deleterious allele from their mother's X chromosome will be homozygous for the deleterious allele at D loci. This intensified selection drives the frequency of deleterious alleles at D loci on X chromosomes (q D Xf ) down to lower levels relative to the equilibrium prior to the inversion. Nevertheless, the deleterious load due to D loci persists because new mutations continually arise on X and non-inverted Y chromosomes. When the transient benefit of the inversion due to W loci can no longer compensate for this load, the inversion becomes deleterious and declines in frequency, at which point the deleterious allele frequencies return to the pre-inversion equilibrium. These dynamics are especially evident when deleterious mutations are strongly recessive (Fig. 1c,f,i).
An overview of the deterministic dynamics for inversions of different lengths initially loaded with different numbers of deleterious alleles is presented in Appendix C, Figures S1-S3. Of particular note in the Supporting Information figures is that lightly loaded large inversions have higher initial fitness than smaller ones, but are still not expected to increase to high frequencies unless deleterious mutations are strongly recessive. The time dependent dynamics can also result in faster accumulation of strongly recessive deleterious mutations on SLR-expanding inversions (Fig. S4).

WRIGHT-FISHER SIMULATIONS
While the deterministic frequency dynamics presented above clearly illustrate the shifting balance between the time-dependent selection processes on inversions due to W and D loci, they do not consider either the stochastic process of gamete sampling present in finite populations that can result in loss of rare inversions, nor the likelihood that a new inversion captures a given number of mutations. In reality, larger inversions are more likely to capture a greater number of deleterious mutations (see Figs. S1-S3), and the combined effects of W and D loci on inversion relative fitness will depend on allele frequency dynamics at each of the n loci. To begin tackling the effects of drift, we use Wright-Fisher simulations carried out in R (R Core Team 2020) to estimate the fixation probability for a single Y chromosome bearing an SLR-expanding inversion as a function of inversion length. We make the (strong) simplifying assumption that deleterious allele frequencies at the n loci spanned by the inversion change deterministically, while the inversion itself is subject to genetic drift due to random gamete sampling. This approach qualitatively captures the effects of time-dependent selection on inversion fixation probability in large populations, where ephemeral indirect selection effects are most likely to be important. However, it also makes our simulation results not applicable to smaller populations where deleterious allele frequency dynamics are dominated by drift rather than selection. Unfortunately, relaxing this assumption using an individual based model becomes computationally intractable for relevant values of selection and dominance coefficients due to the high replication necessary to reliably estimate fixation probabilities (see below).
In each generation, haplotype frequencies are censused among gametes prior to fertilization, following mutation and indirect selection due to segregating deleterious alleles at each locus, with a total effective population size of N. We use our exact deterministic recursions (Appendix B) to generate predictions for the gene frequencies at each of the n loci after indirect selection for each locus × chromosome class. The realized frequency of the inversion among Y chromosomes in each generation was then simulated using pseudo-random binomial sampling during the gametic phase, with the number of Y chromosomes representing the number of trials (N/2), and the deterministic frequency predictions representing the probabilities of sampling the inversion among pollen/sperm. For each replicate simulation, the number of loci spanned by each new inversion was n = n tot x and the number of deleterious alleles initially captured by a new inversion of length x was drawn from a Poisson distribution with mean and variance U x/hs. Fixation probabilities were estimated from the outcomes of at least 100 × N/2 replicate simulations (this was increased to 500 × N/2 for larger inversions because of their low fixation probabilities). For comparison, we also estimated the fixation probability of autosomal inversions using a similar procedure but implementing multinomial pseudo-random sampling of adult genotypes (e.g., Charlesworth & Charlesworth 2010, pp. 229-230).

INVERSION FIXATION PROBABILITIES
We focus on the effect of inversion length on fixation probability because the fitness benefits and costs scale differently with inversion length for autosomal and SLR-expanding inversions. In both cases, the initial fitness benefit of capturing wild-type alleles at more loci increases with inversion length, but so does the probability of capturing more deleterious alleles and carrying a larger permanent deleterious load. For autosomal inversions, these countervailing effects of inversion size cancel out when deleterious mutations are not strongly recessive and population sizes are sufficiently large, resulting in expected fixation probabilities that are independent of inversion length and approximately equal to the initial frequency of the inversion (Connallon & Olito 2020). With substantial (albeit partially recessive) fitness effects in heterozygotes and small population sizes, there is a fixation bias toward small autosomal inversions, which have a maximum fixation probability approximately equal to that of a typical neutral allele, 1/2N (Fig. 2a). In larger populations, the fixation probability becomes increasingly independent of inversion size, particularly when there are fewer average deleterious mutations per standard chromosome arm (Fig. 2b). We focus our analysis on deleterious alleles with equal dominance coefficients of h i = h = 0.25, which corresponds roughly to the average dominance coefficient of deleterious mutations estimated from empirical studies (Manna et al. 2011, Agrawal & Whitlock 2012reviewed by Billiard et al. 2021), but explore the effects of strong recessivity in Appendix D (see Figs. S5 and S6).
For inversions expanding the SLR on Y chromosomes, the fitness costs scale differently with inversion size for reasons that were highlighted by our deterministic model: inversions capturing fewer than the average number of deleterious alleles are initially beneficial and can therefore rise in frequency until they accumulate enough deleterious mutations that they eventually become deleterious. However, because the capture of even a single deleterious allele imparts a permanent deleterious load on a new inversion, initially unloaded inversions benefit most from the time-dependent selection process. Consequently, there is a strong fixation bias toward smaller SLR-expanding inversions (Fig. 2c,d; up to two orders of magnitude larger than the neutral expectation of 2/N) because these are most likely to initially capture no deleterious alleles (see Figs. S1-S3). Nevertheless, larger SLR-expanding inversions that are lightly loaded are still more likely to go to fixation than similarly sized autosomal ones.
The strength of mutation relative to selection, which influences the average number of deleterious mutations per standard (non-inverted) arrangement chromosome, has a similar effect on the fixation probability for autosomal versus SLR-expanding inversions. A higher average number of deleterious mutations on non-inverted chromosomes (higher U/hs) decreases the fixation probability of larger inversions (Fig. 2) because they are more likely to initially capture deleterious mutations and will accumulate new mutations more rapidly, despite having a greater initial fitness benefit over non-inverted chromosomes. The av-erage deleterious load also influences the threshold inversion size where the fixation probability of SLR-expanding inversions drops below 2/N, with higher average numbers of deleterious mutations per standard-arrangement chromosome (large U/hs) corresponding to smaller threshold sizes (Fig. 2c,d).
Despite the involvement of rather complicated timedependent selection processes, the different evolutionary dynamics of autosomal vs. SLR-expanding inversions in our models can be explained rather simply: the homozygous expression of partially recessive deleterious mutations initially captured by autosomal inversions increases with inversion frequency, while for SLR-expanding inversions they are expressed at a rate equal to their frequency on X chromosomes in ovules/eggs, which remains small (≤q i ) and even decreases with inversion frequency. Due to their lower, temporally dynamic, permanent deleterious load, lightly loaded SLR-expanding inversions have a window of time during which they can fix before their fitness benefit decays and they become deleterious. Initially unloaded inversions benefit most from the time-dependent selection process, resulting in elevated fixation probabilities for smaller inversions that can be significantly higher than that of a neutral allele.

Discussion
Our model predictions have several important implications. First, our deterministic model clarifies that the notion proposed by several earlier sheltering hypotheses-that linkage to the permanently heterozygous male-determining factor will prevent or reduce the homozygous expression of deleterious mutations, thereby favoring recombination suppression-is an oversimplification. Deleterious alleles initially captured by an SLRexpanding inversion are not prevented from being expressed as homozygotes, as suggested by earlier verbal and mathematical arguments (Ironside 2010, Jay et al. 2021. Rather, they are expressed as homozygotes at a rate equal to their frequencies on X chromosomes in ovules/eggs, which change with inversion frequency. The deleterious load carried by these inversions is temporally dynamic but still permanent: initially beneficial inversions capturing even one mutation will eventually become deleterious. This is the critical difference between SLR-expanding and autosomal inversions, for which the homozygous expression of any captured deleterious alleles increases with inversion frequency because more individuals are homozygous for the inverted karyotype, which effectively prevents even lightly loaded inversions from ever fixing (Connallon & Olito 2020, Nei et al. 1967. Second, our simulation results suggest that small inversions on proto-Y chromosomes are more likely to contribute to recombination suppression when selection is indirect and attributable to associations between inversions and deleterious variation. The elevated fixation probabilities of small inversions predicted by our models suggest that a relatively simple alternative to the sexually antagonistic selection hypothesis-a suitably located neutral inversion combined with segregating deleterious variationcan drive the evolution of recombination suppression between sex chromosomes in small steps. Whether our model predictions are consistent with empirical patterns of evolutionary strata size remains to be seen, but our results suggest that the length of SLR-expanding inversions may offer insight into the selective processes driving their fixation. In particular, large evolutionary strata generally appear inconsistent with the inversion fixation probabilites expected under the sheltering hypothesis (Connallon & Olito 2020, Santos 1986, Van Valen & Levins 1968).
It remains difficult to determine the relative importance of different hypotheses for suppressed recombination between sex chromosomes for several reasons. First, knowledge about the mutation rate and length distribution of new inversions in natural populations is limited, as is the number of well-described evolutionary strata (although genomic data are starting to shed light on these; Connallon & Olito 2020). Second, our models make several simplifying assumptions. Generalizing to other cases is beyond the scope of the present paper, but important questions remain for future work. For example, allowing non-uniform distributions of dominance and selection coefficients will likely further reduce the size of SLR-expanding inversions that ultimately become fixed, since this is the case for autosomal inversions (Connallon & Olito 2020). Yet, a few loci segregating for highly recessive deleterious alleles (e.g., h i ≈ 0.01) could have the opposite effect, although these may generally have strongly deleterious effects that would oppose this (Huber et al. 2018, Billiard et al. 2021; see Fig. S6). Lastly, our W-F simulations made the strong simplifying assumption that deleterious allele frequencies changed deterministically, while inversions were subject to genetic drift. Although the time-dependent indirect selection effects in our model will be most important in large populations, additional simulation studies are needed to determine at what population sizes they become relevant. Yet, a recent simulation study of Y recombination arrest and sex-specific regulatory evolution sheds some light on the robustness of our model predictions (Lenormand & Roze 2022). In the absence of regulatory evolution (with only partially recessive deleterious mutations, similar to the senario we investigate here), Lenormand & Roze (2022) found higher fitness variance of large SLR-expanding inversions, but faster degradation of large inversions and a fixation bias towards small inversions when N = 10 4 , results that are consistent with our model predictions.
Most notably, we have limited our attention to neutral inversions capturing deleterious mutations, without considering those with beneficial effects or capturing loci under sexually antagonistic selection that could cause linkage disequilibrium with the SLR. However, the presence of partially recessive deleterious ge-netic variation should influence the fixation probabilities of differently sized inversions under these other selection scenarios (e.g., Olito & Abbott 2020) by making the conditions for fixation of small inversions easier to satisfy. We have also focused our analysis on single randomly mating populations and selected loci initially at linkage equilibrium with the SLR. However, both inbreeding and prior linkage disequilibrium with the SLR feature prominently in previous sheltering hypotheses (Branco et al. 2017, Charlesworth & Wall 1999; see Appendix A for further details). We briefly address these possibilities using mathematical models in Appendices D and E, where we show that they do not result in indirect positive selection for recombination-suppressing inversions. However, a more complete treatment of these problems is probably warranted.
Overall, our results highlight the importance of confronting seemingly intuitive hypotheses and verbal arguments with mathematical models. Despite the intuitive appeal of the sexually antagonistic selection hypothesis for recombination suppression between sex chromosomes, our predictions suggest that a simple alternative-the presence of partially recessive deleterious genetic variation on proto-sex chromosomes and a suitably located inversion-is plausible under some conditions and deserves more careful consideration. Yet, our models also show that apparently simple verbal arguments, like the earlier sheltering hypotheses for recombination suppression, also deserve careful scrutiny because the details of the genetic system and emergent selection processes are anything but intuitive.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Figure S1: Overview of deterministic fitness and frequency dynamics for initially beneficial inversions of different sizes initially loaded with different numbers of deleterious alleles Figure S2: Overview of deterministic fitness and frequency dynamics for initially beneficial inversions of different sizes initially loaded with different numbers of deleterious alleles. Figure S3: Overview of deterministic fitness and frequency dynamics for initially beneficial inversions of different sizes initially loaded with different numbers of strongly recessive deleterious alleles (h = 0.01). Figure S4: Comparison of deterministic trajectories of deleterious allele frequencies at W loci on descendent copies of autosomal (q Y I,Aut o W ; black lines) and SLR-expanding (q W Y I , red lines) inversions, divided by their corresponding equilibrium frequency (q). Figure S5: Fixation probabilities of new SLR-expanding inversions on Y chromosomes plotted as a function of inversion length when deleterious mutations are more strongly recessive (hi = 0.1) than the average dominance coefficient of deleterious mutations estimated from empirical studies (Manna et al. 2011;Agrawal and Whitlock 2012;Huber et al. 2018;reviewed in Billiard et al. 2021). Figure S6: Fixation probabilities of new SLR-expanding inversions on Y chromosomes plotted as a function of inversion length when deleterious mutations are strongly recessive (hi = 0.01). Figure D1: Reproduction of Fig. 1A from Charlesworth and Wall (1999) using the 3-locus model of an SLR-expanding inversion Figure