Across many species where inversions have been implicated in local adaptation, genomes often evolve to contain multiple, large inversions that arise early in divergence. Why this occurs has yet to be resolved. To address this gap, we built forward-time simulations in which inversions have flexible characteristics and can invade a metapopulation undergoing spatially divergent selection for a highly polygenic trait. In our simulations, inversions typically arose early in divergence, captured standing genetic variation upon mutation, and then accumulated many small-effect loci over time. Under special conditions, inversions could also arise late in adaptation and capture locally adapted alleles. Polygenic inversions behaved similarly to a single supergene of large effect and were detectable by genome scans. Our results show that characteristics of adaptive inversions found in empirical studies (e.g., multiple large, old inversions that are FST outliers, sometimes overlapping with other inversions) are consistent with a highly polygenic architecture, and inversions do not need to contain any large-effect genes to play an important role in local adaptation. By combining a population and quantitative genetic framework, our results give a deeper understanding of the specific conditions needed for inversions to be involved in adaptation when the genetic architecture is polygenic.

We simulated a two-patch, Wright–Fisher, forward-time simulation using SLiM [v 3.6; Haller & Messer 2019]. We simulated 21 linkage groups (LGs), each 10 centimorgans (cM) long, with a resolution of 0.0001 cM between proximate bases in SLiM. In 20 of the LGs QTN mutations could arise, and the 21st was neutral. Our demes consisted of N = 1000 diploid individuals in each deme, and the population-scaled mutation rate N_eμ for QTNs whose effect size was drawn from a normal distribution with a mean of 0 and standard deviation of ?_m depended on the parameter level: for “polygenic” N_eμ = 2e-5 and ?_m = 0.2 (lower mutation rate and larger effect size of mutations); for “highly polygenic” N_eμ = 2e-4 and ?_m = 0.002 (higher mutation rate and smaller effect size of mutations); which are on the scale that lead to evolution of polygenic architectures. This was a quantitative genetic model where individuals were evolving to a phenotype optimum in each deme with individual fitness determined by a guassian fitness function. The inversion mutation rate per genome was drawn from the binomial distribution with a single draw and a probability of an inversion mutation of ?_inv = 0.001, which equated to, on average, two inversions in the metapopulation per generation. The length of the inversion mutation was drawn from a discrete uniform distribution [100, 50,000] (i.e., up to half the linkage group length) and the location was drawn randomly from the 20 linkage groups that were undergoing selection (inversions did not arise on the neutral linkage group). In inversion heterozygotes, crossovers in the inversion region were suppressed (r = 0); individuals homozygous for the inversion underwent recombination at rate r as usual. We simulated multiple migration rates, selection strengths, and with and without environmental variance added to individual phenotypes. Results from these simulations were used in a custom R script to identify adaptive inversions and test genome scan outlier methods on their ability to accurately detect inversions. In addition, this custom R script was used to create a number of plots provided in this data archive.

README files were created for understanding the output files from SLiM and for each custom R script.