The fitness consequences of genetic divergence between polymorphic gene arrangements

Abstract Inversions restrict recombination when heterozygous with standard arrangements, but often have few noticeable phenotypic effects. Nevertheless, there are several examples of inversions that can be maintained polymorphic by strong selection under laboratory conditions. A long-standing model for the source of such selection is divergence between arrangements with respect to recessive or partially recessive deleterious mutations, resulting in a selective advantage to heterokaryotypic individuals over homokaryotypes. This paper uses a combination of analytical and numerical methods to investigate this model, for the simple case of an autosomal inversion with multiple independent nucleotide sites subject to mildly deleterious mutations. A complete lack of recombination in heterokaryotypes is assumed, as well as constancy of the frequency of the inversion over space and time. It is shown that a significantly higher mutational load will develop for the less frequent arrangement. A selective advantage to heterokaryotypes is only expected when the two alternative arrangements are nearly equal in frequency, so that their mutational loads are very similar in size. The effects of some Drosophila pseudoobscura polymorphic inversions on fitness traits seem to be too large to be explained by this process, although it may contribute to some of the observed effects. Several population genomic statistics can provide evidence for signatures of a reduced efficacy of selection associated with the rarer of two arrangements, but there is currently little published data that are relevant to the theoretical predictions.


Introduction
obtained evidence that some natural inversion polymorphisms in Drosophila pseudoobscura are associated with major differences in fitness among karyotypes, which can lead to their stable maintenance within a single population under constant environmental conditions.There have subsequently been many other experimental studies documenting strong effects of inversion karyotypes on fitness components in several Drosophila species (reviewed in Krimbas and Powell 1992;Kapun and Flatt 2019), and in some other species such as the seaweed fly Coelopa frigida (Butlin et al. 1984;Mérot et al. 2020).The startling observations of Wright and Dobzhansky (1946) raised the question of the causes of fitness differences between apparently functionally insignificant chromosomal variants.This question is still the subject of ongoing inquiry, stimulated by the new evidence from genome sequencing that inversion polymorphisms are abundant in natural populations of many species (Wellenreuther and Bernatchez 2018;Faria et al. 2019;Berdan et al. 2023).
Well before the work of Wright and Dobzhansky, Sturtevant and Mather (1938) had proposed a process that could cause fitness differences between inversion karyotypes, with a fitness advantage to heterokaryotypes over homokaryotypes.In their words: … if a chromosome exists, in a population, in two sequences, differing by an inversion, it will in effect show two distinct lines of descent.There is free exchange of material within any line (i.e.sequence), but none between the sequences.Therefore, fluctuations of the genic contents must occur almost independently in the two sequences.Under such conditions, it is inevitable that in time the gene content of the two sequences will become different.It must be supposed that each sequence is susceptible to the same mutations, and with the same frequencies, but, as a result of the relative rarity in the population of a given mutant allelomorph at any one moment, certain genes will be present in one sequence but not the other.…..It is thus clear, considering two sequences A and B, that the homozygotes AA and BB are more likely to be homozygous for deleterious recessive mutations than is the heterozygote AB. ….
Thus in general the sequence heterozygote AB will be at a selective advantage with respect to either of the corresponding homozygotes.Sturtevant and Mather (1938) did not attempt a quantitative model of this process, simply noting that "there can be no stability in the exact relations of the sequences with each other", and that "…. a single gene difference can not in general cause such heterosis.The simplest effective condition is that in which each sequence contains a deleterious recessive not present in the https://doi.org/10.1093/genetics/iyad218Advance Access Publication Date: 26 December 2023 Investigation other."This proposal raises the question of what strength of selection on the two arrangements might be expected on its basis.Ohta (1971) developed a mathematical model of a closely related, but not identical, process.This was based on the concept of associative overdominance, first outlined by Frydenberg (1963).Here, a polymorphic neutral locus can acquire an apparent heterozygote advantage, because of linkage disequilibrium (LD) generated by genetic drift with a locus subject to selection in favor of heterozygotes or to selection against recessive/partially recessive deleterious alleles maintained by mutation pressure.Unlike the model of Sturtevant and Mather (1938), this process does not require the generation of heterosis by multiple selected loci, often referred to as pseudo-overdominance (Waller 2021).Ohta (1971) summed the effects of individual selected loci that were completely linked to a dialellic neutral locus (equivalent to an inversion polymorphism) over a large segment of genome, and generated expressions for the apparent fitness advantage to heterozygotes at the neutral locus.
It was, however, shown by Zhao and Charlesworth (2016) that Ohta's formulae for the relative fitnesses at the neutral locus induced by LD with the selected locus do not predict any change in allele frequency at the neutral locus.Using a different approach, they found that an induced selection pressure in favor of increased variability at the neutral locus only exists when the product of population size and selection coefficient is of the order of 1. Otherwise, variability is reduced by background selection effects, even for recessive or partially recessive deleterious mutations (see also Charlesworth 2022).Ohta's results therefore do not solve the quantitative problem of whether the strength of selection on inversions revealed by the experiments cited above can be explained by this process.Nonetheless, they are still invoked as a potential contributor to the selective maintenance of inversion polymorphisms, e.g.Faria et al. (2019), Berdan et al. (2021), Jay et al. (2021), and Matschiner et al. (2023).
A different perspective was developed by Nei et al. (1967), who examined a purely deterministic model involving the balance between mutation and selection at numerous autosomal loci.This process results in an equilibrium frequency distribution of the number of mutant alleles per haploid genome.A new autosomal inversion then has a chance of arising on a haplotype with a lower number of mutations than average and acquiring a selective advantage.But, as time goes on, the mutant-free loci on the inverted haplotypes will accumulate mutations (Nei et al. 1967).Unless the inversion goes to fixation, the loci in the inversion subpopulation will eventually acquire the same frequencies of mutant alleles as the corresponding loci in the standard arrangement subpopulation.Nei et al. (1967) interpreted this as implying that the inversion would then be selectively neutral.However, if reverse mutations from mutant to wild-type alleles do not occur, the loci at which mutations were present in the original inversion haplotype will all be homozygous in inversion homokaryotypes, causing these to have a reduced fitness compared to homokaryotypes for the standard arrangement.A reanalysis and extension of this model by Connallon and Olito (2021) suggested that it could produce a net heterozygote advantage to an autosomal inversion in a randomly mating population, resulting in the maintenance of the inversion, but only if the inversion has a sufficiently large direct selective (but nonheterotic) advantage that is independent of the deleterious mutations that it carries.Berdan et al. (2021) conducted simulations of a finite population with multiple loci experiencing mutations to deleterious and completely recessive alleles.They found that a selective advantage to heterokaryotypes could develop, provided that recombinational exchange between arrangements in heterokaryotypes was sufficiently infrequent, and there was a small additional selective advantage to heterokaryotypes that kept the inversion in the population long enough for mutation accumulation to occur.Studies of the effects of deleterious mutations have shown that complete recessivity is unlikely to be frequent, especially for mildly deleterious mutations (Muller 1950;Crow 1993;Manna et al. 2011).Multilocus computer simulations with less extreme assumptions about the degree of recessivity of deleterious mutations showed that autosomal inversion polymorphisms are unlikely to be established at higher than neutral rates (Jay et al. 2022).
These theoretical studies therefore suggest that deleterious mutations in themselves are unlikely to create an initial selective advantage to autosomal inversions in a randomly mating population.Indeed, if a new inversion arises on a unique haplotype, the process of accumulation of new mutational load within the inversion subpopulation will take a long time, and cannot contribute to any initial selective effect of the inversion.It is, however, an open question as to whether the fitness differences among karyotypes mentioned above could have a significant component resulting from the process proposed by Sturtevant and Mather (1938), whereby genetic drift causes the inverted and standard arrangements to differ in their genetic content.There is a strong, but not perfect, analogy with population subdivision, where genetic drift can cause local populations to diverge at weakly selected loci subject to mutation to deleterious variants.In the case of inversions, however, the existence of heterokaryotypes means that selection does not act independently on the two subpopulations, as described below in the section on the mathematical model.
Provided that mutations are at least partially recessive with respect to their fitness effects, interpopulation crosses may show heterosis, due to different loci having accumulated different deleterious mutations in different populations (Whitlock et al. 2000;Glémin et al. 2003;Roze and Rousset 2004;Spigler et al. 2017;Charlesworth 2018).Such heterosis has been observed in populations of animals and plants (reviewed in Charlesworth 2018).In addition, theory predicts that populations with a smaller effective size (N e ) should allow the accumulation of deleterious mutations and hence acquire larger mutational loads than populations with a larger N e , unless mutations have strong selective effects relative to drift and are highly recessive (Wright 1931;Kimura et al. 1963;Nei 1968;Bataillon and Kirkpatrick 2000;Charlesworth 2018).
Such differences in N e can arise either from differences in the adult population size itself, differences in the mating system (especially the frequency of inbreeding), differences in recombination rates associated with different levels of genetic hitchhiking effects, or a combination of all 3 factors.Again, there is empirical support for this predicted effect of N e , both from measurements of fitness components (Leimu et al. 2006;Lohr and Haag 2015;Charlesworth 2018) and from population genomic indicators of the efficacy of selection against deleterious mutations, e.g.Robinson et al. (2022) (population size), Glémin et al. (2006) (mating system), and Campos et al. (2014) (recombination rate).
The purpose of the present paper is to investigate the properties of a population genetic model of mutation, selection, and drift acting on a polymorphism for an autosomal inversion and a standard arrangement, in which the inversion polymorphism is maintained by sufficiently strong selection that the frequency of the inversion is constant over time and space.Similar assumptions were used in an earlier paper that examined neutral differentiation between the two arrangements (Charlesworth 2023).In order to simplify the calculations, and to maximize the effects of drift and mutation within arrangements, no recombinational exchange between arrangements in heterokaryotypes is allowed.The results should therefore provide upper bounds on the likely size of effects.
Models are developed of both a single, randomly mating, population and of a population divided into a large number of local populations.It might be expected that population subdivision, with its greater opportunities for drift, would enhance divergence among arrangements at sites under selection.The model assumes a large number of freely recombining sites subject to selection and mutation, with a wide distribution of selection coefficients over sites.No attempt is made to investigate the consequences of enhanced Hill-Robertson interference among sites due to restricted recombination in heterokaryotypes, which was included in the simulations of the fates of autosomal inversions in Berdan et al. (2021) and Jay et al. (2022).Unless an autosomal inversion is very rare or very common, there should be sufficient recombination within homokaryotypes to prevent major Hill-Robertson effects; very low recombination rates are sufficient to prevent the operation of Muller's ratchet (Charlesworth et al. 1993).
The overall conclusion is that a low-frequency arrangement will have a higher mutational load and exhibit weaker population genomic signals of purifying selection than its counterpart.Heterokaryotypic superiority in fitness is, however, unlikely to be observed unless the inverted and standard sequences are approximately equal in frequency, and it is likely to be small in magnitude unless the inversion contains millions of sites under selection.Population subdivision has only a small effect on the load and population genomic statistics.

