Natural selection drives genome‐wide evolution via chance genetic associations

Understanding selection's impact on the genome is a major theme in biology. Functionally neutral genetic regions can be affected indirectly by natural selection, via their statistical association with genes under direct selection. The genomic extent of such indirect selection, particularly across loci not physically linked to those under direct selection, remains poorly understood, as does the time scale at which indirect selection occurs. Here, we use field experiments and genomic data in stick insects, deer mice and stickleback fish to show that widespread statistical associations with genes known to affect fitness cause many genetic loci across the genome to be impacted indirectly by selection. This includes regions physically distant from those directly under selection. Then, focusing on the stick insect system, we show that statistical associations between SNPs and other unknown, causal variants result in additional indirect selection in general and specifically within genomic regions of physically linked loci. This widespread indirect selection necessarily makes aspects of evolution more predictable. Thus, natural selection combines with chance genetic associations to affect genome‐wide evolution across linked and unlinked loci and even in modest‐sized populations. This process has implications for the application of evolutionary principles in basic and applied science.


| INTRODUC TI ON
Evolution is best understood for single traits and genes, such as the Ectodysplasin gene controlling lateral plates in stickleback fish (Colosimo et al., 2005), the Agouti gene controlling coat colour in mice (Linnen et al., 2013), and Cortex, which controls wing coloration in moths and butterflies (Nadeau et al., 2016). The effects of selection on such traits and genes can, in theory, be transmitted to other (i.e. not directly selected) correlated genetic regions across the genome (i.e. linked selection; Maynard Smith & Haigh, 1974;Barton, 1983Barton, , 2000Flaxman et al., 2013). Thus, even functionally neutral genetic regions can experience selection indirectly due to statistical association (i.e. linkage disequilibrium, LD) with selected regions.
There is increasing evidence that many traits are polygenic (Boyle et al., 2017;Rockman, 2012) such that adaptation often involves selection on many regions of the genome, as documented in rabbits, cichlid fish, great-tits, Rhagoletis flies and Arabidopsis (Bosse et al., 2017;Carneiro et al., 2014;Dowle et al., 2020;Exposito-Alonso et al., 2019;Kautt et al., 2020). In such cases, the indirect effects of selection could be even more pervasive.
Studies of molecular evolution in flies, humans and other organisms have shown that much of the genome can be impacted by indirect selection (i.e. linked selection) (Cai et al., 2009;Elyashiv et al., 2016;Sella et al., 2009). For example, correlations across the genome between genetic diversity and recombination provide evidence of widespread background selection (Begun & Aquadro, 1992;McGaugh et al., 2012;Roselius et al., 2005), and frequent selective sweeps have been documented in Drosophila (Campos et al., 2017;Jensen et al., 2008). However, such analyses of molecular evolution average over long periods of time, and thus, it remains unclear whether indirect selection is pervasive on shorter time scales, such as generations or decades (but see, e.g., Bergland et al., 2014;Buffalo & Coop, 2020). Our best evidence for short-term genetic consequences of selection comes from studies of ecologically important traits determined by one or a few genes, such as in finches, stickleback and mice (Barrett et al., 2008(Barrett et al., , 2019Grant & Grant, 2002, 2006. Strong selection on such traits can have important immediate evolutionary and ecological consequences (Hendry, 2020;Thompson, 2013), but it is unclear if this affects few or many regions of the genome. Thus, the extent and genome-wide consequences of indirect selection on short to moderate time scales are poorly understood. This makes our understanding of evolution incomplete because the genomic distribution of selection affects the dynamics of adaptation and speciation, and the predictability of evolution (e.g. Hahn, 2018;. Indeed, indirect selection could have notable effects on the short-term predictability of evolution (e.g. the direction of allele frequency change), even if it does not increase the magnitude of per-generation change much beyond genetic drift alone. Finally, the genomic consequences of selection have implications for many evolutionary applications, such as demographic and phylogenetic inference, which often assume neutrality (Pouyet et al., 2018;Schrider et al., 2016).
Elucidating the consequences of selection on the genome is challenging because most genetic regions are affected by multiple factors, whose contributions are difficult to isolate without experimental approaches (Buerkle et al., 2011). For example, the LD that mediates indirect selection is affected by past selection, recombination, chance in finite populations (i.e. genetic drift) and population structure (Barton, 2011). Specifically, genetic drift and population structure can generate widespread LD, even among unlinked loci on different chromosomes. In terms of drift, LD arises when combinations of alleles are over-represented by chance. The effect of structure is commonly considered in genome-wide association mapping studies, where even structure within populations sometimes generates sufficient LD to confound genotype-phenotype associations if not properly controlled for (Aranzana et al., 2005;Balding, 2006). Thus, although facilitated by chromosomal linkage, we stress that LD (a statistical association) is distinct from physical linkage (a physical feature of genomes). Consequently, the effects of indirect selection may extend beyond loci physically linked to those under direct selection.
Studies using whole genomes from many individuals and experimental approaches are thus now necessary to understand the potentially complex interplay of factors that shape genome-wide evolution. Here, we determine the genome-wide effects of selection by analysing 491 whole-genome sequences from experimental populations of stick insects. We specifically test two competing hypotheses: (i) indirect selection caused by LD with causal variants is limited to a few regions of the genome, and (ii) indirect selection caused by LD with causal variants is widespread across the genome.
The basis for the relatively untested widespread indirect selection hypothesis is that, especially in small-or moderate-sized populations experiencing polygenic selection, sufficient stochastic LD arises by drift for many or even all single nucleotide polymorphisms (SNPs) to exhibit nontrivial LD with at least one causal variant under selection ( Figure 1). If this occurs, direct selection will be indirectly transmitted to many SNPs across the genome via LD. In reality, the hypotheses posed here form a continuum where a larger number of indirectly selected loci lend greater support to the widespread indirect selection hypothesis. We focus on the Timema stick insect system, but show the generality of our results by additional re-analysis of field transplant experiments in mice and stickleback fish (Barrett et al., 2019;Schluter et al., 2021).

| Study system and field release-recapture experiment
Our focal organism is the wingless, herbivorous stick insect Timema cristinae. This species is endemic to southern California and exhibits green, green-striped and melanistic body colour morphs (Figure 2a).
This system has two key features that make it amenable to testing hypotheses concerning the genomic impact of selection. First, the three morphs are genetically controlled and known to be subject to strong natural selection (Nosil, 2004;Nosil et al., 2018). This includes selection for crypsis on the leaves and stems of the hostplant Adenostoma fasciculatum and Ceanothus spinosus, and negative frequency-dependent selection. Moreover, the 10 megabase (Mbp) region on linkage group 8 (LG8) that controls colour (Mel-Stripe locus hereafter) exhibits heterosis . Thus, past evidence for selection sets up a clear expectation for selection on Mel-Stripe, predicted to impact the genetic loci in LD with this locus. There is potential at least for selection on other traits and genomic regions as well, given past evidence for population divergence in body size, host preference, cuticular hydrocarbons and subtle genome-wide differentiation (Nosil, 2007;Riesch et al., 2017). Second, controlled mark-recapture experiments are possible to estimate survival and thus the effects of selection on Mel-Stripe and genome wide.
We leverage such an experiment here, in which past work reported evidence for selection on genetic regions associated with colour and also other loci, with details that varied with elevation (i.e. microclimate) differences among the transplant sites . This past work thus provides evidence for selection but considered only a small subset of the genome using genotypingby-sequencing (GBS) data and lacked a well-assembled reference genome. Consequently, how much of the genome was indirectly affected by selection because of LD and the extent to which LD occurred between linked vs. unlinked loci is unknown. Here, we use whole-genome sequences to estimate indirect selection in the experiment. Experimental details are as follows. In 2011, we collected 500 T. cristinae from a single population (4.51753°N, 119.80125°W, named FHA), where the dominant host is Adenostoma and all individuals were collected from this host . This site harbours at least thousands of T. cristinae (P. Nosil, personal observation).
A tissue sample was taken from each insect. These individuals were then deployed in a mark-recapture experiment, with 50 insects transplanted onto each of 10 experimental plants (five Adenostoma and five Ceanothus). Our experiment focused on short-term survival in adult and penultimate instar Timema, which tend to live for ~1-3 weeks in the field  and on which predation by birds is a major source of selective mortality (Nosil, 2004;Sandoval, 1994aSandoval, ,1994b. Thus, after 8 days, we recaptured surviving stick insects to quantify genome-wide allele frequency change between release and recapture, and thus test for selection. Notably, past work has shown that recapture is a good proxy for survival, due to the limited dispersal ability of these wingless insects (Nosil, 2004).
The results presented in the main text focus on ~7 million single nucleotide polymorphism (SNPs) from whole-genome sequence data from 246 T. cristinae transplanted to Adenostoma, which is their natal host (Nosil et al., 2018). This approach focuses on survival in a natural habitat, as exemplified by studies of Darwin's finches and other birds (Grant & Grant, 2014;Sepil et al., 2013), and thus avoids effects of transplantation to a novel host. The results in the main text are thus based on allele-frequency change and LD in this experimental population of 246 individuals, which is representative of deme sizes of T. cristinae in nature (Farkas et al., 2013), but also sufficiently finite to allow for stochastic LD to arise among unlinked loci via random sampling and pairing of gametes. Thus, our design realistically embodied how selection acts in this species, and others living in heterogeneous or patchy environments. We present complementary and largely parallel results based on 245 T. cristinae transplanted to Ceanothus, which is their non-natal host, in the Supporting information (see 'Supplementary Analyses and Results' and Figures S1-S3).

| DNA sequence alignment, variant calling and genotype estimation for the experiment
Whole-genome DNA sequence data for individuals from this T. cristinae experiment were described in (Nosil et al., 2018;Riesch et al., 2017), where data were successfully obtained for 491 of the 500 transplanted insects. For the current study, we aligned the F I G U R E 1 Schematic illustration of direct vs. indirect selection. Panel (a) shows how sampling error (i.e. genetic drift) can generate long-range linkage disequilibrium (LD). Here, random sampling from a gene pool results in some multilocus genotype combinations being over-represented (i.e. statistically associated) in a population. Boxes denote the locations of three genetic loci, with filled-in boxes indicating the allele favoured by selection (red) or in LD with the favoured, causal variant (orange). Panel (b) illustrates how direct selection and LD cause indirect selection. A population is shown before and after selection (and, for simplicity, before any reproduction). The 'X's thus denote individuals that died. Selective mortality results in fixation of the favoured allele, and also an increase in the frequency of the associated, functionally neutral alleles. In this example, the effects of association include one allele on a different chromosome than the causal variant [Colour figure can be viewed at wileyonlinelibrary.com] whole-genome DNA sequence data from each of these 491 T. cristinae to the T. cristinae reference genome (version 1.3) using the bwa (version 07.10-r789) mem algorithm with a band width of 100, a 20 bp seed length and a minimum score for output of 30 (Li & Durbin, 2010). We then used samtools (version 1.5) to compress, sort and index the alignments, and to remove PCR duplicates (Li et al., 2009).
We then used the GatK HaplotypeCaller and GenotypeGVCFs modules (version 3.5) to call variants and calculate genotype likelihoods (McKenna et al., 2010). We required a minimum base quality of 30, set the prior probability of heterozygosity to 0.001 and only called variants with a minimum phred-scaled confidence of 50.
The following filters were then applied using custom Perl scripts: minimum coverage of 1× per individual (i.e. 491× coverage across all individuals), a minimum ratio of variant confidence to nonreference read depth of 2, a minimum mapping quality of 40, a maximum phred-scaled p-value of Fisher's exact test for strand bias of 60 and a minimum minor allele frequency of 0.01. Further, we only retained SNPs mapped to one of the 13 T. cristinae linkage groups. This resulted in 7,243,463 SNPs, which were used in subsequent analyses.
Next, we obtained maximum likelihood estimates of allele frequencies for all experimental samples using an expectationmaximization (EM) algorithm, as described in (Li, 2011) and implemented in our own C++ program, estpem (version 0.1) . For this, we used a convergence tolerance of 0.001 and allowed for a maximum of 30 EM iterations. We then used these allele frequency estimates and the genotype likelihoods from GatK to calculate empirical Bayesian genotype estimates. In particular, we computed Pr(g ij = k|l ijk , p i ) ∝ l ijk Pr(g ij = k|p i ) where g ij ∈ (0, 1, 2) denotes the genotype (i.e. count of nonreference alleles) for SNP i and stick insect j, l ijk is the likelihood of genotype k, p i is the nonreference allele frequency, and Pr(g ij = k|p i ) is given by the binomial probability mass function with n = 2. Point estimates (posterior means) were then obtained as ĝ ij = ∑ k∈(0,1,2) k × Pr(g ij = k�l ijk , p i ). These point estimates range from zero to two and are not constrained to be integer values. Lastly, we estimated nucleotide diversity ( ) as an additional metric of genetic variability in the T. cristinae population; this was done using anGsd (0.933-71-g604e1a4) with estimates based directly on the bam alignment files without using the called SNPs (see 'Supplementary Analyses and Results' in the Supporting information for details) (Korneliussen et al., 2014).

Mel-Stripe
We begin with an illustrative example of the potential for indirect selection focused on the Mel-Stripe locus, which exhibits multiple previous lines of evidence for selection including in the specific experiment analysed here Lindtke et al., 2017;Nosil et al., 2018). This serves as a proof-of-concept of indirect selection using a single directly selected locus on which we subsequently expand to consider genome-wide polygenic selection. In F I G U R E 2 Timema cristinae morphs and empirical evidence of selection in the transplant experiment. Panel (a) depicts green unstriped, green striped and melanistic morphs, which are determined by the Mel-Stripe locus. Panel (b) shows linkage disequilibrium (LD) (measured as r 2 ) between genome-wide SNPs and the Mel-Stripe locus. Only the 1% of SNPs exhibiting the greatest LD are shown. Colours denote whether SNPs fall within (red) or outside of (dark grey) the Mel-Stripe locus [Colour figure can be viewed at wileyonlinelibrary.com] this context, we focus on whether substantial LD with Mel-Stripe is limited to SNPs near the locus on linkage group 8 (LG8), or extends to distant SNPs including those on other chromosomes. Following (Nosil et al., 2018), we defined the Mel-Stripe colour locus as spanning ~4.1 Mbps on scaffold 702.1 and ~6.4 Mbps on the adjacent scaffold 128, both on LG8. This is a large, complex structural variant that includes a ~10 Mbp putative inversion and ~1 Mbp deletion Villoutreix et al., 2020). We thus estimated LD as the coefficient of determination between Mel-Stripe genotype (determined from a principal components analysis of SNPs within the Mel-Stripe locus) and genotype estimates for each of ~7 million genome-wide SNPs.
As in past work , we inferred each individual's Mel-Stripe genotype from patterns of clustering in principle components analysis (PCA) space (Li & Ralph, 2019). We first performed a PCA for all 491 T. cristinae based on 79,014 SNPs within the Mel-Stripe locus. This was done on the centred, but not standardized genotype matrix in R (version 3.5.1). We then used k-means clustering to group individuals into clusters based on their PC scores for the first PC axis; we assumed three groups. We delineated Mel-Stripe genotypes based on the k-means cluster (i.e. group) assignment. We then measured pairwise LD between each of the ~7 million SNPs and the Mel-Stripe locus as the square of the correlation between genotypes (Sved, 2009). We wrote a computer program in C++ to efficiently perform this calculation for each SNP (see https:// github.com/zgomp ert/Timem aPoly genic Selec tion.git). The program relies heavily on functions from the Gnu Scientific Library (GSL) (Galassi et al., 2009).

| LD between genome-wide SNPs and Agouti in a mouse experiment or a fitness QTL in a stickleback experiment
To complement our primary analyses, we tested for genome-wide LD between SNPs and a causal variant in two additional, independent data sets. We did this to verify that our general findings were not unique to our specific stick insect data set.
First, we used genomic data from a field enclosure experiment involving cryptically coloured populations of deer mice (Peromyscus maniculatus) transplanted to light-or dark-soil enclosures (Barrett et al., 2019). In this experiment, 481 deer mice were collected and released in light-soil (N = 217) or dark-soil (N = 229) field-enclosures where they were subject to natural selection. Previous analyses of this experiment combined with genetic manipulations showed that a deletion (ΔSer) within the Agouti gene affected the colour of these mice and causally affected fitness (i.e. was subject to direct selection). Here, we tested for LD (measured by r 2 ) between this causal variant and 53,507 genome-wide SNPs that were sequenced for the same experimental individuals. We estimated LD separately for the light-soil and dark-soil experimental populations. This was done using the same analytical approach and software as was used for the analysis of the Mel-Stripe locus in T. cristinae.
Second, we tested for LD between SNP markers and a QTL marker for fitness in an experimental field population of threespine stickleback (Gasterosteus aculeatus) (Schluter et al., 2020(Schluter et al., , 2021. In this experiment, a single F0 intercross was made between a 'marine ecotype' female a 'freshwater ecotype' male. Six F1 females and four F1 males were then crossed (with two males mated twice) to form the F2 generation. Fitness was then measured as the number of surviving F3 offspring attributable to each 632 F2 female fish released in an experimental pond. The experiment uncovered a single QTL for fitness (offspring number), with the peak signal for a SNP on chromosome 4 near the Eda locus (position 12,815,024), a gene known to affect bony armour and fitness (Barrett et al., 2008;Colosimo et al., 2005;O'Brown et al., 2015). We treated this SNP marker as the causal locus for the purpose of our analysis and tested for LD with each of 387 additional genome-wide SNPs in the set of 224 F2 females with genotype data. This was done using the same analytical approach used for T. cristinae and the deer mouse data set.

| Estimating total selection
Hereafter, we return our focus to T. cristinae stick insects. Having documented patterns of LD that would lead to indirect selection (see T. cristinae analyses described above and the Results for details), we next quantified the total extent of selection across the genome during the T. cristinae field experiment. To do this, we used a Bayesian model to estimate the total strength of selection (i.e. direct plus indirect) experienced by each of the ~7 million SNPs during the T. cristinae experiment. Here, total selection reflects the combined effects of direct and indirect selection as captured by differences in marginal fitness (expected fitness averaged across genetic backgrounds) of genotypes at each SNP locus. We assume most SNPs do not causally effect fitness in this short-term experiment (but see Boyle et al., 2017) and thus equate total selection with indirect selection for simplicity. We recognize that this might not be true in a strict sense, but with ~7 million SNPs the average contribution of each SNP to 8-day survival is necessarily small (i.e. on average, each can explain no more than 1 7,000,000 th of the variation in survival), and thus, this assumption should be at least approximately true. We return to this issue of direct vs. indirect selection in the Discussion. We used a Bayesian approach to estimate total selection on each SNP while accounting for genetic drift. Information for this inference comes from mortality in the experiment and patterns of allele frequency change.
We first describe the model and then document details pertaining to how the model was fit.
We assume viability selection such that each individual's absolute fitness (i.e. survival probability) is given by w j ∈ (0, 1). Here, we treat each SNP locus in a separate analysis and assume absolute fitness is a function of an individual's genotype at a SNP and experimental block (i.e. experimental bush). Thus, we allow survival to depend on genotype and also to vary in space. We consider each SNP's association with survival independently with the expectation that any such association mostly arises from LD with (unknown) causal variants not a direct affect on fitness.
We model w as Here, the numerator denotes the expected fitness of individual j based on its genotype at locus i and the selection coefficient (differential) s. The denominator converts this to an absolute fitness by dividing by the product of the sum of the relative fitness values for all individuals in block m ij (i.e. the same block as individual j) and the survival proportion for block m ij (denoted v mj ). This ensures that the mean absolute fitness in each block is equal to the survival proportion for that block. We then assume that y j ∼ Bernoulli(w ij ), where y j is a binary outcome variable denoting whether the individual lived (y j = 1) or died (y j = 0). We assume a beta prior on the strength of total selection experienced by each SNP, such that s i ∼ beta( = 1, = 3). This

| Quantifying LD between SNPs and expected fitness
We hypothesized that the genomically widespread indirect selection documented by the analysis described in the preceding section (see Results for details) was not solely caused by LD with Mel-Stripe, but rather also by LD with other, unknown causal variants, that is by indirect selection arising from oligogenic or polygenic selection.
To test this hypothesis, we asked whether LD between SNPs and genomic predictions of expected fitness (i.e., survival probability) predicted the total selection experienced by each SNP in general, and more so than did LD with Mel-Stripe. The analytical procedure we used to address this question involved three steps: (1) estimate expected fitness of each individual stick insect using genomic prediction and a model of oligogenic/polygenic selection, (2) quantify how strongly each of the ~7 million SNPs was statistically associated with expected fitness at the organism level and (3) test whether the strength of the expected fitness relationship from (2) (an estimate of indirect selection) predicts independent estimates of total selection on each SNP (as described in the previous section). Such predictive power would be expected if our estimates of total selection mostly reflect indirect selection, and thus, this procedure also serves as a semi-independent validation of our estimates of total selection above.
We generated genomic predictions of expected fitness using approach is similar to polygenic genome-wide association methods commonly used for genomic prediction of phenotypes Zhou et al., 2013). Our specific method is tailored to our data and thus has the added benefits of allowing for variation in absolute fitness among transplant sites, directly modelling the experimental design and treating survival explicitly as a binary phenotype. We first describe the model and then second our procedure for fitting the model.
We assume viability selection such that each individual's absolute fitness (i.e. survival probability) is given by w j ∈ (0, 1). Absolute fitness is a function of an individual's multilocus genotype and experimental block (i.e. experimental bush). Thus, we allow survival to depend on genotype but also to vary in space. We model w as follows. We specifically assume multiplicative fitness such that, Here, the numerator denotes the expected fitness of individual j based on its multilocus genotype and the vector of selection coefficients s. The denominator converts this to an absolute fitness by dividing by the product of the sum of the relative fitness values for all individuals in block m j (i.e. the same block as individual j) and the survival proportion for block m j (denoted v m j ). This ensures that the mean absolute fitness in each block is equal to the survival proportion for that block. We then assume that y j ∼ Bernoulli(w j ), where y j is a binary outcome variable denoting whether the individual lived (y j = 1) or died (y j = 0).
We specify a hierarchical model for the selection coefficients.
We first specify a spike-and-slab prior for the number of causal vari- (1) analysis' in the Supporting information). Next, we sample a value for the expected effect size for each casual variant from beta( , ), with α = 1.0 and β = 9.0 (this distribution has an expected value of 0.1).
Effects (selection coefficients) for individual causal variants are then sampled from an exponential distribution with the mean effect set as the expectation; the sign of individual effects is specified at random and with equal probabilities for the nonreference allele increasing or decreasing fitness.
Our primary goal with this model is to use an ABC approach

| Cross-validation to test performance of the genomic prediction model
We used cross-validation to quantify the predictive power of the genomic prediction model described above given our data set. To do this, we ran an additional 400 million simulations with the computer program described in the preceding section. However, in each simulation 80% of individuals were assigned to a training set, with the other 20% assigned to a test set. The training set was used to fit the model, and the test set allowed us to test the model on data not used for model fitting. Assignments to test vs. training sets were random and differed across simulations.
We chose the set of samples to retain to form the posterior as before (i.e., based on the prediction error), but based only on the discrepancy between observed and simulated survival for the 80% of individuals in the training set (i.e., model fit was based solely on the training set). We then measured predictive performance based on the discrepancy between observed and simulated survival for the 20% of individuals in the test set, and considering only the subset of 1000 simulation outputs retained as part of the posterior (based on the training set). We then compared this to the expected error (i.e., discrepancy between observed and simulated) under a null model with no predictive power. To do this, we computed the test set discrepancy for random sets of 1000 simulations (i.e., these were not selected based on the match of the training set to the observed data). We repeated this procedure 100 times to obtain a null distribution of test set errors expected from a model with no actual predictive power.

| Consequences of indirect selection compared to random mortality
We asked whether the patterns we have attributed to indirect selection (see Results) were distinct from patterns of short-term genomic change than would be expected with random mortality (pure drift).
This was done to both assess the genomic consequences of indirect selection and to determine whether our results were expected under a null model of random mortality. We did this by generating and analysing 10 data sets identical to our own, but with survival randomized across individuals within each experimental block (host plant) (we limited our analyses to 10 randomized data sets because of the very large computational burden of running these analyses). Permutations were only done within bushes (blocks) to retain bush-specific mortality rates. We then conducted analyses of selection during the experiment for each of the 10 random mortality data sets exactly as described for the observed data. Thus, these analyses used the same genotypic data, but simply differed in the specific individuals that lived and died.

| Measuring LD and selection in windows
Results from the analyses described in the preceding sections showed LD among genome-wide SNPs likely resulted in SNPs across the genome being impacted by selection (see Results for details), which necessarily makes the short-term direction of evolutionary change for these SNPs more predictable. However, long-range LD is readily broken down by the independent assortment of chromosomes and recombination (i.e. with independent assortment, LD decays by half each generation). Thus, long-range LD is not expected to predict evolution over longer time scales. In contrast, LD among physically linked loci (local LD hereafter) can be more persistent and thus lead to predictable patterns of indirect selection over time. To capture local (rather than long range) patterns of LD and selection, we calculated the average LD and total (mostly indirect) selection in contiguous regions (windows) of the genome. This allowed us to ask whether local LD predicted the intensity of indirect selection during the experiment.
For LD, we calculated the mean LD (measured by the squared genotypic correlation, r 2 ) for all pairs of SNPs in 10 kilo base pair (kbp) windows along each genome scaffold (i.e. windows never spanned multiple scaffolds). This was done using a computer program written in C++ (see https://github.com/zgomp ert/Timem aPoly genic Selec tion.git). For total selection, we instead calculated the mean log-odds of selection in 50 SNP windows. We chose to use windows defined by the number of SNPs here as we found that estimates for physical distance windows were very sensitive to SNP density (in contrast, because LD decays with physically distance, using windows based on physical distance for average LD was necessary). We computed the window-level summaries of total selection in r (version 3.5.1) from the estimates of total selection on SNPs obtained from our total selection model described above.
We fit discrete state, homogeneous hidden Markov models (HMMs) to delineate large-scale, contiguous regions of the genome with elevated LD or total selection (as in Soria-Carrasco et al., 2014).
We used our window-based summaries of LD and total selection (see previous paragraph) as input for the HMM analyses. For LD, we modelled logit mean LD (mean r 2 ) for each window assuming two hidden states, a background LD state (mean set to the empirical mean, µ = −4.49 on the logit scale) and an elevated LD state (mean set to the 99th percentile of the empirical distribution µ = −2.62).
We assumed a Gaussian error distribution in each case with a standard deviation equal to the empirical standard deviation. The Baum-Welch algorithm was used to estimate the transition rate matrix between states (Baum et al., 1970), and the Viterbi algorithm was then used to predict the most likely sequence of hidden states. We modified the r package HiddenMarkov (version 1.8.11) to allow for fixed means and standard deviations for each hidden state and then used the modified version of this package for our analysis (Harte, 2017). We fit the model with a tolerance of 1e −4 and allowing for a maximum of 500 iterations (models converged before this maximum was reached).
We fit the HMM for total selection similarly. We fit the model based on the mean log-odds selection for each window. We again used the empirical mean as the mean for the background state (µ = −0.576), but instead defined the high state as having a mean 0, which corresponds to more evidence for selection than against it when averaged across all SNPs in a window. In each case, we used the empirical standard deviation ( = 0.558). We again used our modified version of the HiddenMarkov package and again fit the model with a tolerance of 1e −4 and a maximum of 500 iterations (convergence was achieved prior to this cutoff). We then tested for an association between window-based estimates of local-LD and total selection.

| Many genetic regions are affected by indirect selection via the Mel-Stripe locus
We LG8 showed LD > 0.1, and LD was small but decidedly nonzero (i.e.

| Generality of results
To verify that such patterns were not unique to our specific stick insect data set, we tested for genome-wide LD between SNPs and a putative causal variant in two additional, independent data sets from two other organisms, deer mice and stickleback fish. We found evidence that SNPs across the genome were associated with the ΔSer genotype in cryptically colored populations of deer mice (Peromyscus maniculatus) transplanted to light-or dark-soil enclosures (Figure 3a,b). LD exceeded ~0.04 for the 1% of SNPs most associated with ΔSer in each population, and these SNPs were distributed across all 23 autosomes and the X chromosome. Moreover, about 20% of SNPs were in LD with ΔSer to a nontrivial extent (i.e.

| The genome is affected by indirect selection because of an independent from Mel-Stripe
Having documented patterns of LD that would lead to indirect selection, we next quantified the total extent of selection across the genome during the T. cristinae field experiment. Our estimates of total selection revealed that different alleles at loci were on average associated with a 11% difference in relative survival probabilities (i.e. a mean difference in marginal fitness or selection differential of 11%, i.e. |ŝ| =0.11, 1st quartile =0.049, 3rd quartile =0.17) ( Figure S5). Notably, because of genetic drift, there was substantial uncertainty in these estimates (i.e. 90% credible in-

| Genome-wide LD predicts indirect selection
We found that the strength of association between SNP genotypes and genomic predictions of expected fitness was strongly explained by our estimates of total selection (Pearson r for log-odds of selection = 0.447, 95% CIs = 0.447-0.448, p < .0001) (Figure 4) (Figure 4b).

| Consequences of indirect selection compared to random mortality
We found that the magnitude of genome-wide allele frequency change for the observed and randomized data sets were very similar, with a mean absolute change of ~0.009 and 95th percentile change of ~0.026 in all cases. Likewise, the distribution of inferred total selection differentials were nearly identical (i.e. |ŝ| = 0.11 for the observed and each randomized data set, with credible evidence of selection impacting ~130,000 SNPs in each data set).
However, the observed and randomized data sets differed in a notable way. Specifically, the extent to which estimates of total selection were correlated with the strength of statistical association between each SNP and genomic predictions of expected fitness was higher for the observed data (r = .447 and .466 for the Adenostoma and Ceanothus host treatments, respectively) than for any of the randomized data sets (r = .249, minimum r = .214, maximum r = .316, P = 2 11 × 1 10 = 0.018). Consequently, these associations with genomic predictions of expected fitness explained 22% of the variation in total selection for the observed data, but only 4%-9% of the variation in the random mortality data sets. This provides evidence that the documented change reflects a contribution from indirect selection, even if selection did not increase the magnitude of shortterm genomic change compared to drift. Because selection adds a deterministic component to evolution, this necessarily increases the extent to which short-term evolutionary change is predictable. In contrast to random drift, this includes prediction of the direction of allele frequency change at each SNP (i.e., neutral alleles positively associated with beneficial alleles or expected fitness can be expected to increase in frequency).

| Physical linkage and effects of local LD in the experiment
Using a HMM approach, we delimited numerous windows with elevated LD relative to the genomic background and did the same for elevated total selection (Figure 5a,b). As predicted by effects of local LD, genomic regions (windows) with higher local LD at the onset of the experiment exhibited greater evidence of selection (Pearson correlation between logit mean r 2 LD and the mean log-odds of selection for each window, r = .23, CI = 0.23-0.24, p < .0001, n = 143,800 windows) ( Figure S6). Notably, the strength of this correlation was greater than for any of the 10 data sets with randomized mortality (mean Pearson correlation r = .18, maximum = 0.21, p = 1/11 = .09).
Thus, both local and long-range LD contributed to indirect selection during the T. cristinae experiment.

| DISCUSS ION
Understanding the impact of selection on the genome is important for understanding adaptation, speciation, genome evolution and eco-evolutionary dynamics. Our results, and related work on linked selection discussed in the introduction (Begun & Aquadro, 1992;Cai et al., 2009;Campos et al., 2017;Elyashiv et al., 2016;Hahn, 2008;Jensen et al., 2008;McGaugh et al., 2012;Pouyet et al., 2018;Roselius et al., 2005;Sella et al., 2009), show that functionally neutral loci in many regions of the genome can be affected by selection.
Specifically, many of the ~7 million SNPs we analysed were likely affected by indirect selection during our relatively brief experiment.
Indeed, our results demonstrate indirect selection affecting SNPs on all 13 T. cristinae chromosomes even when considering only selection on Mel-Stripe, a major locus known to affect colour pattern and fitness (Nosil et al., 2018;Villoutreix et al., 2020). Such widespread indirect selection is not unique to the T. cristinae system, but is a near mathematical certainty in any finite population subject to selection.
Indeed, we conducted additional new genomic analyses of field experiments involving field mice and stickleback fish that demonstrate widespread indirect selection similar to that seen in the stick insects.
Thus, widespread indirect selection is not limited to linked loci and the long time scales of molecular evolution (e.g. Cai et al., 2009;Elyashiv et al., 2016;Sella et al., 2009), but is also prevalent among unlinked loci and on very short time scales, operating even within a single generation.
Nonetheless, our results also show that genetic drift caused by random mortality can mirror many of the genomic patterns detected in our analyses. For example, simulated random mortality caused nearly identical average allele frequency change and associated estimates of total selection. Perhaps this is not surprising, as the intensity of genomic change with ~70% mortality (as occurred in our experiment) should be similar regardless of the specific individuals that live or die, and genotypes at many loci will be statistically associated with realized fitness (survival) even when mortality is random. This conclusion assumes that overall mortality is not higher when selection occurs (i.e. that selection is 'soft') (Wade, 1985;Wallace, 1975). If instead selection results in increased mortality (not just genotype-specific mortality), total genomic change and estimates of selection should be higher than they would be in the absence of selection (Charlesworth, 2009). Moreover, we have convincing evidence of phenotypic selection in the T. cristinae experiment, and consequently, the genomic change we documented necessarily includes indirect selection (i.e. mortality was not random but depended on colour pattern) , and our results did differ from expectations based on random mortality in a few ways. Specifically, we found (i) an increased association between estimates of total selection and the extent to which SNPs were statistically associated with genomic predictions of expected fitness, and (ii) an enhanced association between local LD and genomic windows of elevated total selection. These findings suggest that it is easier to explain or predict the targets of indirect selection than to explain or predict which loci will change the most by drift alone. Still, fully realizing this potential to predict contemporary evolutionary change is not a trivial task, but rather would require a good understanding of the nature of selection, patterns of LD and the genetic basis of relevant trait variation (Nosil, Flaxman, et al., 2020).

F I G U R E 4
Evidence of widespread indirect selection. (a) Conceptual overview and illustration of our method for connecting LD with indirect/total selection. The schematic shows that causal variants determine each organism's expected fitness (e.g. survival probabilities), and statistical association between SNP genotypes and organismal expected fitness determine indirect selection based on differences in the marginal fitness values of (potentially neutral) SNP genotypes. Example plots on the right illustrate our method. First, we estimate expected fitness from our polygenic, genomic prediction model. We then estimate the association between each SNP genotype on expected fitness using linear regression (bottom plot). The r 2 values for each SNP are then used as the independent variables in the next analysis where we ask whether they explain variation in the total selection experienced by each SNP. Panel (b) gives the distribution of error rates for genomic prediction of individual survival over a series of null models (grey region), which is compared to the error rate for the fitted model (orange vertical line). The fitted model outperformed the null models (p <.001). Panel (C) shows that SNPs more strongly associated with expected fitness have a higher posterior probability of having experienced (indirect) selection. The heatmap depicts the relationship across the ~7 million SNPs with individual black points shown in the regions of lowest density. A best fit line is shown (Pearson correlation = 0.45, r 2 = .2, p < .001). The dashed line corresponds with an average posterior probability of 0.5 for selection affecting SNPs in a window [Colour figure can be viewed at wileyonlinelibrary.com] A few additional caveats and limitations concerning our analyses and results should be noted. First, specific details of the systems we considered might make them especially prone to indirect selection.
Specifically, LD was likely elevated in the P. maniculatus and G. aculeatus experimental populations because individuals were collected across an environmental gradient and generated from crosses between divergent ecotypes, respectively (Barrett et al., 2019;Schluter et al., 2020). Moreover, while our results did not rely on specific estimates of selection from these experiments, selection was likely enhanced as the experimental environments were novel or stressful for at least some individuals. However, neither of these concerns applies to the T. cristinae experiment analysed in the main text, which involved a same-host transplant. Second, the first set of analyses we presented for all three systems concerned indirect selection caused by LD with a locus that has a major effect on an ecologically important trait; the presence of a major effect locus affecting fitness is unlikely to be representative of many biological systems. Nonetheless, our results from T. cristinae showed that much more (22×) of the variation in total selection across loci was explained by statistical associations with genomic predictions of expected fitness (and thus LD with many, unknown causal loci) than by LD with the major locus  Thurman & Barrett, 2016). Thus, we cannot rule out a contribution of direct selection to total selection at any specific locus (e.g. Boyle et al., 2017). Still, the contribution of direct selection for most loci is necessarily small, at least in terms of additive effects on fitness, as the average contribution of a locus to the variation in survival over the 8-day experiment cannot be greater than the reciprocal of the number of loci (i.e. ∼ 1 7,000,000 th of the variation in fitness each), which is considerably smaller than the estimates of total selection. Thus, our approach of neglecting direct selection represents an approximation, but likely a very good one. In contrast, no such limits apply to indirect selection; indeed, coupling in hybrid zones can result in all or most loci in the genome experiencing very strong indirect selection because of LD with barrier loci (e.g. Barton, 1983;Barton & Bengtsson, 1986;Butlin & Smadja, 2018;Flaxman et al., 2013).
With the above caveats and considerations in mind, our results have implications for three general themes in biology: (1) the role of deterministic vs. stochastic processes in evolution, (2) the predictability of evolution and (3) the dynamics of adaptation and speciation. In terms of determinism and chance, indirect selection via long-range LD reflects a combination of each. There is a stochastic component to indirect selection that is similar to the chance role F I G U R E 5 Local LD predicts selection and evolution in natural populations. Panel (a) depicts estimates of selection (total selection, which we equate with indirect selection) across the genome. Points denote 50 SNP window averages of the log posterior odds that each SNP experienced selection. The dashed line corresponds to an odds ratio of 1:1. Windows of high indirect selection based on a hidden Markov model (HMM) are shown in red. Panel (b) gives mean pairwise LD (measured by r 2 ) in 10 kb windows along the genome. Windows of high LD based on HMM are shown in red [Colour figure can be viewed at wileyonlinelibrary.com] mutation plays in theories of hitchhiking and genetic draft (Gillespie, 2000). However, conditional on patterns of LD, indirect selection is a deterministic process. Thus, our results support a role for conditional determinism in short-term evolution and provide a demonstration that evolution is simultaneously stochastic and deterministic.
Concerning predictability, the direction and magnitude of change at functionally neutral alleles can be predicted from LD and direct selection, although local LD is more likely to generate longer-term predictions about genetic differentiation, as demonstrated in flycatchers, butterflies and other taxa (Burri et al., 2015;Chaturvedi et al., 2020). This is because local LD decays much more slowly than long-range genome-wide LD. Still, even LD between loci on different chromosomes only decays by 50% each generation, making predictions across several generations possible provided selection pressures and gene interactions remain consistent.
In terms of adaptation, theory predicts that selection is more efficient and selective sweeps more frequent in large populations.
Our results do not contradict this, but do argue that chance longrange LD will cause more of the genome to be affected by indirect selection in smaller populations. Thus, selection in large populations is focused and efficient, whereas in small populations the effects of selection are more widespread but diffuse. Indeed, one way to view the effects of indirect selection for physically unlinked genes in small populations is that across generations it generates a genome-wide reduction in effective population size that can add to that of genetic drift alone. This is true especially if the direction of selection fluctuates in space or time (as selection is then less likely to rapidly exhaust genetic variation), and this effect should be enhanced if selection is accompanied by increased mortality or reduced fecundity (i.e. when selection is hard and maladaptation reduces population size). In sum, our collective results suggest that indirect selection is not simply a statistical nuisance that frustrates attempts to find individual genes under direct selection. Rather, indirect selection is a component of both short-term and long-term processes that shape genome-wide evolution and generate biodiversity.

ACK N OWLED G EM ENTS
We thank A. Buerkle, C. Nice, T. Parchman, R. Pineau, T. Reimchen and the Gompert and Nosil laboratory groups for discussion and comments on previous versions of the manuscript. R. Barrett and S. Laurent graciously provided access to genomic data for Peromyscus.
The support and resources from the Center for High Performance Computing at the University of Utah are gratefully acknowledged. The work was funded by a grant from the National Science Foundation of the United States (grant no. NSF DEB 1844941); this study is part of a project that has received funding from the European Research Council (ERC), to PN, under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 770826 EE-Dynamics).

AUTH O R CO NTR I B UTI O N S
ZG, JF and PN conceived of the project. PN collected the data. ZG analysed the data. ZG, JF and PN wrote and revised the manuscript.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA A N D CO D E AVA I L A B I LIT Y S TATE M E NT
The T. cristinae DNA sequence data have been archived on NCBIs SRA (PRJNA356801). Phenotypic data and genotype estimates for T. cristinae are available on Dryad (https://doi.org/10.5061/dryad. m905q fv26).