Assessing the role of inversions in maintaining genomic differentiation after secondary contact: local adaptation, genetic incompatibilities, and drift

Due to their effects on reducing recombination, chromosomal inversions may play an important role in speciation by establishing and/or maintaining linked blocks of genes causing reproductive isolation (RI) between populations. These views fit empirical data indicating that inversions typically harbour loci involved in RI. However, previous computer simulations of infinite populations with 2-4 loci involved in RI implied that, even with gene flux as low as 10−8 between alternative arrangements, inversions may not have large, qualitative advantages over collinear regions in maintaining population differentiation after secondary contact. Here, we report that finite population sizes can help counteract the homogenizing consequences of gene flux, especially when several fitness-related loci reside within the inversion. In these cases, the persistence time of differentiation after secondary contact can be similar to when gene flux is absent, and notably longer than the persistence time without inversions. Thus, despite gene flux, population differentiation may be maintained for up to 100,000 generations, during which time new incompatibilities and/or local adaptations might accumulate and facilitate progress towards speciation. How often these conditions are met in nature remains to be determined.

initially rare inversion, assuming that the inversion polymorphism captures two or 85 more adaptive loci undergoing divergent selection. A general finding of these studies 86 is that, by suppressing recombination, an inversion protects locally favorable 87 combinations of alleles from mixing with maladapted alleles introduced by migration 88 (except when migration is too strong, see Eq. A9 in Charlesworth and Barton, 2018). 89 In addition, Nei et al. (1967)  Our goal is to test how the effects of more than two loci and of finite population sizes 155 may influence the erosion of differentiation within an inversion polymorphism after 156 secondary contact, paying particular attention to the recent empirical results implying 157 higher rates of gene flux between chromosomal rearrangements than previously 158 thought (Korunes and Noor 2019). 159 160 We show that gene flux between alternative chromosomal rearrangements can be 161 retarded by the finite size of populations (with up to 100,000 individuals in some 162 cases) coming into secondary contact. This is due to a longer waiting time for actual 163 gene flux events to be realized, and also, to some extent, due to the increased 164 probability of the loss of favored allele combinations by random genetic drift. In 165 particular, we find that under these conditions, inversions can retain longer-lived 166 differentiation even with relatively high levels of gene flux (2 • 10 !$ per gamete per 167 generation) when many loci of small fitness effects reside within the inversion 168 polymorphism. Importantly, such inversion polymorphisms can protect population 169 differentiation for tens of thousands of generations (up to 100,000, or longer in some 170 cases) after secondary contact. This time span may provide a window of opportunity 171 for additional differences to accumulate (further local adaptation, more 172 incompatibilities, reinforcement, etc.), facilitating further progress towards speciation. 173 How often these conditions are met for hybridizing populations in nature remains to 174 be determined. Nevertheless, our results do offer a region of parameter space where 175 inversions not containing loci causing fitness trade-offs (divergent selection) in 176 alternate habitats may play an enhanced role in facilitating speciation than collinear 177 regions following secondary contact. 178