General considerations and notation
The main symbols used in this paper are defined in Table 1.Consider first a single randomly mating, discrete generation population of size N, assuming a Wright-Fisher model of reproduction such that N is equal to the effective population size.The frequencies of the two karyotypes, the inverted (In) and standard (St) arrangements are denoted by x and y = 1 -x; designation as In vs St is purely arbitrary in the situation considered here, so the convention x ≤ ½ is used.Balancing selection on the inversion is assumed to be sufficiently strong that x can be treated as constant over time.Let q i and p i = 1 -q i be the respective frequencies of the mutant (A 2 ) and wild-type allele (A 1 ) at a given nucleotide site within haplotypes carrying a type i karyotype, where i = 1 corresponds to In and i = 2 to St.In a given generation, there will be random drift as well as selection within the populations of In and St karyotypes, so that in general q 1 ≠ q 2 .Drift occurs independently within karyotypes under the assumption of constant sizes of the two subpopulations, so that the effective population sizes of carriers of In and St are N 1 = Nx and N 2 = Ny, respectively.
The expectation of q i is denoted by 〈q i 〉, with 〈p i 〉 = 1 -〈q i 〉, where the angle brackets denote an expectation taken over the probability distribution of q i .The variance in q i in a given generation over this probability is denoted by V q i = 〈q 2 i 〉 − 〈q i 〉 2 , with the corresponding F statistic (Wright 1951) given by F i = V qi 〈p i 〉〈q i 〉.In the absence of recombination but the presence of selection, there will be a negative covariance C 12 between q 1 and q 2 , with a correspondingly negative correlation coefficient R 12 .This is because a higher frequency of the mutant allele in one karyotype results in a higher frequency of mutant homozygotes in In/St individuals, enhancing the strength of selection against the mutation in the other karyotype.For neutral sites in the absence of recombination between the arrangements, C 12 = R 12 = 0.
Following Kimura et al. (1963), equations can be written for the genetic loads at a single diallelic autosomal locus, assuming a homozygous selection coefficient s and dominance coefficient h.The fitness of mutant homozygotes relative to wild-type is 1 -s and the fitness of heterozygotes is 1 -hs; s may vary across loci, but h is treated as a constant, although it is easy to relax this assumption.L i is the genetic load for individuals homozygous for karyotype i produced by random mating within the population, defined as the expected reduction below 1 of their mean fitness relative to wild-type homozygotes.L 12 is the corresponding load for heterokaryotypes.The homozygous load H i is the reduction below 1 in the expected relative fitness of individuals made homozygous for gametes with karyotype i, with probability 〈q i 〉 of being homozygous for the mutant allele.B i is the inbreeding load for karyotype i, defined as H i -L i (Charlesworth and Charlesworth 2010, p. 173).
Following Charlesworth (2018), some simple algebra yields the following expressions for these quantities: Inversion polymorphisms and fitness | 3 where The selective difference between heterokaryotypes and a homokaryotype of class i is measured by t i = L i -L 12 .Heterozygote advantage exists if both t i s are positive.To obtain conditions for heterozygote advantage, it is useful to rearrange Equations (1) by writing 〈q 1 〉 = 〈q〉 + 〈δq〉, 〈q 2 〉 = 〈q〉 − 〈δq〉, where 2〈δq〉 is the expected difference in the frequency of A 2 between In and St. 〈δq〉 will be nonnegative if x < ½, due to the greater effectiveness of drift relative to selection against nonrecessive, mildly deleterious mutations in a smaller population (Kimura et al 1963).Substituting these expressions for the q i into Equations ( 1a) and (1b), the following expressions are obtained: Equation (2c) shows that when h < ½ the mean fitness of heterokaryotypes is increased by a difference in the expected frequencies of deleterious mutations between the two karyotypes, as expected intuitively.These expressions yield the following results for the selective differences between heterokaryotypes and homokaryotypes: It is easily seen that, for an inversion with x < ½ and 〈δq〉 > 0 (see above), we have t 1 > 0 provided that h ≤ ½ and C 12 ≤ 0, so that In/In then has a lower fitness than In/St.The expectation for a lowfrequency inversion is thus that t 1 > 0, although the magnitude of the effect is likely to be small for a single, large population.Given that F 1 > 0 and h < ½, this is also the case even when 〈δq〉 = 0, reflecting the fact that genetic drift causes a reduction in mean fitness by increasing the frequencies of homozygotes for recessive or partially recessive mutations; this effect is not experienced by In/St individuals unless C 12 > 0, which can be ruled out by the argument given above.
If 〈δq〉 > 0 and h ≤ ½, the condition for t 2 > 0 is more stringent than for t 1 > 0, due to the opposite sign of the term in 〈δq〉 in Equations ( 3a) and (3b), especially if F 2 is close to 0. It may therefore be difficult to find conditions in which there is an advantage to In/St over both homokaryotypes, unless the inversion frequency is close to ½, so that 〈δq〉 ≈ 0.
These results can be generalized to the case of a subdivided population with a constant frequency of the inversion across all local populations (demes), by taking expectations of within-deme allele frequencies across populations, as described in Section 5 of Supplementary File 1.

A single population: modeling drift and selection
In order to obtain numerical results for the load statistics described above, expressions for the means and variances of the q i , as well their covariance, are needed.Recombination is assumed to be absent in heterokaryotypes.We first consider the expected changes in allele frequencies due to selection within each karyotype.For In karyotypes, the marginal fitness of haplotypes carrying the wild-type allele (A 1 ) at a locus is easily seen to be Similarly, the marginal fitness of In haplotypes carrying the mutant allele is The net expected change in the frequency of A 2 within In karyotypes due to selection (neglecting second-order terms in s) is thus Similarly, the net expected change in the frequency of A 2 within St karyotypes due to selection is Equations ( 4) and (5) bring out the interdependence between the evolutionary processes in the two karyotypes when h ≠ ½.To proceed further, the effects of mutation and drift also need to be analyzed, such that an expression for the stationary joint probability density function (p.d.f.) for q 1 and q 2 , ϕ(q 1 , q 2 ), can be obtained.Kimura (1964, p. 41) derived a pair of coupled forward diffusion equations describing the joint stationary p.d.f. for two variables, using the principle that a zero flux of the probability density of each variable across all their values guarantees a stationary joint distribution.This method can be applied to the above selection equations, together with the terms arising from mutation.
Let the rates of mutation from A 1 to A 2 and vice versa be u and v, respectively, with scaled mutation rates The mutational bias toward deleterious mutations, κ, is equal to u/v.As shown by Kimura (1964), in order to analyze this type of situation it is most convenient to use the natural logarithm of ϕ , ψ = ln ϕ.If drift affects q 1 and q 2 independently, as assumed here, the conditions for a stationary distribution are where V δq1 = p 1 q 1 /(2N 1 ) and V δq2 = p 2 q 2 /(2N 2 ) are the variances in the changes in allele frequencies per generation due to drift within In and St, respectively.Correspondingly, ∂ q1 V δq1 = (1-2q 1 )/2Nx; . The Δq i are given by the selection equations ( 4) and ( 5) together with the relevant mutation terms.For a meaningful solution of Equation ( 6) to exist, (Kimura 1964, p. 41).Kimura showed that this is the case for the mutation terms, so we need only consider the selection terms contributed by Equations ( 4) and ( 5).The following expression satisfies both this condition and Equation ( 6) for the selection contribution to ψ: where γ is the scaled selection coefficient, 2Ns.
The full solution to Equation ( 6) is thus where C is a normalization constant, given by the inverse of the double integral of ϕ over the closed intervals (0, 1) for q 1 and q 2 .

A single population: obtaining the mutational load and population genomic statistics
To obtain the single-locus load statistics for a given h, numerical integration of Equation ( 8) and its product with powers and crossproducts of the q i can be performed, for the purpose of determining the expectations and variances of the q i and their covariance C 12 .The corresponding F i statistics can be obtained from Equation (1e).The means of the single-locus load statistics, with the terms in s omitted, can then be obtained using Equations (1).Details of the integration procedures are given in Section 1 of Supplementary File 1.
In order to calculate the load statistics themselves under reasonably realistic assumptions, a gamma distribution of the scaled selection coefficient γ = 2Ns is assumed, with a p.d.f.given by where a is the shape parameter, b = ̅ γ/a is the location parameter, and Γ(a) is the gamma function.This distribution has been widely used in population genomic methods for estimating the distribution of fitness effects of deleterious mutations (e.g.Charlesworth 2015;Booker et al. 2017).
The values of ̅ γ and the shape parameter a are chosen to correspond to estimates from the population genomics studies just mentioned, which indicate ̅ γ values of hundreds or thousands for nonsynonymous mutations and shape parameters of approximately 0.3, implying a wide distribution of selection coefficients.It is assumed that fitness effects are multiplicative across sites, so that the products of the number of sites, n s , and the expectations of Equations ( 1a)-(1c) over the distribution of γ correspond to the natural logarithms of the corresponding multisite load statistics.The exponentials of the negatives of these expressions then yield the mean fitnesses of the karyotypes concerned, relative to that of mutant-free individuals.The exponential of the negative of the product of n s and the expectation of Equation (1c) yields the fitness of totally inbred individuals of karyotype i relative to outbred individuals that are homozygous for karyotype i, i.e. a measure of the inbreeding depression experienced by carriers of karyotype i.
The net selection coefficient against homozygosity for karyotype i relative to heterokaryotypes for a given n s is measured by where E{} indicates the expectation over the distribution of γ (as opposed to an expectation over the distribution of q, denoted by angle brackets).
It is also of interest to summarize the expected patterns of variation at the loci themselves.To do this, the p.d.f. for q 1 and q 2 are used to calculate the expected nucleotide site diversities within each karyotype subpopulation for a given selection coefficient, π i = 2〈p i 〉〈q i 〉(1 − F i ), which are then averaged over the distribution of s.In addition, the expected proportion of segregating sites for a sample of n alleles, S ni , can be determined as described in Section 2 of Supplementary File 1. Division by the sum of the harmonic series, a n =  j=n−1 j=1 (1/j), yields the expected values of Watterson's theta (θ wni ) for each subpopulation (Watterson 1975).The skew of the distribution of segregating variants toward rare variants for subpopulation i can conveniently be measured by Δθ wni = 1 − π i /θ wi (Campos and Charlesworth 2019).
For practical purposes of computation, it is convenient to divide the range of values of γ into several zones according to the strength of selection and to compute the integrals of the load statistics over each zone separately.The overall values of the load statistics are then given by summing the results over all zones.The details are given in Section 1 of the Appendix, and the computer code for generating the numerical results for this case is given in Supplementary File 2.

A finite island model metapopulation
In this case, a metapopulation of total size N T is divided into a large number d of subpopulations (demes), each of size N = N T /d.A Wright-Fisher model is assumed to apply to each deme.N is assumed to be sufficiently large that the frequency of the inversion is held at the same frequency x in all demes.A fraction m of each deme is derived by migration from a pool with equal contributions from all demes.Let the current mean frequency across all demes of the mutant allele A 2 at a locus be ̅ q i for karyotype i, so that migrants contribute m ̅ q i to the new frequency of A 2 among type i haplotypes within a deme and m ̅ p i to the new frequency of A 1 .This model poses the problem that the evolutionary processes within demes change the values of the ̅ q i , so that they cannot realistically be treated as fixed quantities.Following previous treatments of this problem, it is assumed here that the process of change in the ̅ q i can be described by a pair of coupled diffusion equations, using the expectations of p i q i /(2N T ) over all d demes as the drift variance terms together with the corresponding expectations of the expressions for the deterministic changes in ̅ q i (Whitlock 2002;Cherry and Wakeley 2003;Roze and Rousset 2003;Wakeley 2003).The mutational contributions to the latter are simply obtained by substituting ̅ q i for q i in the standard formulae for the deterministic changes in the ̅ q i , since the mutational changes are linear in q i .The expectation of p i q i over demes conditioned on ̅ q i can be written as , where F STi is the variance among demes in q i divided by ̅ p i ̅ q i .The total effective population size for the metapopulation for karyotype i, N emi , is thus equal to N T /(1 -F STi ) (Wakeley and Aliacar 2001).
The nonlinearity with respect to q i of the selection terms for allele frequency changes within demes (after division by p i q i ) when h ≠ 0.5 means that an exact closed expression for their contributions to the expected changes in the ̅ q i cannot be obtained, except in the absence of dominance (Whitlock 2002;Roze and Rousset 2003;Wakeley 2003).However, a useful approximation can be obtained by neglecting the 3rd moments about their means of the within-deme allele frequencies; these moments are necessarily smaller than the variances of the ̅ q i and will be considerably smaller when selection is sufficiently strong in relation to drift that the ̅ q i are close to 0. When selection is sufficiently weak, a neutral approximation for the 3rd moment can be used (Whitlock 2002), but this will in general overcorrect when selection is strong and so is Inversion polymorphisms and fitness | 5 not employed here.The error introduced by ignoring this correction affects only the small portion of the distribution of selection coefficients where s is O(1/2N T ) and is thus unlikely to be important for the load and population genomic statistics calculated here.
The following expressions are obtained after some algebra: We also have the following expression for the variance in ̅ q i due to drift, e.g.Whitlock (2002): where N T1 = N T x and N T2 = N T y.
Combining Equations ( 11) and ( 12), and using the approach that led to Equation ( 8) and carrying out some rearrangements of terms, we obtain the equivalent of the ψ function for the panmictic case, which describes the selection contribution to the logarithm of the p.d.f. for the metapopulation.The details are given in Section 2 of the Appendix.Section 5 of Supplementary File 1 describes how the p.d.f.can be used to obtain the load statistics.The computer code for generating the numerical results is given in Supplementary File 3.

General considerations
Intuitively, the mutational load associated with each homokaryotype in an inversion polymorphism would be expected to be strongly affected by the dominance coefficient (h), the scaled strength of selection (2Ns = γ), the mutation rate to deleterious mutations (u), the number of selected sites that it contains (n s ), and its frequency (x).Dominance coefficients less than one-half are well known to be required for inbreeding depression and heterosis, whose magnitudes are inversely related to h (Charlesworth and Charlesworth 2010, Chap. 4).It would thus be expected that the reduction in fitness associated with an arrangement, and any fitness advantage to In/St heterokaryotypes, would decrease with h but increase with n s and u, if γ is sufficiently small that deleterious mutations are significantly affected by drift.
It is less clear how these fitness effects are related to γ, since stronger selection reduces the frequencies of deleterious mutations but also reduces their effects on fitness if they rise to high frequencies.However, if selection is so strong in relation to drift that allele frequencies are at mutation-selection equilibrium, no differentiation in allele frequencies between In and St will occur, removing any possibility of a selective advantage to heterokaryotypes (Sturtevant and Mather 1938;Nei et al. 1967).The mean fitnesses of both homokaryotypes with n s selected sites will then be equal to the deterministic value, exp(−2n s u) ≈ 1 −2n s u when fitnesses are multiplicative and h is > 0 (Haldane 1937).
The rarer of the two arrangements experiences more genetic drift than its counterpart, so that x < 0.5 means that inversion homokaryotypes should have a lower overall mean fitness than standard homokaryotypes.It is less clear when an advantage to heterokaryotypes can be generated; the increased load associated with the rarer arrangement may simply generate a net selective advantage to its counterpart, especially when mutations are only partially recessive.Finally, it is likely that population subdivision will increase the magnitude of the mutational loads for each homokaryotype and the selective differences among the three karyotypes, because drift within demes in a metapopulation occurs at a faster rate than in a single population with the same size as the metapopulation.However, this is counterbalanced by a slower rate of drift for the metapopulation as a whole (Whitlock 2002;Cherry and Wakeley 2003;Wakeley 2003), so the net effect is hard to predict intuitively.