Model and Methods 179
To assess the role of multiple loci and finite population sizes in maintaining 180 population differentiation after secondary contact, we used individual-based computer 181 simulations of two diploid populations with finite, constant population size (denoted 182 by N below) and discrete, non-overlapping generations. At the start of each 183 simulation, the two populations (referred to as Population 1 and Population 2) were 184 fixed for different chromosomal arrangements. For convenience, we refer to the 185 arrangement in Population 1 as standard and to the other one as inverted. We 186 modelled 2L bi-allelic loci within the region polymorphic for the inversion (see 187 below) and assumed that at the start of each simulation the two populations were fixed 188 for alternative alleles at each locus. At generation 0, both populations consisted of 189 adult virgin individuals coming into secondary contact. Thereafter, adult virgin 190 individuals underwent migration, followed by soft selection, recombination and 191 mating (always locally in each population). We neglected mutations throughout. In homokaryotypes (either for the standard or for the inverted chromosomal 201 arrangement), consecutive loci in the simulated genomic region underwent 202 recombination with a per-gamete per-generation probability of %& . We modelled the 203 recombination rate %& in two different ways, that we explain below. 204 205 First, we focused on inversions of the same size in terms of their total map distance. 206 Throughout, we refer to this version of the model as the conserved-size model. Here, 207 we chose the recombination rate %& in such a way that it was smaller when the 208 number of fitness-related loci (i.e., 2L) in the inversion polymorphism was larger, and 209 that the total recombination distance (denoted by below) spanned by the inversion 210 was independent of the number of loci it contained. In particular, and in line with 211 Feder and Nosil (2009), when 2 = 2, we set % = = 0.1 , whereas for 2 > 2, we 212 used %& = 1 − (1 − ) '/(%&!') . In addition, due to our focus on the inversion, gene 213 flux between the alternative arrangements in heterokaryotypes occurred with a per-214 gamete per-generation probability Inv , that was independent of the number of loci 2 . 215 Gene flux was implemented essentially as in Feder and Nosil (2009), except that we 216 assumed that gene flux occurred 50% of the time by gene conversion and 50% by 217 double crossover (whereas the corresponding ratio in Feder and Nosil (2009) was 218 30:70). This difference takes into account a recent study by Korunes and Noor (2019)  219 showing that double crossover in hybrids between Drosophila pseudoobscura and D. 220 persimilis result in slightly lower gene flux per site and generation than gene 221 conversion; but the reverse is observed for double crossovers spanning a region larger 222 than 1 Mb. Since in our simulations the ratio between the number of exchanged alleles 223 due to double crossover and gene conversion varies depending on how many loci are 224 exchanged in each double crossover (see below), we assumed equal rates of double 225 crossover and gene conversion to avoid any bias towards one or the other process. In 226 our model, a double crossover was assumed to occur within a region containing any 227 number between 1 and 2L of consecutive loci, each combination being equally likely. 228 Gene conversion events were modelled so that an allele at a given locus at a given 229 chromosome was altered to the allele at the same locus residing at the homologous 230 chromosome. If the two alleles were identical prior to gene conversion, they remained 231 identical also afterwards. Gene conversion events occurred at a single locus, each 232 locus being equally likely. 233