A single population: load statistics
Numerical results for a single randomly mating population are presented here.These are based on Equations (1) for individual selected sites together with the procedures for combining the effects of mutation, drift, and selection for all sites that were described above.The number of sites was set to 10 5 , corresponding to an inversion containing 100 genes with a mean of 1,000 nonsynonymous sites per gene.The expectations of the single-locus mutational loads (L i and L 12 ) and the inbreeding loads (B i ), defined by Equations ( 1), are then multiplied by 10 5 to obtain their net values.If multiplicative fitnesses are assumed, these quantities are equivalent to the negatives of the natural logarithms of the corresponding mean fitnesses (for the Ls) or differences in log mean fitnesses (for the Bs).The values of the corresponding mean fitnesses for a different number of sites, n s1 , can be found by taking the exponentials of their negatives multiplied by n s1 /10 5 .The selection coefficients t i against the 2 homokaryotypes are obtained by exponentiation of the product of n s1 and the expectation of L i -L 12 (Equation 10).If these products are small, as is mostly the case in practice, t i ≈ n s E{L i -L 12 }.Heterokaryotype advantage requires both of the t i to be positive; if one is positive and the other negative, there is directional selection in favor of the karyotype with the negative t i .
Figure 1 displays the values of the relevant load statistics and the selection coefficients against homokaryotypes as a function of the dominance coefficient h, for three different inversion frequencies and two different mean scaled selection coefficients (̅ γ = 1,000 and 4,000), using a gamma distribution of selection coefficients with shape parameter a = 0.3, and a mutation rate to deleterious alleles of u = 5 × 10 −9 .These values are broadly consistent with population genomic estimates from Drosophila melanogaster (Charlesworth 2015).A strong mutational bias to deleterious mutations of κ = 1.5 was assumed, in order to maximize the magnitude of the mutational loads and selection coefficients.
If there were no effects of drift, the net load (L) for both subpopulations would be approximately 2n s u = 0.001, and the net inbreeding load (B) would equal n s u(h −1 -2)-see Charlesworth and Charlesworth (2010, Chap. 4).This gives B = 0.002 with h = 0.25 and 0.036 with h = 0.05, with zero selection coefficients for both homokaryotypes.The numerical results show that the loads for both In and St subpopulations are always substantially larger than the deterministic values, even with h = ½ (Supplementary Table 1).With x = 0.1 and ̅ γ = 1,000, L 1 , the load for In is more than 20-fold greater than the deterministic value for h between 0.05 and 0.5.Both L values are higher for the lower ̅ γ values, as might be expected from the larger effects of drift in this case.Similarly, with x < ½, L 1 is always larger than L 2 , and decreases with x for a given h and ̅ γ .Furthermore, L 1 strongly decreases with h for x = 0.1 and 0.3, whereas L 2 is only slightly affected by h.The fact that L 1 < L 2 with h = ½ and x < ½ shows that the increase in load due to smaller population size is partly caused by increased mean frequencies of deleterious mutations, not simply by increased frequencies of homozygotes.
In contrast, the B i values are always lower than the deterministic values, reflecting the reduction in variability caused by drift, and are decreasing functions of h (vanishing when h = 0.5), although this is hard to see from Fig. 1 due to their overall very small values.These effects are especially noticeable for B 1 when x = 0.1 and ̅ γ = 1,000.The small B i values reflect the fact that, except for the lowest dominance coefficient (h = 0.05), the L i are always quite close to the corresponding homozygous loads (the H i of Equation 1c).For example, with ̅ γ = 1,000 and x = 0.1, for h = 0.05, we have There can be substantial selection against the inversion homokaryotypes at the lower h values, especially for the smaller x values: t 1 reaches 0.19 with h = 0.05 and ̅ γ = 1,000, but decreases sharply with x and h; it is only 0.04 for h = 0.25, x = 0.1, and ̅ γ = 1,000.The values for ̅ γ = 4,000 are somewhat smaller, consistent with smaller effects of drift in causing divergence between the 2 subpopulations.In contrast, t 2 is mostly negative and quite small (approximately −0.01 for a wide range of h values when x = 0.1 and ̅ γ = 1,000), indicating weak directional selection against the inversion.Only when x is close to 0.5 does the heterokaryotype experience a slight advantage over both homokaryotypes (a maximal value of t 1 = t 2 ≈ 0.004 for h = 0.05 and ̅ γ = 1,000).With x = 0.5, the selective advantage to heterokaryotypes declines sharply with h, and is only 0.0016 for h = 0.25 and ̅ γ = 1,000; it vanishes when h = 0.5.With 10 6 instead of 10 5 selected sites, the heterokaryotype advantage would be approximately 10 times as large.With x = 0.5, symmetry implies that t 1 = t 2 , and the equality of effective population sizes means that the mean fitness of the heterokaryotypes is superior to that of the homokaryotypes purely because they have a lower frequency of mutant homozygotes than either of the homokaryotypic subpopulations.
It is also of interest to examine the effect of differences in population size and selection coefficients on the load statistics for a constant scaled selection strength.Figure 2 is similar to Fig. 1, but with population sizes of 2 × 10 6 and 10 5 , and a mean selection coefficient (0.001) that is one-half of that in Fig. 1.A comparison of the top panels of Figs. 1 and 2 shows that reducing the mean selection coefficient by one-half, but keeping ̅ γ constant, results in a reduction in the L values, by a factor of close to 2 for the case of the In subpopulation with x = 0.1, but by considerably less for the St population, where L 2 is only slightly greater than the deterministic value of 0.001, even for h = 0.05.The relative effect on a subpopulation is reduced as it increases in size, as might be expected intuitively.These effects reflect the fact that the load caused by drift is greater when the selection coefficients involved are larger (for fixed ̅ γ ), in contrast to the deterministic formula L = 2u.The two selection coefficients on homokaryotypes relative to heterokaryotypes, t i , show a similar pattern.In contrast, the inbreeding loads are only slightly reduced by a reduction in the strength of   The effect of a reduced population size while holding mean s constant can be seen by comparing the upper and lower panels of Fig. 2. As might be expected from the greater effects of drift in a smaller population, the L i and t i are considerable larger when the population size is reduced by a factor of 20, especially for the inversion subpopulation with x = 0.1, for which N 1 = 10,000 in the lower panel.However, with x = 0.5 and h = 0.25, the selective advantage to the heterokaryotypes is only 0.0003 for the smaller population size.The inbreeding loads are barely affected by the population size difference.

A single population: population genomics statistics
It is also of interest to examine the effects of an inversion polymorphism on population genomics statistics that can be used to assess the effects of the differences in effective population size between the In and St subpopulations, and between these and the part of the genome that is independent of the inversion polymorphism.Figure 3 shows the results for the St population using the same parameter values as in Fig. 1 for 3 different values of x; the case with x = 0.1 is nearly equivalent to the situation for the rest of the genome.It can be seen, somewhat surprisingly, that the mean frequency of deleterious alleles in the St subpopulation (̅ q 2 , half-filled triangles), which includes both fixed and segregating sites, is nearly independent of h.It increases slightly as the size of the St population, given by N 2 = N(1 -x), decreases as x changes from 0.1 to 0.5; the effect of N 2 is somewhat greater for ̅ γ = 1,000, as would be expected from the greater effects of drift in this case.This difference partly reflects the higher proportion of sites that can become fixed for deleterious mutations with weaker selection, as well as their higher mean frequency at segregating sites.Its behavior implies that the effect of population size on the load for a subpopulation when the mean strength of selection is strong is largely caused by changes in the variance of allele frequency rather than its mean.
The mean nucleotide site diversity at the selected sites for St (π sel2 , half-filled lozenges) decreases somewhat with increasing h and decreasing N 2 , but is always between 0.0011 and 0.0025 for ̅ γ = 4,000 and between 0.0017 and 0.0037 for ̅ γ = 1,000.The ratio of π sel2 to the mean nucleotide site diversity at neutral sites (π sel2 /π neut2 , filled triangles) decreases with N 2 , especially when h is small.The measure of skew toward low-frequency variants at selected sites, Δθ w2 , for a sample of 20 haploid genomes (filled lozenges) increases with N 2 , but the effect is weak unless h is small; the effect of h on Δθ w2 is smaller than on π sel2 /π neut2 , and plateaus around h = 0.25 (for the definition of Δθ w , see A single population: modeling drift and selection).
These effects of h and N 2 can also be seen in the right-hand panels of Fig. 4 (labeled as x = 0.5 vs x = 0.1) which plot the ratios of ̅ q 2 , π sel2 , π sel2 /π neut2 , and Δθ w2 with x = 0.5 to their values with x = 0.1, a 1.8-fold difference in N 2 .The plots for π sel2 /π neut2 bring out clearly that a lower subpopulation size is associated with larger π sel2 /π neut2 , especially when h is small.For Δθ w2 , the ratio is approximately 0.85 for h = 0.05, but increases rapidly toward 1 as h increases.The 4-fold difference in ̅ γ between the upper and lower panels has a remarkably small effect on the ratios for all 4 statistics.
The other 2 panels of Fig. 4 display plots against h of the ratios of these statistics for the In subpopulation to their values for the St subpopulation for x = 0.1 and 0.3.The ratios of population sizes for St vs In are 9-fold for x = 0.1 and 2.3-fold for x = 0.3, flanking the ratio N 2 /N 1 for the right-hand panels.N 2 for x = 0.3 is 1.4-fold greater than N 2 for x = 0.5, which enhances the contrast between In and St. Accordingly, the patterns are more marked than for the right-hand-most panels, and are strongest for the case with x = 0.1, with a maximal ratio of π sel1 /π neut1 to π sel2 /π neut2 of 2.4 at h = 0.05 and ̅ γ = 4,000.This is still much less than the ratio of 9 for neutral diversity.The ratio Δθ w1 /Δθ w2 with x = 0.1 is close to 0.5 for both ̅ γ values when h = 0.05, but rapidly becomes close to 1 as h increases-for x = 0.1, h = 0.25, and ̅ γ = 4,000, we have Δθ w1 /Δθ w2 = 0.915.The skew in the site frequency spectrum at selected sites is thus unlikely to be a powerful statistic for detecting a reduced efficacy of selection on a low-frequency inversion, given that very small h values are unlikely for mildly deleterious mutations (Crow 1993;Manna et al. 2011).

A subdivided population
As described in General considerations, A finite island model metapopulation, this case assumes that a metapopulation of total size N T is divided into a large number d of demes, each of size N = N T /d.A Wright-Fisher model applies to each deme, and the deme size is assumed to be sufficiently large that the frequency of the inversion is held at the same frequency x in all demes.A fraction m of each deme is derived by migration from a pool with equal contributions from all demes.For sites that are independent of the inversion, the level of neutral genetic differentiation between demes is described by F STn ≈ 1/(1 + M), where M is the scaled mutation rate 4Nm.Different levels of subdivision are characterized by different values of F STn ; for simplicity, the subscript n is dropped in what follows.The scaled selection parameter γ is now defined as 2N T s, and the scaled mutation parameters are α = 4N T u and β = 4N T v.