234
In a second version of the model, we assumed that the size of the inversion increased 235 with the number of fitness-related loci within it. We refer to this version of the model 236 as the increasing-size model. Here, we set the recombination rate %& between each 237 pair of consecutive loci in homokaryotypes to be independent of the number of 238 fitness-related loci within the inversion, i.e., %& = = 0.1. In this version of the 239 model, the total rate of gene flux in heterokaryotypes (denoted by Inv,2& below) was 240 assumed to depend linearly on the number of fitness-related loci, i.e., Inv,2& = Inv , 241 where Inv had the same value as in the conserved-size model explained above. Note 242 that this means that when the inversion contained exactly two fitness-related loci (i.e., 243 2 = 2), the conserved-size model and the increasing-size model were equivalent, 244 whereas this was not true for 2 > 2. This setting allowed us to assess the impact of 245 increasing the number of fitness-related loci within the inversion beyond two while 246 either keeping the total map distance of the inversion constant or increasing it 247 (roughly linearly) with the number of fitness-related loci within the inversion. 248 249 We further assumed that mating was random within each local population, that is, the 250 pool of 2N gametes obtained after selection and recombination was randomly divided 251 into N pairs, thus producing N offspring in each population (this corresponds to 252 sampling without replacement). After reproduction, all adult individuals were replaced 253 by their offspring, which were then treated as the next generation of adults. 254 255 As mentioned above, we studied two out of the five models analyzed by Feder and 256 Nosil (2009), specifically their Model 2 that involves alleles that are locally favored in 257 one habitat and neutral in the other, and Model 4 that involves universally beneficial 258 alleles with negative fitness epistasis (i.e., incompatibilities). In both models and in 259 line with Feder and Nosil (2009), we considered two types of loci, denoted here by 0 260 and 0 for = 1,2, … , (the difference between the two types of loci is explained 261 further below). Initially, all loci were bi-allelic: the two alleles at locus 0 ( 0 ) are 262 denoted by 0 and 0 (that is, 0 and 0 for locus 0 ). Furthermore, the alleles denoted 263 by uppercase letters at each locus were initially fixed in Population 1, whereas the 264 alternative alleles (denoted by lowercase letters) were initially fixed in Population 2. 265 266 Selection was implemented in the two models as follows. In Model 2, the alleles 267 denoted by uppercase letters at loci 0 ( = 1,2, … , ) were advantageous over the 268 alternative alleles at these loci in the first population, whereas both allelic types at 269 these loci were selectively neutral in the second population. The opposite was true for 270 loci 0 ( = 1,2, … , ): whereas at these loci both allelic types were selectively neutral 271 in the first population, the alleles denoted by lowercase letters were advantageous over 272 the alternative alleles in the second population. We considered two different models 273 with respect to how loci 0 and 0 were ordered in the genomic region simulated. Genetic incompatibilities were modelled as occurring between universally favored 303 alleles 0 and 0 at pairs of loci 0 and 0 (note the same index ). We accounted for 304 incompatibilities slightly differently from Feder and Nosil (2009). Namely, for each 305 individual we counted the number of incompatibilities at each pair of loci 0 and 0 as 306 follows. (1) If an individual is a homozygote for allele 0 and a heterozygote at locus 307 0 , we accounted for two incompatibilities at this pair of loci; (2) for a homozygote 308 for allele 0 , and a homozygote for allele 0 , we accounted for four incompatibilities 309 for this pair of loci; (3) for an individual that is a heterozygote at both loci 0 and 0 , 310 we counted one incompatibility; and (4) for an individual that is a heterozygote at 311 locus 0 but a homozygote for allele 0 at locus 0 , we accounted for two 312 incompatibilities at this pair of loci. All other allele combinations at these two pairs of 313 loci gave zero incompatibilities. For each individual we thus computed the number of 314 incompatibilities at each pair of loci 0 and 0 ( = 1,2, … , ), and summed them up 315 to obtain the total number of incompatibilities (denoted by inc ) that the individual 316 carried. To each individual that carried at least one incompatibility (i.e. inc ≥ 1), we 317 assigned a fitness of (1 − %&) 3 inc , with %& denoting the strength of the negative 318 epistasis per incompatibility. Note that, because we assumed that incompatibilities 319 may arise between pairs of loci, the maximum number of incompatibilities in the 320 model with 2 loci is equal to 4 , and this is obtained when an individual is a 321 homozygote for alleles 0 and for alleles 0 at all loci 0 and 0 . To assure meaningful 322 comparison between the models with different number of fitness-related loci within 323 the inversion, we set %& = 1 − (1 − ) '/& , so that when 2 = 2, we have % = 324 . Note that our implementation of incompatibilities differs from that in Feder and 325 Nosil (2009), where it was assumed that an individual carrying any combination of 326 incompatible alleles has a fitness disadvantage equal to 1 − (independently of how 327 many incompatibilities it has and of their linkage phase on chromosomes, as long as 328 there is at least one). Similar to the model with locally favored and neutral alleles, 329 here we also analyzed two different orderings of the loci 0 and 0 , that is, we 330 considered the half-half array ordering, and the consecutive array ordering (see 331 above). 332 333 The parameter values we tested corresponded largely to those tested by Feder and 334 Nosil (2009), but here we focused on a much larger rate of gene flux (by four orders 335 of magnitude) than that used by them, and even larger in the increasing-size version of 336 the models with more than two loci (see above). In addition, and unlike  Table I, as well as  343   Tables S1-S2.
(1) 365 366 Here, the factor 1/(2 ) serves to average over all fitness-related loci, and index 367 accounts for all sampling time points when the recorded allele-frequency difference 368 Δ 9 ( 0 ) (or Δ 9 ( 0 )), at a given locus 0 (resp. 0 ) was nonzero (here, Δ 9 ( 0 ) 369 represents the frequency difference of allele 0 between Population 1 and Population 370 2; and similarly for Δ 9 ( 0 )). The last such sampling point for locus 0 ( 0 ) is 371 denoted < ' (resp. = ' ). Note that in the absence of mutations, which we assumed, once 372 the allele-frequency difference between the two populations becomes zero, it stays 373 zero infinitely. However, we did not observe fixation in many cases we simulated 374