Inversion polymorphisms and fitness | 9
Figure 5 shows the effects of population subdivision on the load statistics of most interest, together with the ratio of the diversities at selected sites for the In vs St subpopulations.There is an inversion frequency of 0.1 in a metapopulation of 200 demes with a total size N T = 10 6 (N = 5,000); the selection and mutation parameters are the same as in Fig. 1.The results for F ST = 0 were obtained from the single population calculations described above.The results can be summarized very simply: there is a remarkably small effect of population subdivision as F ST changes from 0 to 0.25, with the most marked effect occurring over the change from F ST = 0 to F ST = 0.05, especially for the smallest dominance coefficient (h = 0.05).For example, with ̅ γ = 1,000 and h = 0.05, the selection coefficients against the In and St homokaryotypes relative to the heterokaryotype change from t 1 = 0.1868 and t 2 −0.0107 with F ST = 0 to 0.211 and −0.0120 with F ST = 0.05, reaching 0.250 and −0.0142 at F ST = 0.25.With the more plausible value of h = 0.25, the changes are much smaller: t 1 = 0.0405 and t 2 −0.0136 at F ST = 0, and 0.0442 and −0.0146 at F ST = 0.25.Both the absolute values of the loads and selection coefficients and their dependence on F ST are much smaller with ̅ γ = 4,000 than 1,000.As in the case of a single population, a selective advantage to the heterokaryotype is not found unless there are nearly equal frequencies of In and St (Supplementary Table 3).The magnitude of such an advantage is not greatly increased by subdivision; for example, with x = 0.5, ̅ γ = 1,000 and h = 0.05, t 1 = t 2 = 0.0039 at F ST = 0.0, and 0.00860 at F ST = 0.25; with h = 0.25, the corresponding values are 0.00160 and 0.00196.
The relation with F ST of the ratio of diversities at selected sites for In vs St with x = 0.1 is also shown in Fig. 5, and is similarly rather weak.This finding also applies to the population genomic statistics that are not displayed in Fig. 5-see Supplementary Table 3 for details (the Δθ w statistic was not calculated, as this statistic is hard to evaluate with population subdivision and was not very informative in the single population case).For example, with x = 0.1, ̅ γ = 1,000 and h = 0.05, the ratio of π sel /π neut for In vs St decreases from 2.32 with F ST = 0 to 1.80 with F ST = 0.25, indicating a greater efficacy of selection on the smaller subpopulation when there is greater subdivision.With h = 0.25, the change is much smaller, from 1.97 to 1.90.Similarly, the ratio of the mean frequency of A 2 for In vs St with h = 0.05 changes from 3.02 with F ST = 0 to 3.78 with F ST = 0.25, but only from 2.21 to 2.32 with h = 0.25.
These results are all for a relatively large deme size of N = 5,000.Intuitively, it might seem that reducing the deme size would enhance the effects of drift within demes, and lead to larger loads and magnitudes of the selection coefficients, as well as reducing the values of such signatures of purifying selection as the mean frequencies of mutant alleles and the ratios of diversities at selected sites to neutral sites.The example in the upper part of Table 2, where deme sizes of 500 and 5,000 are compared for the same neutral F ST value, shows that this expectation is met for Fig. 4. The first 2 panels in each row show the ratios of the population genomics statistics for the In subpopulation to those for the St subpopulation, plotted against the dominance coefficient h.The third panel shows the ratios of these statistics for the St subpopulation with an inversion frequency of 0.5 to their values for a frequency of 0.1, for which the corresponding ratio of population sizes is 0.556.The evolutionary parameters are the same as in Fig. 1.The half-filled triangles are the ratios of the mean frequencies of A 2 ; the half-filled lozenges are the ratios of the mean diversities at selected sites.The filled triangles are the ratios of mean diversities at selected sites to mean diversities at neutral sites.The filled lozenges are the ratios of Δθ w values for the selected sites.
the case of equal frequencies of In and St.For the mean fitnesses of the homokaryotypes and the selection coefficients, the effects are small and are only visible in the table in a few cases.However, the mean frequencies of mutations and the ratios of diversity at selected sites vs neutral sites show a clear pattern of reduced efficacy of selection with the smaller deme size, which is reduced in magnitude by larger h and F ST .Further results for cases with smaller deme sizes are shown in Supplementary Table 4.
The patterns are somewhat more complex, however, when there is a large difference in frequency between arrangements, as shown in Table 2 for the case of x = 0.1.Here, the differences in the load statistics between the large and small deme size cases are negligible.There is, however, a signature of a slightly enhanced efficacy of selection with the smaller deme size in the In subpopulation relative to St, indicated by a consistently smaller value of ̅ q 1 /̅ q 2 when N = 500.Confusingly, π sel1 /π neut1 is greater than π sel2 /π neut2 when N = 500, indicating the opposite pattern.The first of these results is explained by the fact that the much smaller size of the In subpopulation means that mutations are behaving nearly neutrally within demes, so that a lower deme size will not have much of an effect, whereas it will reduce the overall efficacy of selection on the St subpopulation, resulting in an increase in ̅ q 2 .For example, with h = 0.05 and F ST = 0.05, the ratio of ̅ q 2 for N = 500 vs ̅ q 2 for N = 5,000 is 1.024, whereas the value for ̅ q 1 is 1.006.The corresponding ratio for π sel1 /π neut1 is 1.072 whereas that for π sel2 /π neut2 is 1.005.This is entirely due to a higher ratio of π sel1 for N = 500 vs N = 5,000, since there is no effect of local deme size on neutral diversity for a fixed neutral F ST .It is not entirely clear how to interpret this pattern, but one factor is likely to be the fact that the relaxation of the efficacy of selection in St subpopulation leads to an increase in F ST at selected sites, which would work against any increase in π sel2 due to a reduced efficacy of selection.

Discussion
It should be borne in mind that the results described above relate to a situation in which an inversion polymorphism is maintained by balancing selection that is invariant over space and sufficiently strong that the frequency of the arrangements is constant over time and space.This obviously does not apply to many natural situations, e.g. when there are clinal patterns of variation in inversion frequencies or temporal fluctuations in frequencies, as is often the case (Krimbas and Powell 1992;Kapun and Flatt 2019).Nevertheless, the theoretical results described above address several questions that can in principle be answered by comparisons with empirical studies.First, does a low-frequency arrangement accumulate a larger mutational load than its counterpart?Fig. 5.The mutational load and polymorphism statistics for an inversion polymorphism with inversion frequency x = 0.1 in a subdivided island population of total size N T = 10 6 with 200 demes and 10 5 selected sites, plotted against F ST for neutral sites independent of the inversion.The results for 3 different values of the dominance coefficient, h, are shown.The selection coefficients follow a gamma distribution with a mean of 5 × 10 −4 and a shape parameter of 0.3.The mutation rate to deleterious alleles is 5 × 10 -9 per basepair, with a mutational bias toward deleterious variants of 1.5.The upper and lower panels have products of 2N T and mean s of 4,000 and 1,000, respectively.The filled and half-filled symbols denote values for the In and St subpopulations, respectively.The lozenges are the net mutational loads within the respective karyotypes, and the triangles are the selection coefficients against homokaryotypes relative to the heterokaryotype.The crosses are the ratios of the nucleotide site diversities for In vsSt subpopulations.

Inversion polymorphisms and fitness | 11
Second, can mutational load contribute a significant selective advantage to heterokaryotypes, which might help to stabilize the polymorphism?Third, can differences in population genomic statistics between common and rare arrangements shed light on differences in mutational load between arrangements?These questions are each discussed in turn below.The focus will be on the results for a single population, since the analysis of the properties of a subdivided population showed that there were no major differences from those of a single population over the range of F ST values examined, although subdivision tends to slightly magnify the differences between rare and common arrangements.

Are the relative values of the fitness components of carriers of different gene arrangements consistent with mutational load?
The relevant variables with respect to the relative fitnesses of the In and St subpopulations are the L i , the mutational loads associated with individuals produced by random mating among outbred individuals homozygous for karyotype i, where i = 1 for the rarer arrangement (In) and i = 2 for the more common one (St).These quantities are, course not directly observable, since the fitness of mutation-free individuals is unknown, but they can be used to predict the relative values of the mean fitnesses (or fitness components) for In and St. On the assumption of multiplicative fitnesses, the ratio of the mean fitnesses of type 1 and type 2 homokaryotypes is ̅ w 1 / ̅ w 2 = exp(L 2 − L 1 ), corresponding to a difference L 2 -L 1 in the natural logarithm of mean fitness.The mutational and selection parameters used here (a mutation rate of 5 × 10 -9 per site and a gamma distribution of selection coefficients with shape parameter 0.3) are consistent with the results of population genomic studies of D. melanogaster, e.g.Kousathanas and Keightley (2013) and Wang et al. (2023).
In order to relate the theoretical predictions to data on inversion polymorphisms, it is necessary to have information about h and the parameters of the distribution of s, which is made difficult by the fact that there is much uncertainty about the relation between h and s.The evidence from Drosophila studies of the effects of mutations on fitness components suggest that only very strongly selected deleterious mutations such as homozygous lethals have h values low as 0.05, whereas the much more abundant deleterious mutations with s ≤ 0.02 have h values of the order of 0.25 or more (Crow 1993;Manna et al. 2011).It therefore seems safe to use h = 0.25 as a working value for comparing theory with data, since mutations like lethals that are strongly selected against when heterozygous will be held close to their deterministic equilibrium values in both the In and St subpopulations (Nei 1968), and hence will not contribute to the genetic differences between them.It is worth noting that a higher load for the rarer arrangement does not require mutations to have h < ½, contrary to what is often stated (e.g.Jay et al. 2021).If drift is sufficiently strong in relation to selection that there is a higher expected frequency of deleterious mutations in the smaller subpopulation, a higher load can arise even for completely dominant mutations.
Estimates of the distribution of γ = 2N e s from population genomics studies reflect the more weakly selected part of the distribution of selection coefficients and are thus the most useful source of information for the present purpose.With a shape parameter of 0.3 and h = 0.25, π sel2 /π neut2 = 0.15 with ̅ γ = 4,000 when x = 0.1.This is only slightly higher than the observed ratios of nonsynonymous to silent site diversities in normally recombining regions of the genome in putatively ancestral range populations of D. melanogaster (e.g.Campos et al. 2014), suggesting that this value of ̅ γ can be used as a working estimate for nonsynonymous mutations.
One of the best-studied D. melanogaster inversions is In(3R)P.This has a frequency of 0.1 in populations in its ancestral range in Africa (Kapun and Flatt 2019), which are the most relevant populations for comparisons with the theoretical predictions.It is in the middle of the size range for polymorphic inversions in this species, covering approximately 8 Mb of sequence (Kapun et al. 2023).This corresponds to approximately 1,000 protein-coding sequences, i.e. 10 6 nonsynonymous sites rather than the 10 5 selected sites illustrated in the figures.Functional noncoding sequences may also contribute to the mutational load, and these appear to be under weaker selective constraints than nonsynonymous mutations (Andolfatto 2005;Casillas et al. 2007;Campos et al. 2017).As a rough estimate of their contribution to the load and population genomic statistics, it is plausible to assume that functional noncoding sites are three times as abundant as nonsynonymous sites (Halligan and Keightley 2006) but have ̅ γ = 1,000 rather than 4,000.With x = 0.1 and h = 0.25, the results in Supplementary Table 1 imply that, by combining the effects of 10 6 nonsynonymous sites with ̅ γ = 4,000 and 3 × 10 6 noncoding sites with ̅ γ = 1,000 and h = 0.25, we would have L 1 = 2.09 and L 2 = 0.08 for an inversion with frequency 0.1, giving a predicted ̅ w 1 / ̅ w 2 = 0.13.Predictions of this kind are highly sensitive to the frequency of the inversion, and to assumptions concerning the abundance of functional sites and the strength and mode of selection.With the model just described and an inversion frequency of 0.3, L 1 = 0.42 and L 2 = 0.11, giving ̅ w 1 / ̅ w 2 = 0.73.If only nonsynonymous sites are considered, ̅ w 1 / ̅ w 2 rises to 0.69 for x = 0.1 and 0.95 for x = 0.3.Synergistic epistasis among deleterious mutations Table 2. Some load and population genomics statistics for a subdivided population of total size 10 6 for 2 different local deme sizes.Results for N = 500 and 5,000 are shown in the upper and lower parts of each cell, respectively.̅ w1 and ̅ w2 are the mean fitnesses of the products of random mating among In and St karyotypes, respectively.The selection and mutation parameters are the same as in Fig. 1. is expected to reduce the L i by a factor of approximately 2 (Kondrashov 1995;Charlesworth 2013), so that the model of a mixture of 10 6 nonsynonymous and 3 × 10 6 functional noncoding sites would generate ̅ w 1 / ̅ w 2 values of 0.37 and 0.86 for x = 0.1 and x = 0.3, respectively.These complexities means that it is difficult to make rigorous comparisons between theory and data, and the estimate of ̅ w 1 / ̅ w 2 = 0.13 for In(3R)P is probably at the lower end of what is to be expected.
A large number of studies of the effects of inversions on fitness and components of fitness under laboratory conditions have been published, with D. melanogaster and D. pseudoobscura being been the most intensively studied species-see reviews by Sperlich and Pfriem (1986), Krimbas and Powell (1992), and Kapun and Flatt (2019).Many of these studies do not, however, provide useful information for the present purpose, either because they involve traits like male mating success that are hard to quantify in terms of relative fitness measures, because only small numbers of independently sampled arrangement haplotypes were used, or because samples from laboratory rather than natural populations were used.The most informative available estimates are shown in Supplementary File 4, with their standard errors where available.
For the two large-scale studies of egg-to-adult viability in D. melanogaster (Mukai and Yamaguchi 1974;Watanabe et al. 1976), where balancer crosses were used to extract 2nd or 3rd chromosomes from wild flies, the crosses involving pairs of independently extracted chromosomes were divided into cases where either each member of a pair was inversion-free or at least one member of each pair carried an inverted chromosome.There is no evidence for a difference between these categories either study.One possible explanation for this discrepancy is the fact that a single component of fitness such as viability reflects only a portion of the net effect of deleterious mutations on fitness.Table S8 in Charlesworth (2015) presents estimates of the ratios of the fitness effects of deleterious mutation on various fitness components to their effects on net fitness (the α parameters).The estimates of α are subject to considerable uncertainty, but even the lowest estimate for viability (0.09) cannot explain the absence of effects on viability in these experiments.
Studies such as these that involve populations of D. melanogaster outside the ancestral range of the species in south-eastern Africa are difficult to interpret, because of the effects on variability of the severe population bottlenecks associated with the spread of the species out of Africa (e.g.Haddrill et al. 2005).Indeed, a curious feature of the data of Mukai and Yamaguchi (1974) on the 2nd chromosome is the much higher frequency of homozygous lethal In chromosomes than St chromosomes (54 vs 36%); this is associated with a higher frequency of pairs of crosses in which combinations of different lethal-bearing chromosomes were lethal (19 vs 11%), suggesting that there may have been a recent population size bottleneck that affected the In chromosomes more severely than the more abundant St chromosomes.Another possibility, proposed by a reviewer, is that high frequency lethals associated with the inversions were also associated with the Sd segregation distorter haplotype, as has been found in samples from US natural populations (reviewed by Larracuente and Presgraves 2012).Unfortunately, no information on the fitness effects of inversion genotypes appears to be available for samples from the ancestral range of D. melanogaster.
Species that are less subject to the problem of recent range expansion, such as D. pseudoobscura, are thus more favorable material.Estimates of the net relative fitnesses and of two measures of viability for D. pseudoobscura 3rd chromosome arrangements in population cages are also given in Supplementary File 4 (the horizontal lines separate data from different experiments).AR and CH cover approximately 30 and 20% of chromosome 3 (the homolog of chromosome arm 2R of D. melanogaster), respectively (Powell 1992), and are thus suitable for comparison with the theoretical prediction for In(3R)P given above.For the viability experiments of Dobzhansky (1947), involving CH and ST from a California population, CH had a frequency of between 0.20 and 0.35 in the original population, but which varied substantially over the year (Wright and Dobzhansky 1946, Fig. 2).The larval viability of CH/ CH individuals was 86% of that of ST/ST individuals, which is reasonably consistent with expectations for an inversion with a frequency of around 0.3.The relative viabilities of CH/CH and AR/ AR individuals were similar, as expected from the similar frequencies of CH and AR in this population.
The results for egg-to-adult viability measurements that used a balancer chromosome to extract 3rd chromosomes from a natural population of D. pseudoobscura (Crumpacker and Salceda 1968) show smaller effects than those predicted from mutational load; only the results for the 2 most common arrangements are shown in the table, with AR and CH having frequencies of approximately 0.5 and 0.28 in the population; the ratio of viabilities of CH/CH vs AR/AR in crosses between carriers of independently extracted chromosomes is 0.98, which is unlikely to differ significantly from 1.
The net fitness estimates obtained from population cage experiments on the natural population used for Dobzhansky's viability estimates, but which segregated for ST, AR, and CH (Wright and Dobzhansky 1946), show patterns that are inconsistent with the mutational load predictions; in particular, AR/AR has a much lower fitness than CH/CH despite their similar frequencies.In contrast, the ratios of net female fitnesses for CH/CH to AR/AR estimated by Anderson and Watanabe (1997) for a laboratory population derived from the same population but segregating only for CH and AR is 0.85, which is consistent with the fact that the population equilibrated at about 25% CH.However, the ratio for ST/ST vs AR/AR in another experimental population was 0.70, despite the fact that ST is usually at least as frequent as AR (Powell 1992).
The overall conclusion from these analyses of the relative performance of rare vs common arrangements is that some measurements fit the expectation of a larger equilibrium mutational load for the less frequent Drosophila inversions, but that the overall patterns imply that other factors obscure the contribution of load to homokaryotype fitnesses.

Can mutational load create a selective advantage to heterokaryotypes?
The theoretical results described earlier make it clear that the mutational load model used here can only create a heterokaryotypic advantage when In and St are present at nearly equal frequencies (see Figs. 1 and 2, and Supplementary Tables 1 and  2).Furthermore, the size of such an advantage is likely to be small, even in subdivided populations, except for large inversions.For example, with x = 0.5, h = 0.25 and 10 5 selected sites, the selection coefficient against both homokaryotypes vs the heterokaryotyope, as given by Equations (3), is 0.0016 with ̅ γ = 1,000, and 0.0010 with ̅ γ = 4,000.However, if the above model of 3 × 10 6 sites with ̅ γ = 1,000 and 10 6 sites with ̅ γ = 4,000 is used, the selection coefficient becomes 0.057.With x = 0.4, the selection coefficients under this model are t 1 = 0.13 for In/In and t 2 = 0.0013 for St/St.For lower values of x, the selective advantage to