Results 388
In Model 2 involving neutral and universally beneficial alleles, our simulations 389 showed that when the fitness-related loci were captured by an inversion 390 polymorphism (either with or without gene flux), the persistence time of 391 differentiation weighted by the differentiation after secondary contact ( w ; see Model 392 and Methods) was longer than in the absence of inversion polymorphism: the 393 difference between the two ranged from thousands to tens of thousands of generations 394 (Figs. 1-2, and see Figs. S1-S10; but note that the actual persistence time of 395 differentiation after secondary contact was larger than w , due to the weighting we 396 made when defining w ). This was true both when gene flux between the alternative 397 arrangements was absent, as well when it was relatively high (in particular, four 398 orders of magnitude higher than that assumed in Feder and Nosil, 2009  with increasing the number of fitness-related loci within the inversion, with w usually 406 reaching a plateau after six, or more loci involved. This was true for all parameter 407 values and variations of the model we considered, but the advantage was typically 408 stronger for lower than for higher migration rates, and in the conserved-size than in 409 the increasing-size variation of the model (compare Fig. 1 and Fig. 2). Our 410 simulations showed that there were also some differences between the model 411 variations with respect to the ordering (half-half array or consecutive array) of the 412 two types of loci ( 0 that carried locally beneficial alleles in Population 1, and 0 that 413 carried locally beneficial alleles in Population 2), but these differences were subtle. between consecutive loci was less than the migration timescale). 500

501
The patterns described for the model involving universally beneficial alleles with 502 genetic incompatibilities were also reflected in ST (and 67 ) patterns from individual 503 realizations of the model (Figs. S12-S13, S15-S16, S18-S19, S21-S22). As in the 504 model with locally beneficial and neutral alleles, stochastic fluctuations between 505 individual realizations were high (not shown). 506 507 The findings outlined above were retained when considering populations of larger size 508 ( = 5,000; Figs. S23-S24), but the advantage of an inversion polymorphism in 509 maintaining population differentiation after secondary contact was slightly less than 510 when = 1,000 (compare Fig. 1a-b to Fig. S23a-b, and Fig. 3 to Fig. S24). 511 To further assess how population size impacts on the advantage of inversions over 512 collinear regions in maintaining population differentiation after secondary contact, we 513 performed simulations with population sizes of up to 100,000 individuals for a subset 514 of parameter values, focusing on cases with 20 loci (Figs. 7-10). We found that the 515 advantage of inversions over collinear regions was either fully retained when we 516 increased the population size ( Fig. 8a-d, 9a, 10a, 10c-d) or was less but still 517 noticeable. However, in the consecutive-array, conserved-size version of the model 518 involving universally beneficial alleles with genetic incompatibilities, population 519 differentiation was strong and maintained throughout the entire simulated time period 520 after secondary contact (100,000 generations) either with or without inversions ( One limitation of our study shared with earlier work is that we considered a two-612 population system because such a simple scenario is convenient for modelling. 613 However, if testable quantitative predictions applicable to empirical systems are to be 614 obtained, it will be necessary to use spatially-explicit models involving many local Our results support the hypothesis that inversions may be a strong mechanism 627 contributing to maintenance of population differentiation, and potentially speciation, 628 when divergence has a polygenic basis (Feder et al. 2005       Butlin for valuable comments on the manuscript. This study was supported by the 795 European