Inversion polymorphisms and fitness | 13
In/St heterokaryotypes over St/St is replaced by a selective disadvantage, as shown in Figs. 1 and 2.
The reason for this behavior is that a rarer arrangement accumulates a larger mutational load than its counterpart, due to the lower efficacy of selection with smaller N e for mildly deleterious, partially recessive mutations (see Introduction).When a haplotype from the In population is made heterozygous with a haplotype from the St subpopulation, there is a smaller expected number of heterozygous mutations in the In/St individuals than in the In subpopulation.This means that t 1 > 0 if h < ½, as is evident from Equation (3a), where F 1 , 1 -F 1 , and −C 12 are all positive, as is the difference in the expected frequency of mutations, 〈δq〉, between the In and St subpopulations.In contrast, In/St individuals have a larger expected number of mutations than the St subpopulation.As can be seen from Equation (3b), if h < ½ and the magnitude of 〈δq〉 is sufficiently large, t 2 is negative.But if the 2 arrangements are equally frequent, 〈δq〉 = 0, and the remaining terms guarantee that t 1 = t 2 > 0 when h < ½.In all cases, t 1 = t 2 = 0 if h = ½.The only surprising aspect of the theoretical results is that t 2 is so sensitive to the effect of the relative subpopulation sizes on 〈δq〉.
These theoretical predictions can be compared with the data in Supplementary File 4. The 2 studies of viability in D. melanogaster showed no difference in viability between chromosomal heterozygotes that were free of inversions and chromosomal heterozygotes where at least 1 of the pairs of chromosomes involved carried an inversion.Since the inversions were rare, most of the latter cases will have involved heterokaryotypes, so the lack of any difference is consistent with an absence of heterokaryotypic superiority, as expected for rare inversions.For D. pseudoobscura, the results on net fitness and on viability from population cage experiments indicate strong heterokaryotypic advantages, much larger than the theoretical predictions for inversions of the size involved here.In contrast, the balancer cross data on viability showed a small heterokaryotypic advantage (Crumpacker and Salceda 1968), consistent with the theoretical predictions for equally frequent arrangements.Mérot et al. (2020), Huang et al. (2022), andPei et al. (2022) found no evidence for heterokaryotypic superiority for several fitness components in seaweed flies, sunflowers, and zebra finches, respectively.Overall, it seems likely that the mutational load model may contribute modestly to heterokaryotypic superiority for inversions that are at intermediate frequencies, but cannot explain the large net fitness effects seen in the D. pseudoobscura or C. frigida inversions.This parallels the finding that mutational load is unlikely to provide a selective advantage to new autosomal inversions in randomly mating populations (Nei et al. 1967;Connallon and Olito 2021;Jay et al. 2022).

Population genomic indicators of a reduced efficacy of selection on low-frequency arrangements
As described in A single population: population genomics statistics, the two most informative population genomic statistics concerning a reduced efficacy of selection on a low-frequency arrangement are the mean frequencies of variants at the selected sites (the ̅ q i ) and the ratios of diversities at selected and neutral sites (the π seli /π neuti ).The measure of skew toward low-frequency variants relative to neutral expectation (Δθ wi ) is much less sensitive to subpopulation size, unless the dominance coefficient is implausibly small.In addition, measures of skew are sensitive to population size changes (Tajima 1989), and must be treated with caution when making inferences about selection.For this reason, only the other two statistics will be considered here.
A problem with using the ̅ q i is that these are not directly observable unless bioinformatic methods for classifying variants as deleterious are used; simply using the mean frequency of derived nonsynonymous variants in a sample as a proxy (cf., Campos et al. 2014) is not necessarily adequate when selection is weak, since it does not take into account fixed sites.Stenløk et al. (2022) used PROVEAN scores to estimate the mean numbers of deleterious missense mutations in inverted and standard arrangements of Atlantic salmon, but found no significant differences; an enrichment of small indels in the large (3.09Mb) Chr18 inversion was, however, detected.The frequencies of this inversion are, however, highly variable between populations, so it is not clear how to interpret this difference.
There are, however, theoretical reasons for expecting the ratio ̅ q 1 /̅ q 2 to behave similarly to the ratio R π = (π sel1 /π neut1 )/ (π sel2 /π neut2 ), at least when h is not too small.Welch et al. (2008) analyzed the properties of the ratio of selected site to neutral site diversity under a similar model to that used here, assuming h = 0.5 and a gamma distribution of selection coefficients (see also James et al. 2017).With ̅ γ >> 1/a, where a is the shape parameter of the gamma distribution in Equation ( 9), the following approximate relation holds: where k is a constant that depends on the effective population size and the mean strength of selection.
In the present case, this relation implies that While we would not expect this relation to be exact for h < 0.5, it is plausible to assume that it would be a reasonably good approximation when h is not too close to 0. The argument used by Welch et al. (2008) also implies that a similar relation should apply to the mean frequencies of deleterious mutations: The accuracy of these approximations can be tested using the numerical results in Supplementary Tables 1 and 2. Table 3 gives some examples for the case of a single randomly mating population, showing that the two ratios on the left-hand sides of Equations ( 13b) and (13c) behave very similarly as functions of h, with R ln(q) > R ln(π) when h < 0.25, approaching the prediction on the right-hand sides of the equations for h ≥ 0.35.Similar results apply to the case of a subdivided population.These results suggest that R π can used as a conservative proxy for ̅ q 1 /̅ q 2 , at least for the case of a gamma distribution of selection coefficients with ̅ γ >> 1/a and an intermediate dominance coefficient.An objection to using the correlation between R ln(π) and a measure of neutral diversity such as synonymous site diversity to investigate whether the efficacy of purifying selection declines with N e is that π neut is the denominator of π sel /π neut .In addition to the statistical problem of a negative correlation introduced by this relationship, discussed by James et al. (2017), sites under sufficiently strong purifying selection could maintain a constant diversity across different N e values (Campos et al. 2014).If this were the case, the expectation of π sel /π neut would simply be proportional to 1/π neut , and we would then have R π = π neut2 /π neut1 .As described in A single population: population genomics statistics, this is not the case with the model used here; R π with x < 0.5 is always less than y/x.For example, with ̅ γ = 4,000 and h = 0.25, R π = 2.00 for x = 0.1 (y/x = 9) and R π = 1.33 for x = 0.3 (y/x = 2.33).This indicates that R π provides a signal that purifying selection is weakened by smaller subpopulation size.Comparisons of this kind could easily be done using real data.
Overall, therefore, these considerations suggest that, despite the above reservations, R π is quite a useful index of the efficacy of purifying selection, and that one might expect R ln(π) to be approximately equal to a ln(π neut2 /π neut1 ).A related principle was used by James et al. (2017) and Castellano et al. (2018) to test for relations between the efficacy of purifying selection and N e for animal mitochondrial genomes and different regions of the D. melanogaster nuclear genome, respectively.Unfortunately, there appears to be relatively little relevant information for autosomal inversions, other than cases such as the mimicry supergene in Heliconius numata (Jay et al. 2021) and the behavioral supergene of the white-throated sparrow Zonotrichia albicolis (Jeong et al. 2021), which are largely maintained as heterozygotes due to negative assortative mating.These systems are analogous to sex chromosomes, where one arrangement is permanently heterozygous and effectively lacks recombination.There is thus likely to be intense Hill-Robertson interference (Charlesworth and Charlesworth 2000), which would greatly reduce the efficacy of selection below the simple effect of a lower subpopulation size.This is consistent with the strongly elevated π sel /π neut values found for the H. numata inversions; Z. albicolis showed, however, only a modest effect.Jay et al. (2021) also found a large increase in the density of transposable elements (TEs) in the H. numata inversions and interpreted this as evidence for an increased mutational load.The accumulation of TEs in low recombination regions of genomes, including low-frequency Drosophila inversions (Sniegowski and Charlesworth 1994), has long been documented (Charlesworth et al. 1994).Most insertions are found in intergenic regions, where direct selective effects are likely to be weak, and where ectopic exchange inducing deleterious chromosome rearrangements is probably a major factor in causing their elimination.It is thus likely that a reduced frequency of ectopic recombination is the major factor in causing higher densities of TE insertions in such cases (Charlesworth et al. 1994), so that this phenomenon cannot be taken as evidence for an increased mutational load.

What strength of selection has the main effect on the load and population genomic statistics?
Another question raised by the theoretical results is: what part of the distribution of selection coefficients contribute to the differences between the In and St subpopulations?Stronger selection reduces the frequencies of deleterious mutations and their chances of fixation within a subpopulation, but also increases the sizes of any resulting loads.The major contribution to the relevant load statistics is thus likely to come from selection coefficients that are neither too large nor too small.This expectation can be tested by examining the contributions from the different zones described in A single population: obtaining the mutational load and population genomic statistics.These are shown in Supplementary Table 1 for the case of a single population.These results show that the major contributions to the L i and t i come from zone 2a, defined by an intermediate intensity of selection (for more details, see Supplementary Table 5).For example, for a randomly mating single population with ̅ γ = 4,000, x = 0.1, and h = 0.25, zones 1, 2a, 2b, and 3 contribute 4, 36, 24, and 36%, respectively, to the distribution of selection coefficients against deleterious mutations.Zone 2a alone gives L 1 = 0.0375, L 2 = 0.0012, t 1 = 0.0270, and t 2 = −0.0090,compared with net values of L 1 = 0.0381, L 2 = 0.0016, t 1 = 0.0271, and t 2 = −0.0090.It covers the interval (0.278, 463) of γ for the whole population.
The finer dissection of the distribution of selection coefficients used in the case of a subdivided population reveals that the so-called quasineutral zone 2 in this case (see Section 5 of Supplementary File 1) contributes most to the load statistics.For example, with ̅ γ = 4,000, x = 0.1, h = 0.25, and F ST = 0.05 (M = 19), zones 1, 2, 3, and 4 for the subdivided case contribute 4, 17, 21, and 58%, respectively.Zone 2 covers the γ interval (0.278, 55.6), and contributes L 1 = 0.0344, L 2 = 0.0011, t 1 = 0.0249, and t 2 = −0.0082,compared with values of L 1 = 0.0401, L 2 = 0.0019, t 1 = 0.0285, and t 2 = −0.0094for the whole distribution.In this case, with 200 demes of size 5,000 each, both the In and St subpopulations behave as effectively neutral within demes for this part of the distribution of selection coefficients.A similar pattern holds even with F ST = 0.25.Nevertheless, the results are very similar to those for the single population, implying that the expectations of the population genetic parameters for an island model with a relatively low level of isolation between demes are mainly controlled by the size of the metapopulation, as was shown to be the case analytically by Wakeley (2003).Intuitively, this reflects the fact that, under these conditions, a mutation spends relatively little of its total sojourn time in the deme in which it arose.

Confounding factors
The present study assumes that selected sites are at statistical equilibrium under mutation, selection and drift, an absence of recombinational exchange between In and St in heterokaryotypes, and complete independence among sites under selection.These Table 3.Values of the ratios with respect to In vs St of the natural logarithms of π sel /π neut (R ln(π) ) and ̅ q (R ln(q) ) for the case of a single population.
̅ γ = 1,000 ̅ γ = 4,000 The predicted values of the ratios with the shape parameter a = 0.3 are 0.659 for x = 0.1 and 0.254 for x = 0.3.The mutational parameters and population size are the same as in Fig. 1.
assumptions are likely to violated in many real-life situations.First, consider the question of departure from equilibrium.If the inversion arises as a unique mutational event, as the available evidence seems to suggest (reviewed in Charlesworth 2023), the In subpopulation will initially completely lack genetic variability, and much time will be needed for equilibrium to be approached.
The St subpopulation can be assumed to be close to equilibrium initially and will thus approach its new equilibrium much faster than the In subpopulation, so that only the latter need be considered here.For this subpopulation, the magnitudes of the L i , B i , and t i will be below their equilibrium values for a long time after the inversion has approached its equilibrium frequency under balancing selection.
It is difficult to make exact predictions about the rate of approach to statistical equilibrium when both drift and selection play a role, which has been shown above to be the situation that contributes the most to the selective differences among karyotypes.For the limiting case of complete neutrality, it is known that the divergence of the expected nucleotide site diversity from its equilibrium value at time t is equal to the product of its initial value and exp(−t/2N e ) in a randomly mating population (Malécot 1969, p. 40), so that the timescale for approach to equilibrium is of the order of 2N e generations.For the other limiting case of fully deterministic evolution, with γ >> 1, the divergence at time t of the mutant allele frequency q from its equilibrium value of u/hs is approximately equal to the product of its initial value and exp(−t/hs).For the intermediate situation when γ is not much greater than 1, the two measures of the rate of approach to equilibrium are not very different, so it is plausible to assume that the true rate lies between them.The time needed to approach equilibrium with respect to mutations that have the largest effect on the load statistics is thus likely to be < 2Nx generations under a Wright-Fisher model.For a population with an N e of 10 6 and x = 0.1, this would correspond to 10 5 generations, i.e. about 10,000 years for a species like D. melanogaster with approximately 10 generations per year.For x = 0.5, about 50,000 years would be required.
A problem with assessing whether these timescales are consistent with data on inversion polymorphisms is that the occurrence of recombinational exchange between In and St due to gene conversion and/or double crossing over (reviewed by Krimbas and Powell 1992;Korunes and Noor 2019), means that estimates of inversion age based on sequence divergence between different arrangements tend to produce underestimates of the age of the derived arrangement (Charlesworth 2023).A variety of lines of evidence suggests, however, that inversions such as In(3R)P are close to selection-drift-recombination equilibrium with respect to neutral variability (Charlesworth 2023); since selection against deleterious mutations will cause a faster approach to equilibrium than with neutrality, it is likely that the load statistics will also be close to their equilibrium values, unless demographic disturbances have caused serious perturbations.
However, the theory developed here for predicting mutational loads has ignored recombination.If the estimate of a typical rate of gene conversion of 10 -5 in female meiosis in inversion heterokaryotypes in Drosophila (Chovnick 1973;Korunes and Noor 2019) is accepted, the effective rate is 0.5 × 10 -5 due to the absence of exchange in males.An hs value that somewhat exceeds 10 -5 would thus be sufficient to overcome the effects of gene flow between In and St; with h = 0.25 and N T = 10 6 , this would correspond to γ > 80, which lies outside the range of γ values that contribute to a noticeable difference in load statistics between In and St with x = 0.1, as discussed in the previous section.It is therefore likely that recombination will significantly reduce such differences, providing another reason for regarding the above estimates as providing upper bounds to the predictions.
In addition, the fitness effects of any deleterious mutations that were present on the initial inversion haplotype have been ignored; however, the final equilibrium state considered here, where reverse mutations have been included, means that such effects will have been removed.During the approach to equilibrium they must, of course, play a role in reducing any selective advantage to a new inversion (Nei et al. 1967;Connallon and Olito 2021;Jay et al. 2022).
The low frequency of recombination in inversion heterokaryotypes for a rare inversion may create Hill-Robertson interference effects, as noted above in Population genomic indicators of a reduced efficacy of selection on low-frequency arrangements.Berdan et al (2021) used simulations to study this possibility, and found strong interference effects under the assumption that mutations were completely recessive.As pointed out previously, this assumption is implausible, even for mutations with large homozygous fitness effects.In addition, recombination between In and St will reduce interference effects.These will tend to increase the load in the less frequent subpopulation and in heterokaryotypes, and so will not contribute to the maintenance of the inversion polymorphism.There is little evidence for any effects of Hill-Robertson interference on the available population genomic statistics for the D. melanogaster inversions (Charlesworth 2023).

Conclusions
The theoretical results described here show that a longmaintained autosomal inversion polymorphism with no recombination in heterokaryotypes may develop a substantially higher mutational load for the less frequent arrangement.The magnitude of the difference between arrangements can be large for rare polymorphic inversions of the size usually encountered in Drosophila populations, but declines quickly as the frequency of the rare arrangement increases.It is also strongly influenced by the abundance of relative weakly selected noncoding sequences, since drift acts more strongly on these than on strongly selected nonsynonymous mutations.A selective advantage to heterokaryotypes is only expected when the two alternative arrangements are nearly equal in frequency, and is likely to be small even in this case.Experiments on the effects of several Drosophila inversion polymorphisms on fitness components give inconsistent results, although mutational load may contribute to some of the effects that have been detected.It should also be possible to detect molecular signatures of an increased load, such as an enhanced ratio of nonsynonymous to synonymous nucleotide site diversities, but the data are currently too scanty to draw firm conclusions.The effects of recombinational exchange in heterokaryotypes and Hill-Robertson interference, which oppose each other, were ignored here, and deserve further study.

Fig. 1 .
Fig. 1.The mutational load statistics for a single population of size N = 10 6 with 10 5 selected sites, plotted against the dominance coefficient h.The results for 3 different frequencies of the inversion are shown.The selection coefficients follow a gamma distribution with a mean of 5 × 10 −4 and a shape parameter of 0.3.The mutation rate to deleterious alleles is 5 × 10 -9 per basepair, with a mutational bias toward deleterious variants of 1.5.The upper and lower panels have products of 2N and mean s of 4,000 and 1,000, respectively.The filled and half-filled symbols denote values for the In and St subpopulations, respectively.The lozenges are the net mutational loads within the respective karyotypes, and the crosses are the corresponding inbreeding loads.The triangles are the selection coefficients against homokaryotypes relative to the heterokaryotype.Only the values for In are shown when x = 0.5. 2 .

Fig. 2 .
Fig. 2. The mutational load statistics for a single populations of size N = 2 × 10 6 (upper panel) and N = 10 5 (lower panel), with the same mean selection coefficient (0.001), plotted against the dominance coefficient h.The results for 3 different frequencies of the inversion are shown.The other parameters are the same as in Fig. 1.The filled and half-filled symbols denote values for the In and St subpopulations, respectively.The lozenges are the net mutational loads within the respective karyotypes, and the crosses are the corresponding inbreeding loads.The triangles are the selection coefficients against homokaryotypes relative to the heterokaryotype.Only the values for In are shown when x = 0.5.

Fig. 3 .
Fig. 3.The population genomics statistics for the St subpopulation, plotted against the dominance coefficient h.The evolutionary parameters are the same as in Fig. 1.The results for 3 different frequencies of the inversion are shown.The half-filled triangles are the mean frequencies of A 2 and the half-filled lozenges are the mean diversities at selected sites.The filled triangles are the ratios of mean diversities at selected sites to mean diversities at neutral sites.The filled lozenges are the values of Δθ w for the selected sites.Note that the distinction between filled and half-filled symbols does not refer to the In vs St subpopulations, as was the case for the previous figures.

Table 1 .
Definitions of the most important symbols used in the text.
T Subscripts i = 1 and i = 2 are used to denote parameters applicable to the In and St subpopulations, respectively.Subscript d is used to denote within-deme parameters.