Natural selection and the predictability of evolution in
Timema stick insects

Patrik Nosil, Romain Villoutreix, Clarissa de Carvalho, Timothy Farkas,
Víctor Soria-Carrasco, Jeffrey Feder, Bernard Crespi, Zach Gompert

► To cite this version:

Patrik Nosil, Romain Villoutreix, Clarissa de Carvalho, Timothy Farkas, Víctor Soria-Carrasco, et al.. Natural selection and the predictability of evolution in Timema stick insects. Science, 2018, 359 (6377), pp.765-770. 10.1126/science.aap9125. hal-05200425

HAL Id: hal-05200425

https://hal.science/hal-05200425v1

Submitted on 5 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

The
University
Of
Sheffield.

This is a repository copy of *Natural selection and the predictability of evolution in Timemastick insects*.

White Rose Research Online URL for this paper:
https://eprints.whiterose.ac.uk/127707/

Version: Accepted Version

**Article:**

Nosil, P. orcid.org/0000-0002-8271-9005, Villoutreix, R. orcid.org/0000-0002-1815-3844, de Carvalho, C.F. orcid.org/0000-0002-5780-5454 et al. (5 more authors) (2018) *Natural selection and the predictability of evolution in Timemastick insects*. *Science*, 359 (6377). pp. 765-770.

https://doi.org/10.1126/science.aap9125

**Reuse**

Items deposited in White Rose Research Online are protected by copyright, with all rights reserved unless indicated otherwise. They may be downloaded and/or printed for private study, or other acts as permitted by national copyright laws. The publisher or other rights holders may allow further reproduction and re-use of the full text version. This is indicated by the licence information on the White Rose Research Online record for the item.

**Takedown**

If you consider content in White Rose Research Online to be in breach of UK law, please notify us by emailing eprints@whiterose.ac.uk including the URL of the record and the reason for the withdrawal request.

<figure><img src="image_3.png" /></figure>

White Rose
university consortium
Universities of Leeds, Sheffield & York

eprints@whiterose.ac.uk
https://eprints.whiterose.ac.uk/
1 Natural Selection And The Predictability Of Evolution In *Timema* Stick

2 Insects

3

4 Patrik Nosil¹, Romain Villoutreix¹, Clarissa F. de Carvalho¹, Timothy E. Farkas², Víctor Soria-

5 Carrasco¹, Jeffrey L. Feder³, Bernard J. Crespi⁴, Zach Gompert⁵

6

7 ¹Department of Animal and Plant Sciences, University of Sheffield, Sheffield, S10 2TN, UK

8 ²Department of Ecology and Evolutionary Biology, University of Connecticut, Storrs,

9 Connecticut 06369, USA

10 ³Department of Biological Sciences, University of Notre Dame, Notre Dame, Indiana 46556,

11 USA

12 ⁴Department of Biological Sciences, Simon Fraser University, Burnaby, BC V5A 1S6, Canada

13 ⁵Department of Biology, Utah State University, Logan, Utah 84322, USA

14

15 Predicting evolution remains difficult. We study the evolution of cryptic body coloration
16 and pattern in a stick insect using 25 years of field data, experiments, and genomics. We
17 find that evolution is more difficult to predict when it involves a balance between multiple
18 selective factors and uncertainty in environmental conditions than when it involves
19 feedback loops that cause consistent back and forth fluctuations. Specifically, changes in
20 color morph frequencies are modestly predictable through time (r² = 0.14), and driven by
21 complex selective regimes and yearly fluctuations in climate. In contrast, temporal changes
22 in pattern morph frequencies are highly predictable due to negative frequency-dependent
23 selection (r² = 0.86). For both traits, however, natural selection drives evolution around a
24 dynamic equilibrium, providing some predictability to the process.

25

26 Introduction:

27

28 Evolutionary biology is often portrayed as a descriptive rather than predictive science (1, 2).

29 Nonetheless, the extent to which past evolution predicts future evolution can be quantified by

30 testing how well early subsets of a time series predict subsequent changes. However,

31 predictability in the form of such temporal autocorrelation does not consider the underlying
mechanisms driving evolutionary change, and thus can be inherently low. Considering the

mechanisms of evolution can lead to increased understanding of evolutionary change and its

predictability, and we study such mechanisms here.

Evolutionary predictability is mediated by several factors. First, evolution can be unpredictable because of random processes, such as genetic drift (3). Second, even when evolution occurs by deterministic natural selection, predictive power can be diminished if multiple, complex forms of selection act simultaneously and by uncertainties in how the ecological conditions that affect selection change through time (1, 2, 4-7). For example, negative frequency-dependent selection (NFDS) favoring rare alleles can enhance predictability by causing increases in allele frequency to be followed predictably by decreases (and vice-versa) (8, 9). However, the extent of predictability will depend on how NFDS interacts with other evolutionary processes (e.g., directional selection stemming from climate change), and on whether the ecological conditions that affect these other processes can themselves be predicted. Third, the interaction of genes within their genomic context (i.e., dominance and epistasis) and with the environment (i.e., plasticity) may affect the anticipated trajectory of evolution (10-12). For example, directional selection is expected to produce a predictable evolutionary response only if the traits affected are reasonably heritable, and even then responses can be complex and nuanced (4-7, 10).

Studying the predictability of evolution across different timescales is particularly challenging (1, 2). For example, the immediate impact of natural selection can be readily measured in short-term field studies or experiments. Such studies suggest that strong selection is not uncommon at the scale of one or a few generations (13, 14), especially when new environments are colonized (15-17). However, short-term changes need not translate into long-term directional trends. Rather, evolution across geologic and phylogenetic time scales may be characterized by periods of relative stasis interspersed between occasional bursts of sustained directional change, consistent with Simpson's fossil-record inspired model of 'adaptive zones' (18, 19). Indeed, strong but fluctuating selection can generate a pattern of little change when averaged over longer time periods (14). Field studies that measure patterns of evolution over many years or decades are somewhat intermediate between immediate and geologic time scales and thus have the potential to illuminate how short-term selection relates to longer-term patterns of evolution (1, 20-22).
1 Such studies are rare because they require long-term temporal monitoring that cannot be sped up
2 with more effort. We here analyze the predictability of evolution in a long-term study of the stick
3 insect *Timema cristinae* (Fig. 1), and bolster our inferences using manipulative experiments and
4 genomic analyses.
5
6

Polymorphism in stick insects

8 *T. cristinae* is a univoltine, wingless, plant-feeding stick insect that exhibits three morphs that are
9 cryptic on different plant species or tissues (Figs. 1, 2)(23-27). A green morph bearing a white
10 dorsal stripe is cryptic on the leaves of *Adenostoma fasciculatum*, a green and unstriped morph is
11 cryptic on the leaves of *Ceanothus spinosus*, and a melanistic (i.e., brownish/grey and unstriped)
12 morph is cryptic on the stems of both hosts (but is conspicuous on leaves). Accordingly, the
13 striped morph is common on *Adenostoma*, the unstriped morph is common on *Ceanothus*, and
14 the melanistic morph is found at ~10% frequency on both hosts. We refer to the variation
15 between green (striped plus unstriped) and melanistic individuals as 'color polymorphism' and
16 that between green-striped and green-unstriped individuals as 'pattern polymorphism' (i.e., color
17 and pattern are different 'traits'). These polymorphisms are highly heritable with strong genetic
18 dominance (melanistic body coloration is recessive to green and stripe pattern is recessive to
19 unstriped; details below)(23, 24).
20

20

21 Several processes maintain color and pattern polymorphism (23-27). A balance between
22 divergent selection and gene flow between populations on *Adenostoma* versus *Ceanothus* helps
23 maintain pattern polymorphism (26, 27). However, other factors likely contribute because even
24 areas dominated by one host rarely fix for a single pattern morph (26, 27). The frequency of
25 melanism does not vary markedly among populations such that variation in color is maintained
26 by balancing selection within populations, potentially involving heterozygote advantage and
27 selection that varies with microhabitat (stems versus leaves)(23, 24). Spatial and host-related
28 aspects of evolution for these morphs are thus reasonably well understood. In contrast, whether
29 and how morph frequencies change through time, and if they do so predictably, is unknown.
30

30 We studied temporal dynamics in *T. cristinae* using 25 years of field data from the mountains
31 around Santa Barbara, California (545 locality-by-host-by-year estimates of morph frequency on

7

8

9

10

11

12

13

14

15

16

17

18

19

21

22

23

24

25

26

27

28

29

30

31
1 the basis of 34,383 individuals collected from 1990 to 2017, mean n per locality-by-host-by-year
2 = 63, s.d. = 141, Database S1). We focus our autocorrelation analyses of predictability (28) on
3 the locality HV (an acronym for Hidden Valley), where we collected data over a continuous 18-
4 year period, with no years of missing data (total n = 3470, mean yearly n = 193, s.d. = 268).

5 Results:

6 Temporal change in allele frequencies

7
8 We tested the hypothesis that temporal change in allele frequencies at the genetic region
9 underlying the morphs is due, in part, to natural selection. Although past studies support
10 selection on the morphs (23-27), these results do not mean that all temporal changes are due to
11 selection, because selection and drift are not mutually exclusive (29, 30). Genomic data from
12 different time points provide a means to test for selection, because strongly selected regions are
13 expected to show greater change through time than the more neutral genomic background (29,
14 31). Genomic data are further required in our specific instance because genetic dominance and
15 heterozygote excess complicate inference of allele frequency change using phenotypic data alone
16 (23, 24).
17

18 We used *de novo* genome sequencing of a melanistic and green morph, with Dovetail hi-rise
19 scaffolding of Illumina reads (N50 = ~16 and 8 megabases, respectively)(32), linkage mapping
20 (25), and genome-wide association (GWA) mapping to explicitly delimit a single, contiguous
21 genomic region (~10.5 megabases in size) associated with color and pattern variation (Figs. 2,
22 S1-2). Consistent with a similar study with a more fragmented reference genome, this region
23 exhibits three core haplotypes (i.e., alleles), one corresponding to each morph, designated *s*, *u*,
24 and *m* for green-striped, green-unstriped, and melanistic, respectively (i.e., in terms of diploid
25 genotypes and phenotypes: *uu*, *us* and *um* = green-unstriped; *ss* and *sm* = green-striped; *mm* =
26 melanistic)(23). We refer to this region as the *Mel-Stripe* locus hereafter.
27

28 We quantified allele frequency changes at *Mel-Stripe* over time within three published data sets:
29 (1) genotyping-by-sequencing (GBS) data collected in a natural population on *Adenostoma*
30 (FHA, acronym for locality Far Hill on *Adenostoma*) in 2011 and 2013 (30, 33), (2) re-sequenced
1 whole genomes from individuals collected in FHA and used in an eight-day (i.e., within-generation) release and recapture field experiment (30), and (3) GBS data in a between-year (i.e., between-generation) field transplant experiment (25).

4
5 In each case, we contrasted change through time at *Mel-Stripe* to that of the remainder of the
6 genome (to all genomic scaffolds, i.e., loci, that harbored as many single nucleotide
7 polymorphisms as *Mel-Stripe*, which was 40, 16, and 39 loci, respectively, for the data sets noted
8 above). Due to our explicit interest in *Mel-Stripe* we did not attempt to delimit other loci under
9 selection. If such loci exist they could upwardly bias our estimates of genome-wide change
10 relative to a case of neutrality, making our results for *Mel-Stripe* conservative.

11
12 We found that *Mel-Stripe* showed the greatest temporal allele frequency change of all genomic
13 regions, in all three data sets (FHA, change = 0.0273, *P* = 0.024; within-generation experiment,
14 change = 0.0340, *P* = 0.059; between-generation experiment, change = 0.0988, *P* = 0.025; exact
15 probabilities; Fisher's combined probability test across data sets: X² = 20.50, d.f. = 6, *P* = 0.0023,
16 Fig. 2). Dispersal alone is unlikely to drive these observed patterns because FHA was sampled
17 over an area that is larger (>10,000 m²) than the dispersal capacity of *T. cristinae* (i.e., one to a
18 few dozen meters per generation)(30, 33, 34). Furthermore, field surveys detected essentially no
19 dispersal off experimental bushes in the recapture study (30), and selection on pattern has been
20 previously observed in the presence, but not the absence of predation (with dispersal possible in
21 both treatments)(35, 36). Thus, selection likely contributed to the genetic change we observed at
22 the *Mel-Stripe* locus. We thus next turned to whether such selection was associated with weakly
23 or strongly predictable patterns of evolution.

24
25 Predictability of the evolution of body color and complex selection regimes
26

26
27 We quantified the predictability of evolution using autoregressive moving average models,
28 ARMA (28). This analysis revealed that that color morphs at Hidden Valley (HV) exhibited
29 subtle and only moderately predictable changes through time (median predictive r² = 0.14,
30 ARMA, Fig. 3, Table S1). This was associated with support for multiple, complex, and
31 counteracting sources of selection (Fig. 4). Across the 25-year study period, the frequency of
1 melanistic morphs increased in years where spring temperatures were warmer (overall effect of temperature = 0.187, 95% equal-tail probability intervals = 0.063-0.309; correlation between observed and predicted frequency from cross-validation = 0.16, 95% CI = 0.04-0.28, P = 0.0102, Bayesian hierarchical linear model, Fig. 4, Table S2). However, lab experiments indicate that melanistic individuals have weaker heat tolerance, relative to green individuals (B = 3.57, 95% CIs = 1.34-9.51, P = 0.0111, Cox proportional hazards regression model using exact likelihood). These results imply that selection for crypsis on dry, brownish plants in warmer years may favor dark colors, but that thermoregulatory selection acts in an opposing direction. However, further work is required to test this hypothesis directly, and to establish how well the laboratory experiments match field conditions. Notably, melanistic individuals exhibit fewer fungal infections and greater mating success than other morphs, further suggesting that selection is multi-faceted (24).

13
14 Given selection appears complex, we used our genomic data to estimate selection coefficients explicitly (for all six diploid genotypes underlying the morphs, with three alleles: s, m, and u). This analysis revealed that viability selection on *Adenostoma* during late life-history stages (i.e., in the within-generation experiment) favored the *s/s* homozygote (posterior probability that fitness of *s/s* > than the following genotypes: *m/u* = 0.93, *u/u* = 0.81, *m/m* = 0.82, *m/s* = 0.84, *u/s* = 0.91). In contrast, the most-fit genotype between years at the FHA locality (also on *Adenostoma*) was the *s/m* heterozygote (posterior probability that fitness of *s/m* > than the following genotypes: *m/u* = 0.90, *s/s* = 0.97, *u/u* = 0.81, *m/m* = 0.92, *u/s* = 0.94). Both *s/s* and *s/m* are green-striped in terms of phenotype, and cryptic on *Adenostoma*. Thus, a fluctuating balance between many factors, potentially including heterozygote advantage (23), may explain why evolution involving the deterministic process of selection was not more highly predictable.
20
21
22
23
24
25
26 Predictability of the evolution of pattern
27
28 Our findings for color raise the question of whether evolution is ever highly predictable? We address this issue by re-analyzing the data from Hidden Valley (HV) considering striped versus unstriped individuals (i.e., pattern, rather than color, polymorphism). This demonstrates that striped morph frequencies at HV exhibited consistent increases followed by decreases (i.e., up

28
29
30
31
1 and down fluctuations) across 18 consecutive years (Binomial sampling probability < 0.0001,
2 Fig. 3). Thus, the evolution of pattern was highly predictable (median predictive r² = 0.86,
3 ARMA). In fact, predictive power at the scale of three to five years was near perfect (>0.95), and
4 remained high even after a decade (>0.80). The observed pattern of predictable of up and down
5 fluctuations could reflect a case where predators have specific search images for common prey,
6 resulting in NFDS selection favoring rare prey phenotypes.
7

8 Morph frequency and selection

10 To test the NFDS hypothesis, we transplanted green-striped and green-unstriped *T. cristinae* to
11 *Adenostoma* in either 1:4 or 4:1 ratios (n = 1000 individuals, 10 replicates per treatment, Fig. 5).
12 The NFDS hypothesis predicts that the striped morph, cryptic on *Adenostoma*, would exhibit a
13 stronger survival advantage when rare. Supporting this prediction, the striped morph experienced
14 strong selection when initially rare (selection coefficient, *s* = 0.70), and increased in frequency in
15 all 10 experimental replicates (posterior probability that change > 0 was >0.999). In contrast, the
16 striped morph showed idiosyncratic changes when initially common (*s* = -0.04, posterior
17 probability that change > 0 = 0.43). Although our results differ from NFDS where the sign of
18 selection reverses very strongly, they are consistent with the strength of selection being
19 dependent on frequency (i.e., directional selection that weakens with increasing allele frequency
20 is akin to frequency-dependent selection). It is possible that selection against the striped morph
21 would be more strongly negative if ratios were manipulated more extremely (e.g., to 10:1).

22 Conclusions:

23
24 We observed complex and fluctuating sources of selection. Together with gene flow (26), these
25 selection pressures likely contribute to relative stability in the difference between hosts in morph
26 frequency over the 25-year study period (Fig. 5). Complex and fluctuating selection may thus
27 help maintain polymorphism, but prevent divergence of sufficient magnitude to strongly drive
28 speciation (24). NFDS in particular may cause evolutionary systems to exhibit resilience, as
29 reported in other complex ecological, social, and physical systems (37-39). Such resilience
30 increases predictability of short-term evolution, because a system returns to its former state
31 following perturbation. However, it can make long-term predictions difficult because substantial
1 evolution only happens when the system reaches a tipping point that pushes it more permanently
2 to a very different state.

3 Finally, our results suggest that the predictability of evolution can depend on the nature of
4 selection, and our understanding of it. Thus, we predicted evolution more accurately for pattern,
5 where selection appears to be strongly associated with frequency, than for color, where a myriad
6 of factors, some poorly understood, affect fitness. As further illustration of this point,
7 predictability in the form of temporal autocorrelation is modest to weak in other well-known
8 studies of contemporary evolution: beak and body size changes in Darwin's finches and morph
9 frequency changes in the scarlet tiger and peppered moth(21, 40, 41)(median predictive power
10 per study system, r², in 3-10 year forecasting analyses 0.03-0.18, ARMA, Figs 6, S3-S6, Table
11 S1). Adding climatic (i.e., rainfall) data to the *Geospiza* finch case, where climate is known to
12 affect seed distributions, improves predictive power in *Geospiza fortis* (e.g., mean r² increase
13 relative to a model without rainfall = 0.08, ARMA; Table S3). Nonetheless, predictive power
14 even considering rainfall is modest and increased only in one of the two finch species examined
15 (Table S3). It is possible that predictability remains limited because seed size itself was not
16 modeled, because the relationship between evolution and rainfall is complex such that only
17 extreme droughts have strong effects, and because some extreme climatic events precede the ten-
18 year period that our forecasting is based upon.

19
20 In conclusion, our constrained understanding of selection and environmental variation (i.e.,
21 limits on data and analysis), rather than inherent randomness, can thus limit ability to predict
22 evolution. In turn, these limitations may affect our understanding of ecological processes,
23 because to the extent that evolution can be predicted, perhaps so can its ecological consequences
24 for population dynamics, community structure, and ecosystem functioning (42-44).
25
26
27
28 The work was funded by a grant from the European Research Council (NatHisGen R/129639) to
29 P. Nosil and a Natural Sciences and Engineering Research Council grant to B.J. Crespi (NSERC
30 06505). V. Soria-Carrasco was supported by a Leverhulme Trust Early Career Fellowship. We
31 thank C. Sandoval for contributing data on morph frequencies. C. Sandoval, A. Comeault, M.

26 Acknowledgements

27
28 The work was funded by a grant from the European Research Council (NatHisGen R/129639) to
29 P. Nosil and a Natural Sciences and Engineering Research Council grant to B.J. Crespi (NSERC
30 06505). V. Soria-Carrasco was supported by a Leverhulme Trust Early Career Fellowship. We
31 thank C. Sandoval for contributing data on morph frequencies. C. Sandoval, A. Comeault, M.
1 Muschick, R. Riesch, and T. Reimchen provided discussion and comments on previous versions
2 of the manuscript. The support and resources from the Center for High Performance Computing
3 at the University of Utah are gratefully acknowledged, as well as access to the High Performance
4 Computing Facilities, particularly to the Iceberg HPC cluster, from the Corporate Information
5 and Computing Services at the University of Sheffield. The genome drafts reported in this paper
6 have been deposited in the Genomes database of the National Center for Biotechnology
7 Information (NCBI) under BioProject Accession PRJNA417530. Other data, including custom
8 code written for analyses, have been archived in Dryad Digital Repository XXXX (number can
9 only be provided by Dryad if the manuscript is accepted). Rosa Ribas drew all the figures.

10

11 Author contributions

12

13 PN and ZG conceived the project. PN, RV, CC, TEF, VS, ZG contributed data and analyses. PN,
14 JLF, and ZG wrote the manuscript, with feedback from all authors.

15

16 **Figure 1.** Drawings of the three morphs of *T. cristinae*.

17

18 **Figure 2.** Genomic change through time at the *Mel-Stripe* locus. (a) Manhattan plots showing
19 results for genome-wide association mapping of color. The y-axes show *P*-values, with red
20 denoting genome-wide significance. The left-hand plot shows results genome wide (LG = linkage
21 group). The right-hand plot is zoomed in on LG8, which shows the bulk of association, and here
22 numbers below the x-axis delimit different genomic scaffolds. The *Mel-Stripe* locus is evident by
23 the block of strong association spanning scaffolds 702.1 and 128. (b) Allele-frequency change
24 through time in the natural population FHA (2011 versus 2013). (c) Allele-frequency change
25 through time in the within-generation experiment. (d) Allele-frequency change through time in
26 the between-generation experiment. In panels b-d the vertical red line shows change at the *Mel-
27 Stripe* locus and the histogram shows the distribution of change across other similar-sized
28 scaffolds in the genome (i.e., the genomic background).

29

30 **Figure 3.** Predicting evolution in *Timema cristinae* stick insects. (a) Schematic of analytical
31 approach for predicting evolution using temporal autocorrelation. Black points represent
observed values in a time series. Some number of these observed points are removed (in this case the six right-hand most points) and the missing values for them are predicted using the remainder of observed points. Predictive power is the strength of association between observed and predicted values. (b) Color morph frequencies through time. (c) Pattern morph frequencies through time. (d) Change in color morph frequencies. (e) Change in pattern morph frequencies. (f) Predicting change in color morph frequencies ($r^2$). (g) Predicting change in pattern morph frequencies ($r^2$). (h) Predicting change in color morph frequencies (r). The difference from panel (f) is that r-values are not squared such that their sign is evident. Shaded areas are 95% confidence intervals. (i) Predicting change in pattern morph frequencies (r). The difference from panel (g) is that r-values are not squared such that their sign is evident. Shaded areas are 95% confidence intervals.

**Figure 4.** Complex patterns of natural selection. (a) Associations between the frequency of melanistic morphs and yearly spring temperature. Positive effects indicate increases in melanistic frequency with increased temperature. Significant effects are shown in red. The left-most data point represents the average effect across populations, and the remaining points are for individual populations. Note that among-population variation is high, but that all significant effects are positive. (b) Morph-specific survival time in thermoregulatory (i.e., heat tolerance) lab experiments. (c) Genotype specific fitness in the within-generation experiment on the host *Adenostoma* (*s/s* is most fit). Shown are the posterior probabilities from estimates of genotype-specific fitness. (d) Genotype specific fitness in the natural population FHA on *Adenostoma* (*s/m* is most fit). Shown are the posterior probabilities from estimates of genotype-specific fitness.

**Figure 5.** Evidence for negative-frequency dependent selection on pattern, and resulting stability in morph frequency differences between hosts. (a) Posterior probability estimates of the selection coefficient in each treatment. Positive values on the x-axis represent selection favoring striped individuals. (b) Changes in the frequency of striped morphs in releases at 20% of the population and 80% of the population during the course of the experiment. (c) Posterior probability distributions of change in the frequency of striped morph per treatment. (d) Yearly differences between host plant species in natural populations in the frequency of striped morphs (lines denote posterior medians and shaded regions given the 95% equal-tail probability intervals).
1

**Figure 6.** Predicting evolution in different systems. (a) Predicting evolution for evolutionary time series in (left to right) *Geospiza fortis* body size, *Geospiza fortis* beak morphology (principle components 1 and 2), *G. scandens* body size, *G. scandens* beak morphology (principle components 1 and 2), *Panaxia dominula* and *Biston betularia* morph frequencies, and *T. cristinae* color and pattern morph frequencies. Boxplots show the distribution of $r^2$ between true and predicted evolutionary change across 3 to 10 year model-based forecasts. (b) Predicting evolution in studies (r). The difference from panel (a) is that r-values are not squared such that their sign is evident.

10

11 **References**

12

13 1. P. R. Grant, B. R. Grant, Unpredictable evolution in a 30-year study of Darwin's finches.
   *Science* **296**, 707-711 (2002).

15 2. M. Lässig, V. Mustonen, A. M. Walczak, Predicting evolution. *Nature Ecology & Evolution* **1**, 0077 doi:0010.1038/s41559-41017-40077 (2017).

17 3. D. L. Hartl, A. G. Clark, *Principles of population genetics*, fourth edition. (Sinauer, 2007).

19 4. A. Ozgul, S. Tuljapurkar, T. G. Benton, J. M. Pemberton, T. H. Clutton-Brock, T. Coulson, The Dynamics of Phenotypic Change and the Shrinking Sheep of St. Kilda. *Science* **325**, 464-467 (2009).

22 5. T. Coulson, D. R. MacNulty, D. R. Stahler, B. Vonholdt, R. K. Wayne, D. W. Smith, Modeling Effects of Environmental Change on Wolf Population Dynamics, Trait Evolution, and Life History. *Science* **334**, 1275-1278 (2011).

25 6. J. Merila, Evolution in response to climate change: In pursuit of the missing evidence.
   *Bioessays* **34**, 811-818 (2012).

27 7. M. Bosse, L. G. Spurgin, V. N. Laine, E. F. Cole, J. A. Firth, P. Gienapp, A. G. Gosler, K. McMahon, J. Poissant, I. Verhagen, M. A. M. Groenen, K. van Oers, B. C. Sheldon, M. E. Visser, J. Slate, Recent natural selection causes adaptive evolution of an avian polygenic trait. *Science* **358**, 365-368 (2017).
1 8. M. Chouteau, V. Llaurens, F. Piron-Prunier, M. Joron, Polymorphism at a mimicry supergene maintained by opposing frequency-dependent selection pressures. *Proc. Natl. Acad. Sci. U. S. A.* **114**, 8325-8329 (2017).

4 9. D. I. Bolnick, W. E. Stutz, Frequency dependence limits divergent evolution by favouring rare immigrants over residents. *Nature* **546**, 285-+ (2017).

6 10. R. Lande, Quantitative genetic analysis of multivariate evolution, applied to brain-body size allometry. *Evolution* **33**, 402-416 (1979).

8 11. C. Rueffler, T. J. M. Van Dooren, O. Leimar, P. A. Abrams, Disruptive selection and then what? *Trends Ecol. Evol.* **21**, 238-245 (2006).

10 12. X. Thibert-Plante, A. P. Hendry, The consequences of phenotypic plasticity for ecological speciation. *J. Evol. Biol.* **24**, 326-342 (2011).

12 13. J. G. Kingsolver, H. E. Hoekstra, J. M. Hoekstra, D. Berrigan, S. N. Vignieri, C. E. Hill, A. Hoang, P. Gibert, P. Beerli, The strength of phenotypic selection in natural populations. *Am. Nat.* **157**, 245-261 (2001).

15 14. G. Bell, Fluctuating selection: the perpetual renewal of adaptation in variable environments. *Philos. Trans. R. Soc. B-Biol. Sci.* **365**, 87-97 (2010).

17 15. P. F. Colosimo, K. E. Hosemann, S. Balabhadra, G. Villarreal, M. Dickson, J. Grimwood, J. Schmutz, R. M. Myers, D. Schlüter, D. M. Kingsley, Widespread parallel evolution in sticklebacks by repeated fixation of ectodysplasin alleles. *Science* **307**, 1928-1933 (2005).

21 16. S. M. Rogers, P. Tamkee, B. Summers, S. Balabahadra, M. Marks, D. M. Kingsley, D. Schlüter, GENETIC SIGNATURE OF ADAPTIVE PEAK SHIFT IN THREESPINE STICKLEBACK. *Evolution* **66**, 2439-2450 (2012).

24 17. C. Heliconius Genome, Butterfly genome reveals promiscuous exchange of mimicry adaptations among species. *Nature* **487**, 94-98 (2012); published online Epub2012-Jul-5 (26 18. G. G. Simpson, *Tempo and mode in evolution*. (Columbia University Press, New York, 1944).

28 19. G. G. Simpson, *The major features of evolution*. (Columbia University Press, New York, 1953).

30 20. A. M. Siepielski, J. D. DiBattista, S. M. Carlson, It's about time: the temporal dynamics of phenotypic selection in the wild. *Ecol. Lett.* **12**, 1261-1276 (2009).
21. P. R. Grant, B. R. Grant, *40 Years of Evolution: Darwin's Finches on Daphne Major Island*. (Princeton University Press, Princeton, 2014).

22. T. Dobzhansky, *Genetics and the Origin of Species*. (Columbia University Press, New York, NY, ed. 3rd, 1951), pp. 364.

23. D. Lindtke, K. Lucek, V. Soria-Carrasco, R. Villoutreix, T. E. Farkas, R. Riesch, S. R. Dennis, Z. Gompert, P. Nosil, Long-term balancing selection on chromosomal variants associated with crypsis in a stick insect. *Mol. Ecol.* **26**, 6189-6205 (2017).

24. A. A. Comeault, S. M. Flaxman, R. Riesch, E. Curran, V. Soria-Carrasco, Z. Gompert, T. E. Farkas, M. Muschick, T. L. Parchman, T. Schwander, J. Slate, P. Nosil, Selection on a Genetic Polymorphism Counteracts Ecological Speciation in a Stick Insect. *Current Biology* **25**, 1-7 (2015).

25. V. Soria-Carrasco, Z. Gompert, A. A. Comeault, T. E. Farkas, T. L. Parchman, J. S. Johnston, C. A. Buerkle, J. L. Feder, J. Bast, T. Schwander, S. P. Egan, B. J. Crespi, P. Nosil, Stick Insect Genomes Reveal Natural Selection's Role in Parallel Speciation. *Science* **344**, 738-742 (2014).

26. C. P. Sandoval, The effects of relative geographical scales of gene flow and selection on morph frequencies in the walking-stick *Timema cristinae*. *Evolution* **48**, 1866-1879 (1994).

27. P. Nosil, Divergent host plant adaptation and reproductive isolation between ecotypes of *Timema cristinae* walking sticks. *Am. Nat.* **169**, 151-162 (2007).

28. N. Cressie, C. K. Wikle, *Statistics for Spatio-Temporal Data*. (John Wiley and Sons, 2011).

29. Z. Gompert, Bayesian inference of selection in a heterogeneous environment from genetic time-series data. *Mol. Ecol.* **25**, 121-134 (2016).

30. Z. Gompert, A. A. Comeault, T. E. Farkas, J. L. Feder, T. L. Parchman, C. A. Buerkle, P. Nosil, Experimental evidence for ecological selection on genome variation in the wild. *Ecol. Lett.* **17**, 369-379 (2014).

31. M. Foll, H. Shim, J. D. Jensen, WFABC: a Wright-Fisher ABC-based approach for inferring effective population sizes and selection coefficients from time-sampled data. *Molecular ecology resources* **15**, 87-98 (2015).
1 32. N. H. Putnam, B. L. O'Connell, J. C. Stites, B. J. Rice, M. Blanchette, R. Calef, C. J.

2 Troll, A. Fields, P. D. Hartley, C. W. Sugnet, D. Haussler, D. S. Rokhsar, R. E. Green,

3 Chromosome-scale shotgun assembly using an in vitro method for long-range linkage.

4 *Genome Research* **26**, 342-350 (2016).

5 33. R. Riesch, M. Muschick, D. Lindtke, R. Villoutreix, A. A. Comeault, T. E. Farkas, K.

6 Lucek, E. Hellen, V. Soria-Carrasco, S. R. Dennis, C. F. de Carvalho, R. J. Safran, C. P.

7 Sandoval, J. L. Feder, R. Gries, B. J. Crespi, G. Gries, Z. Gompert, P. Nosil, Transitions

8 between phases of genomic differentiation during stick-insect speciation. *Nature Ecology*

9 *and Evolution* **1**, 0082 (2017).

10 34. C. Sandoval, Persistence of a walking-stick population (Phasmatoptera : Timematodea)

11 after a wildfire. *Southwestern Naturalist* **45**, 123-127 (2000).

12 35. P. Nosil, Reproductive isolation caused by visual predation on migrants between

13 divergent environments. *Proc. R. Soc. B-Biol. Sci.* **271**, 1521-1528 (2004).

14 36. P. Nosil, B. J. Crespi, Experimental evidence that predation promotes divergence in

15 adaptive radiation. *Proc. Natl. Acad. Sci. U. S. A.* **103**, 9090-9095 (2006).

16 37. M. Scheffer, S. R. Carpenter, T. M. Lenton, J. Bascompte, W. Brock, V. Dakos, J. van de

17 Koppel, I. A. van de Leemput, S. A. Levin, E. H. van Nes, M. Pascual, J. Vandermeer,

18 Anticipating Critical Transitions. *Science* **338**, 344-348 (2012).

19 38. M. Scheffer, S. Carpenter, J. A. Foley, C. Folke, B. Walker, Catastrophic shifts in

20 ecosystems. *Nature* **413**, 591-596 (2001).

21 39. P. Nosil, J. L. Feder, S. M. Flaxman, Z. Gompert, Tipping points in the dynamics of

22 speciation. *Nature Ecology and Evolution* **1**, 0001 (doi:0010.1038/s41559-41016-40001)

23 (2017).

24 40. L. M. Cook, D. A. Jones, The medionigra gene in the moth Panaxia dominula: The case

25 for selection. *Philos. Trans. R. Soc. Lond. Ser. B-Biol. Sci.* **351**, 1623-1634 (1996).

26 41. L. M. Cook, S. L. Sutton, T. J. Crawford, Melanic moth frequencies in Yorkshire, an old

27 English industrial hot spot. *Journal of Heredity* **96**, 522-528 (2005).

28 42. A. P. Hendry, *Eco-evolutionary dynamics*. (Princeton University Press, Princeton, New

29 Jersey, 2017).

30 43. T. W. Schoener, The Newest Synthesis: Understanding the Interplay of Evolutionary and

31 Ecological Dynamics. *Science* **331**, 426-429 (2011).
1 44. T. E. Farkas, T. Mononen, A. A. Comeault, I. Hanski, P. Nosil, Evolution of Camouflage

2 Drives Rapid Ecological Change in an Insect Community. *Current Biology* **23**, 1835-

3 1843 (2013).

4

5 **Supplementary References**

6

7 45. B. Langmead, S. L. Salzberg, Fast gapped-read alignment with Bowtie 2. *Nature*

8 *Methods* **9**, 357-U354 (2012).

9 46. H. Li, B. Handsaker, A. Wysoker, T. Fennell, J. Ruan, N. Homer, G. Marth, G. Abecasis,

10 R. Durbin, The Sequence Alignment/Map format and SAMtools. *Bioinformatics* **25**,

11 2078-2079 (2009).

12 47. Z. Gompert, L. K. Lucas, C. C. Nice, C. A. Buerkle, GENOME DIVERGENCE AND

13 THE GENETIC ARCHITECTURE OF BARRIERS TO GENE FLOW BETWEEN

14 LYCAEIDES IDAS AND L-MELISSA. *Evolution* **67**, 2498-2514 (2013).

15 48. Y. S. Aulchenko, S. Ripke, A. Isaacs, C. M. Van Duijn, GenABEL: an R library for

16 genome-wide association analysis. *Bioinformatics* **23**, 1294-1296 (2007).

17 49. A. L. Price, N. J. Patterson, R. M. Plenge, M. E. Weinblatt, N. A. Shadick, D. Reich,

18 Principal components analysis corrects for stratification in genome-wide association

19 studies. *Nature Genetics* **38**, 904-909 (2006).

20 50. P. Rastas, F. C. F. Calboli, B. C. Guo, T. Shikano, J. Merila, Construction of Ultradense

21 Linkage Maps with Lep-MAP2: Stickleback F-2 Recombinant Crosses as an Example.

22 *Genome Biology and Evolution* **8**, 78-93 (2016).

23 51. H. Li, R. Durbin, Fast and accurate short read alignment with Burrows-Wheeler

24 transform. *Bioinformatics* **25**, 1754-1760 (2009).

25 52. T. Schwander, B. J. Crespi, Multiple Direct Transitions from Sexual Reproduction to

26 Apomictic Parthenogenesis in Timema Stick Insects. *Evolution* **63**, 84-103 (2009).

27 53. S. V. Angiuoli, S. L. Salzberg, Mugsy: fast multiple alignment of closely related whole

28 genomes. *Bioinformatics* **27**, 334-342 (2011).

29 54. A. McKenna, M. Hanna, E. Banks, A. Sivachenko, K. Cibulskis, A. Kernytsky, K.

30 Garimella, D. Altshuler, S. Gabriel, M. Daly, M. A. DePristo, The Genome Analysis
1 Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data.

2 *Genome Research* **20**, 1297-1303 (2010).

3 55. Z. Gompert, J. P. Jahner, C. F. Scholl, J. S. Wilson, L. K. Lucas, V. Soria-Carrasco, J. A.

4 Fordyce, C. C. Nice, C. A. Buerkle, M. L. Forister, The evolution of novel host use is
5 unlikely to be constrained by trade-offs or a lack of genetic variation. *Mol. Ecol.* **24**,
6 2777-2793 (2015).

7 56. O. J. T. Briet, P. H. Amerasinghe, P. Vounatsou, Generalized Seasonal Autoregressive Integrated Moving Average Models for Count Data with Application to Malaria Time Series with Low Case Numbers. *PLoS One* **8**, (2013).

10 57. M. C. Jones, Randomly choosing parameters from the stationary and invertibility region
11 of autoregressive-moving average models. *Applied Statistics* **36**, 134–138 (1987).

12 58. C. P. Sandoval, Differential Visual Predation on Morphs of Timema-Cristinae (Phasmatodeae, Timemidae) and Its Consequences for Host-Range. *Biol. J. Linnean Soc.* **52**, 341-356 (1994).

15 59. Q. D. Team, **QGIS Geographic Information System. Open Source, Geospatial**
**Foundation Project.** http://www.qgis.org/. (2016).

17 60. T. Therneau, P. Grambsch, *Modeling Survival Data: Extending the Cox Model*. Springer-Verlag (2000).

19 61. T. Therneau. A Package for Survival Analysis in S version, 2.38, <URL:
20 http://CRAN.R-project.org/package=survival>. (2015).

21 62. R. D. C. Team. (Vienna, Austria, 2013).

22 63. Z. Gompert, F. J. Messina, Genomic evidence that resource-based trade-offs limit host-range expansion in a seed beetle. *Evolution* **70**, 1249-1264 (2016).

24

25
<figure><img src="image_4.png" /></figure>

Supplementary Materials for

Natural Selection And The Predictability Of Evolution In *Timema* Stick Insects

Patrik Nosil, Romain Villoutreix, Clarissa F. de Carvalho, Timothy E. Farkas, Víctor Soria-Carrasco, Jeffrey L. Feder, Bernard J. Crespi, Zach Gompert

correspondence to: p.nosil@sheffield.ac.uk and zach.gompert@usu.edu

This PDF file includes:

Materials and Methods

Figs. S1 to S6

Tables S1 to S3

Caption for database S1

Other Supplementary Materials for this manuscript includes the following:

Database S1 [Raw population data]
# Materials and Methods

2

3 **Approach for delimiting the genetic region affecting color and color-pattern**

4

5 Previous work showed that each morph in *T. cristinae* is a chromosomal form underlain by
6 a haplotype on a single linkage group (LG8), with restricted recombination between
7 chromosomal forms (23, 24). However, it relied on a fragmented reference genome such that it
8 could not delimit a single, contiguous region (i.e., locus) underlying each morph. We here
9 delimit the locus underlying the morphs and quantify its change through time relative to the rest
10 of the genome.
11

12

13 To do so, we generated higher-quality reference genomes for a melanistic and a green
14 morph of *T. cristinae* using Dovetail hi-rise scaffolding of Illumina reads (N50 = ~16, 8
15 megabases, respectively)(32). Comparison of the reference genomes, linkage mapping (25), and
16 genome-wide association (GWA) mapping allowed us to explicitly delimit a single, contiguous
17 genomic region associated with color and pattern variation (Figure 2, S1, S2). Accordingly, this
18 region exhibits three core haplotypes (i.e., alleles), one corresponding to each morph (with
19 melanistic recessive to green body coloration and stripe recessive to unstriped pattern), and we
20 refer to it as the *Mel-Stripe* locus hereafter. Details are contained below.
21

21

21 **Reference genome with Dovetail**

22

23 We generated reference genomes for a melanistic and a green morph using Dovetail
24 technology (32). For the melanistic morph we used two sequencing runs. The first run (short
25 reads + Chicago library) was done on a melanistic female from FHA caught in 2015 (id:
26 15_0190). The second run (Chicago library only) was done using another melanistic female
27 caught in 2016 in FHA (id 16_0359). For the green reference (short reads + Chicago libraries), a
28 green unstriped female from population PRC caught in 2015 was used (id 15_0802). The 2015
29 samples were flash frozen in liquid nitrogen, shipped to Sheffield and stored in a -80°C freezer.
30 It was de-gutted prior to shipping to Dovetail. The 2016 sample was caught and degutted 'fresh'
31 in California and sent directly to Dovetail.
32

32

33 The Dovetail assembly method relies on building a conventional reference assembly using
34 Meraculous with paired-end Illumina reads and then using Chicago libraries for scaffolding
35 using the HiRise pipeline (32). Chicago libraries are produced by reconstituting chromatin in
36 *vitro* with chaperones and histones, followed by crosslinking (i.e. DNA stabilization by creating
37 covalent bonds among the histones), digestion with restriction enzymes, and ligation. This
38 process results in many chimeric fragments composed from physically distant regions, but
39 ensures they come from the same stabilized large fragment. In theory, the read pairs produced
40 can have separations up to the maximum fragment size of the DNA. A model of insert
41 distribution derived from the distances among the original fragments is then used for scaffolding.
42

42

43 The assembly based on melanistic females (draft 1.3) had a 63.0x sequencing depth with a
44 total length of 953.3 Mb (73.3% of the estimated genome size by flow cytometry)(25). It
45 comprised 4068 scaffolds (N50=16.4 Mb, N90=1.1 Mb, L50=16, L90=135), a significant
improvement relative to the previous draft 0.3 (14,221 scaffolds, N50=312.5 Kb, N90=52 Kb, L50=788, L90=3869; DDBJ/ENA/GenBank accession MSSY00000000.3)(33). We clustered scaffolds in major linkage groups as described in detail below in the section on delimitation of the Mel-Stripe locus, resulting in draft 1.3c2. This Whole Genome Shotgun project was deposited at DDBJ/ENA/GenBank under the accession PGFK00000000. The version described in this paper is version PGFK01000000. The assembly based on the green female had a 42.7x sequencing depth with a total length of 932.1Mb (71.1% of the estimated size). This assembly was poorer than the 1.3, but still significantly better than the previous 0.3 (5653 scaffolds, N50=8.2 Mb, N90=503.2 Kb, L50=22, L90=222). This assembly was labeled as draft 2.1. This Whole Genome Shotgun project was deposited at DDBJ/ENA/GenBank under the accession PGTA00000000. The version described in this paper is version PGTA01000000.

12

13 Genome-wide association (GWA) mapping

14

15 We mapped color and pattern variation using previously published GBS data (33) aligned to 16 the new reference genome 1.3b2. We aligned 96.1% (789,388,267) of reads from 602 individuals 17 using BOWTIE 2.2.9 (45) with the '--very-sensitive-local' preset. We used SAMTOOLS 1.3.1 (46) 18 to sort and index the alignments. We used aligned reads with a mapping quality score of at least 19 20 to call single nucleotide polymorphisms (SNPs) with SAMTOOLS mpileup and BCFTOOLS 20 1.3.1 (46), using the original consensus caller (-c) with a P-value threshold of 0.05. From the 21 1,369,070 variants called, we excluded those with quality score of less than 20, sampling 22 coverage of less than 50%, maximum depth more than 10 times the number of total, minor-allele 23 frequency (MAF) equal or less than 0.01, and more than two alleles. The number of phenotyped 24 individuals was different for color (590) and pattern (536) and we subsequently subset variants 25 and applied filters relative to the respective number of samples. Thus, we retained 418,209 bi- 26 allelic variants for color and 416,405 variants for pattern. Both datasets were very similar, 27 showing the same mean coverage depth per SNP per individual of 5.1x (95%: 0-15; per SNP 28 average: 5.1x, 95%: 1.0-9.5; per individual average: 5.1, 95%: 2.2-7.9). We used custom Perl 29 scripts along with a custom C++ program (alleleEst 0.1b) to co-estimate allele frequencies and 30 genotypes using a Bayesian model (47). Genotype estimates were stored in BIMBAM format as 31 values ranging from 0 to 2 representing minor allele dosage. 32

32

33 Following past work (24), we used GENABEL v1.8.0 (48) to perform single locus GWA 34 mapping analyses. Briefly, we recoded genotype probabilities into genotype values accepted by 35 GENABEL using a custom Perl script as follows: [0-0.5]=homozygote for major allele, [0.5- 36 1.5]=heterozygote, [1.5-2]=homozygote for minor allele. Transformed genetic probabilities were 37 filtered using GENABEL quality control function. SNPs with MAF inferior or equal to 1%, if any, 38 were excluded from analysis. Individuals with extreme heterozygosity at a false discovery rate 39 <1% and too high an identity by state (hereafter IBS >=0.95, calculated on a subset of 2000 40 SNPs), if any, were discarded from analysis. 41

41

42 Analyses were run controlling for population structure using the GENABEL egscore function 43 (49). This function extracts principal components of a kinship matrix (here IBS indices) 44 calculated using a subset of 2000 SNPs. The principal components are then used as covariates in 45 the GWA linear models. The kinship matrix was computed excluding markers on linkage group 46 8 (to avoid over-correcting for genome-wide population structure by including causal variants),
1 and excluding markers that were not assigned to linkage groups. We display results in the form
2 of Manhattan plots. These graphics shows the association score (expressed as -log₁₀(pvalue)) of
3 every SNP tested along their physical position on the genome. Gaps between scaffolds are not
4 represented on these graphics. SNP with a significant P-value after Bonferroni correction
5 (calculated as 0.05/number of tested SNPs) are displayed in red in the Manhattan plots.

6

7 Defining the *Mel-Stripe* locus

8

9
10 We combined results from GWA mapping of color and pattern with whole genome
11 comparative alignments and recombination rate estimates from crosses to define approximate
12 boundaries for the main locus responsible for color and pattern variation in *T. cristinae* (Figures
13 S1, S2). We focused on scaffolds 702 and 128 from the melanistic genome, which contained the
14 overwhelming majority of SNPs significantly associated with color (96%) and pattern
15 (73%)(numbers refer to significance following a strict Bonferroni correction, i.e., *P* < 0.05/(no.
16 of tests)). Our approach included the following steps, which we detail below: (i) generate a
17 linkage map with the genome scaffolds, (ii) split one key scaffold (702) based on inconsistencies
18 in the linkage map, (iii) align the green and melanistic morph genomes to each other, (iv) delimit
19 the *Mel-Stripe* locus based on the total evidence from the mapping results and comparative
20 alignment. These boundaries are meant to serve as a working hypothesis for the region
21 controlling color and pattern (which can then be usefully contrasted to the genomic background),
22 and not as the precise boundaries of the functional variant(s).
23

22

23 *Linkage map*- We used the *LepMap2* software (50) and previously published data from
24 three F1 crosses to construct a linkage map for the *T. cristinae* melanistic morph genome
25 sequence scaffolds (the data, comprising 158 million ~100 base pair, bp, genotyping-by-
26 sequencing reads, are fully described in (25)(NCBI BioProject PRJNA356911). Families
27 consisted of 114 (female melanistic by male green), 48 (female green by male melanistic), and
28 24 (female green striped by male melanistic) full-sib offspring. However, note that the GWA
29 described above used a draft (1.3c2) based on only the largest family. Sequence data for the
30 parents and offspring were aligned to the melanistic morph genome using *bwa aln* and *samse*
31 (version 0.7.10-r789)(51) with a maximum of 4 miss-matches, and not more than 2 miss-matches
32 in a 20 bp seed. We then compressed, sorted and indexed the alignments using *SAMTOOLS*
33 (1.2)(46), and identified variable nucleotides using the *call* variant caller in *BCFTOOLS* (version
34 1.3)(46). We only considered alignments with a mapping quality of 10 or more and bases with a
35 base quality of 15 or more, and we applied a population prior with theta set to 0.001 when calling
36 variants and only considered a SNP if the probability of the data assuming the locus was
37 invariant was less than 0.01. We then applied a variant filter using *vcfutils varFilter* to retain only
38 those SNPs with a total read depth of 464 and that were more than 5 bp from the nearest gaps
39 (insertion-deletions).
40

40

41 We then generated the genotype input data for the mapping program, *LepMap2*. In doing so,
42 we used custom Perl scripts to select the subset of SNPs that were recombination informative for
43 each parent, and then estimated offspring genotype posterior probabilities using the genotype
44 likelihood from *BCFTOOLS* (46) (from the vcf file) with a prior given by Mendelian inheritance
45 expectation. We then only retained genotypes when the posterior probability of the most
46 probable genotype was 0.95 or greater (in other cases the genotype estimate was converted to
missing data). From this, we retained 17,478 SNPs (across all three families) for linkage map construction. As a first step with *LepMap2*, we further filtered the data for each family to retain only markers with missing data from fewer than 10 individuals, and with a P-value for segregation distortion greater than 0.005 (i.e., to remove loci with substantial deviations from Mendelian expectations). We allowed for a data error rate of 0.01. This resulted in a total of 4312 maternally informative SNPs and 5989 paternally informative SNPs.

We next used the *LepMap2 SeparateChromosomes* algorithm with a LOD minimum of 4 and with a minimum linkage group size of 50 SNPs for initial assignment of SNPs to LGs. This resulted in 6873 SNPs being assigned to 12 linkage groups (i.e., autosomes, *T. cristinae* has 13 chromosomes, see below for consideration of the sex chromosome). The *JoinSingles* algorithm was then used to assign additional SNPs to these linkage groups at the lower LOD threshold of 3, if the difference in support between their best and next best possible linkage group differed by 2 LOD units. Next we used a custom Perl script and approach to assign entire scaffolds to linkage groups based on the SNP assignments. Specifically, for a scaffold to be assigned to a linkage group (and thus all of its SNPs to be assigned to a linkage groups) required at least two SNPs (and 10% of all SNPs on a scaffold) to have been assigned to that linkage group, and for fewer than half as many SNPs to have been assigned to the next best linkage group. Based on this, we were able to assign 237 scaffolds (which accounted for 89% of all SNPs) to linkage groups. Finally, the *OrderMarkers* algorithm in *LepMap2* was used to estimate marker/SNP order on each linkage group. We took the median position in cM for all markers on a scaffold as the position for each scaffold in each cross.

As one of the filters applied with *LepMap2* was to remove markers with non-Mendelian patterns of inheritance, we expected to miss the sex (i.e., X) chromosome, and thus to find 12 of the 13 chromosomes, as we did. We thus employed a complementary approach to identify the X-linked scaffolds (in *T. cristinae* males are XO and females are XX)(52). Using SAMTOOLS DEPTH (version 1.2)(51) and custom Perl scripts, we extracted the coverage data from a previously published GBS data set that was used for genome-wide association mapping and comprised 395 female and 197 male *T. cristinae* (data from (24), but aligned to the current genome as described above; we lacked data on offspring sex in the mapping families so used this data instead). We then identified scaffolds where the ratio of read depth for males to females was less than the expected 1:1 ratio expected for autosomal markers (specifically less than 0.75). Twenty-nine scaffolds met this requirement, and also were not assigned to the 12 autosomal scaffolds described above. These included 380 recombination informative markers. Seventeen of these scaffolds were joined into a single linkage group (presumably the X chromosome) using the *SeparateChromosomes* algorithm in *LepMap2* with a LOD limit of 1.5 and a minimum size of 50 SNPs. The 17 scaffolds included 93.4% of the SNPs on the 29 scaffolds we identified as possibly being X-linked based on the coverage ratio. We used *OrderMarkers* to order these markers as described for the X-chromosome.

*Splitting and re-mapping scaffold 702*- Scaffold 702 from the melanistic morph genome showed a strong association with color and pattern in GWA analyses, but was not originally assigned to a linkage group. Upon examining this further we noted that one large chunk (SNPs up to position 14,171,514) of this scaffold was assigned to linkage group 8 (the linkage group where another scaffold, 128, showed a strong association with color and pattern and where we
1 had previously seen associations with these traits) whereas a second large chunk (SNPs after
2 position 14,757,049) was assigned to linkage group 5 (preventing placement of this scaffold).
3 The portion assigned to linkage group 8 showed an association with color and pattern, whereas
4 the other half of the scaffold did not. Based on this evidence we inferred that this scaffold was
5 over-assembled and thus we split scaffold 702 into three new scaffolds: 702.1 (positions 1-
6 14,171,514), 702.2 (starting at position 14,757,049) and 702.3 (the middle ambiguous section
7 lacking an informative SNP from 14,171,414-14,757,049). The new scaffolds 702.1 and 702.2
8 were added to their respective linkage groups and the *OrderMarkers* algorithm in *LepMap2* was
9 re-run for these linkage groups.

10
11 *Whole genome comparative alignment and defining Mel-Stripe* - We aligned the melanistic
12 and green morph genomes to each other using *Mugsy* (v1r2.3)(53). Our goal was twofold: (i) to
13 refine the orientation of scaffolds 702.1 and 128 (the two scaffolds with the greatest association
14 with color and pattern) based on overlap between these and scaffolds from the green morph
15 genome, and (ii) to identify possible structural variants associated with the GWA color and
16 pattern signal. Scaffold 702.1 (from the melanistic genome) partially aligned to green scaffold
17 1575; green scaffold 1575 also aligned to melanistic scaffold 2963 (which was 'left' of scaffold
18 702.1). Melanistic scaffold 2963 showed a negative correlation between SNP map position and
19 physical position, suggesting it was in a reverse orientation. This combined with the overlap of
20 both melanistic scaffolds 2963 and 702.1 with green scaffold 1575 allowed us to also specify
21 (flip) the orientation of 702.1. Melanistic scaffold 128 was in the correct orientation based on the
22 correlation (positive) between SNP physical and cM positions. Many small green morph
23 scaffolds with uncertain orientations span the right side of the re-orientated melanistic scaffold
24 702.1 and melanistic scaffold 128 (> 15 small scaffolds). No green scaffold was found that
25 aligned the portion of melanistic scaffold 128 from approximately 5 to 6.4 mbps. This region
26 also exhibits lower sequence coverage in green individuals, suggesting it might be a large
27 insertion-deletion polymorphism.

28
29 Given this information and the GWA mapping signal, we defined the bounds of a putative
30 Mel-Stripe color and pattern locus as comprising melanistic scaffold 702.1 starting from the edge
31 of the alignment with green scaffold 1575 (702.1 4,139,489 bp) to the edge of 702.1 (bp 1, given
32 the reverse orientation) along with the neighboring melanistic scaffold 128 from bp 1 to right
33 edge of the putative insertion-deletion polymorphism (bp 6,414,835). This specific region (that
34 is, the *Mel-Stripe* locus) contains 70% and 31% of SNPs associated with color and pattern,
35 respectively (59% of color or pattern-associated SNPs). As a comparison, this region only
36 contains about 1% of the sequenced SNPs. Thus there is a 61 and 31-fold enrichment of color
37 and pattern associated SNPs, respectively, in the *Mel-Stripe* locus. As emphasized above, our
38 main goal was to delimit a *Mel-Stripe* locus that could be contrasted to the genomic background,
39 and not to precisely identify causal functional variants affecting color and pattern. A schematic
40 summary of the delimitation of *Mel-Stripe* can be found in Figure S1.

41
42 Genomic change at the *Mel-Stripe* locus

43
44 We quantified changes at *Mel-Stripe* between time periods using three published data sets:
45 (1) genotyping-by-sequencing (GBS) data from 1102 individuals collected in a natural
46 population on *Adenostoma* (FHA) in 2011 and 2013 (n = 500 and 602, respectively)(30, 33), (2)
1 491 re-sequenced whole genomes from an eight-day (i.e., within-generation) release and
2 recapture field experiment (30), and (3) GBS data from 451 individuals in a between-year (i.e.,
3 between-generation) field transplant experiment (25). The within-generation experiment
4 involved releasing 500 *T. cristinae* in a paired-block design and recapturing the survivors (30).
5 We obtained whole genome sequence data from 491 of these individuals (33), allowing us to
6 compare allele frequency changes between release and recapture. As described previously (25),
7 the between-generation experiment involved transplanting 2000 stick insects from a single
8 variable population (OGA) onto 10 host plant bushes in a block design (five blocks each with
9 one *Adenostoma* bush and one *Ceanothus* bush per block; 200 *T. cristinae* were released on each
10 bush). 421 F1 descendants of these individuals were then captured the following year (2011). We
11 compared 30 individuals representative of the founders (collected in 2010) to the 421 F1s.
12 Phenotypic change (proportion at time period two minus proportion at time period one) for each
13 of these three data sets was as follows: FHA, stripe change = 0.06, unstriped change = -0.11,
14 melanistic change = 0.05; within-generation experiment, stripe change = 0.05, unstriped change
15 = -0.04, melanistic change = 0.01; between-generation experiment, stripe change = -0.24,
16 unstriped change = 0.32, melanistic change = -0.07).
17

18 The GBS data were aligned to the *T. cristinae* reference genome with *bwa* (version 07.10-r789)(51) using the *aln* and *samse* algorithms. We allowed 5 miss-matches, 2 miss-matches in an
19 initial 20 bp seed, trimmed bases with phred-scaled quality scores lower than 10, and only placed
20 reads with a single best match. We then used SAMTOOLS (version 1.2)(46) and the BCFTOOLS
21 call algorithm (version 1.3)(46) to identify SNPs and calculate genotype likelihoods. We used
22 the recommended mapping quality adjustment (-C 50), skipped alignments with a mapping
23 quality less than 20 and bases with a base quality less than 30, and used the multi-allelic SNP
24 caller with θ set to 0.001 and a posterior probability of 0.01 or less for the homozygous reference
25 genotype given the data to consider a SNP variable. We then filtered the initial set of SNPs to
26 retain only those with a mean coverage of >= 2X (per individual), total coverage (across all
27 individuals) less then three standard deviations above the mean across all loci, at least 10 reads
28 of the non-reference allele, a mapping quality of 30, sequence data for at least 80% of the
29 individuals, a minimum minor allele frequency of 0.01, less then 1% of reads in the reverse
30 orientation (with our GBS method all reads should be in the same orientation), and separated by
31 at least 5 bps. Filtering was done using custom Perl scripts. Following filtering, we retained
32 178,141 SNPs for the natural FHA population and 249,074 SNPs for the between-generation
33 experiment.
34

35
36 We aligned the whole genome re-sequence data from the within-generation experiment to
37 our reference genome using the *bwa* (version 07.10-r789)(51) mem algorithm with a band width
38 of 100, a 20 bp seed length and a minimum score for output of 30. We then used SAMTOOLS (46)
39 to compress, sort and index the alignments, and *Picard Tools* to mark and remove PCR
40 duplicates (version 2.1.1)(https://broadinstitute.github.io/picard/). We then used the *GATK*
41 HaplotypeCaller and GenotypeGVCFs modules (version 3.7)(54) to call variants and calculate
42 genotype likelihoods. We used a minimum base quality score of 30 for consideration in
43 calculations, a prior probability of heterozygosity of 0.001, and called variants with a minimum
44 phred-scaled confidence of 50. The following filters were then applied using custom Perl scripts:
45 minimum coverage of 1x per individual, a minimum value of the base quality rank sum test of -
46 8, a minimum value of the mapping quality rank sum test of -12.5, a minimum value of the read
1 position rank sum test of -8, a minimum ratio of variant confidence to non-reference read depth
2 of 2, a minimum mapping quality of 40, a maximum phred-scaled P-value of Fisher's exact test
3 for strand bias of 60, and a minimum minor allele frequency of 0.01. The resulting 6,175,495
4 SNPs were used for downstream analyses.
5
6
7 We obtained maximum likelihood estimates of allele frequencies for all populations /
8 experimental samples using an expectation-maximization (EM) algorithm, as described in (55).
9 For this, we used a convergence tolerance of 0.001 and allowed for a maximum of 20 EM
10 iterations.
11

12 Population genomic parameters were then calculated based on the *Mel-Stripe* locus and
13 additional reference loci based on the maximum likelihood allele frequency estimates. Additional
14 loci were defined for all genome scaffolds placed on linkage groups that contained as many
15 SNPs as *Mel-Stripe* and were defined by selecting (at random) a contiguous block of SNPs of the
16 same number as *Mel-Stripe* (FHA: 780 SNPs, 40 reference loci; between-generation experiment:
17 1180 SNPs, 39 reference loci; within-generation experiment: 47,305 SNPs, 16 reference loci).
18
19 We analyzed genomic change based on raw allele frequency changes, allele frequency
20 changes controlling for underlying genetic diversity (i.e., residual change), and using Wright's
21 Fixation Index ($F_{ST}$). Specifically, we calculated nucleotide diversity ($\pi$) within the 2011 FHA
22 sample or the founders of each experiment, allele frequency change between these samples and
23 the 2013 FHA sample (natural FHA population) or recaptured stick insects (both experiments),
24 the residuals from regressing change on diversity, and $F_{ST} = \Sigma(\pi_{total} - \pi_{subpop})/\Sigma(\pi_{total})$. In all cases,
25 $F_{ST}$ was =
26 0.0051 ($P$ = 0.024, Exact probability). For the within-generation experiment, raw change was =
27 0.0340, residual change was = 0.00212, and $F_{ST}$ was = 0.0030 ($P$ = 0.059). For the between-
28 generation experiment, raw change was = 0.0988, residual change was = 0.0595, and $F_{ST}$ was =
29 0.0540 ($P$ = 0.025; Fisher's combined probability test across data sets: $X^2 = 20.50$, d.f. = 6, $P$ =
30 0.0023).
31
32 Autoregressive-moving-average models fit to different long-term evolutionary data sets
33
34 We fit Bayesian autoregressive-moving-average (ARMA) models to 10 evolutionary time-
35 series data sets (details of each data set are given below; two are from *T. cristinae* and the others
36 from published data in other systems). This approach uses past observations as covariates in a
37 model. There are two specific types of terms in these models, autoregressive terms (AR) and
38 moving-average terms (MA). AR terms use the data values from prior years as covariates
39 whereas the MA terms use residuals from prior years as covariates. Different numbers of prior
40 years (i.e., different order models) can be considered.
41
42
43 Specifics of the models are as follows. We first considered models with order 0, 1 or 2 for
44 the auto-regressive and moving-average components of the model; a null model with a constant
45 expectation was included for comparison. As an example, ARMA (1,2) denotes a model with
46 order 1 for the autoregressive component and order 2 for the moving-average component,
meaning that information from the last year is used for the autoregressive component and that

47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146

Autoregressive-moving-average models fit to different long-term evolutionary data sets

We fit Bayesian autoregressive-moving-average (ARMA) models to 10 evolutionary time-
series data sets (details of each data set are given below; two are from *T. cristinae* and the others
from published data in other systems). This approach uses past observations as covariates in a
model. There are two specific types of terms in these models, autoregressive terms (AR) and
moving-average terms (MA). AR terms use the data values from prior years as covariates
whereas the MA terms use residuals from prior years as covariates. Different numbers of prior
years (i.e., different order models) can be considered.

Specifics of the models are as follows. We first considered models with order 0, 1 or 2 for
the auto-regressive and moving-average components of the model; a null model with a constant
expectation was included for comparison. As an example, ARMA (1,2) denotes a model with
order 1 for the autoregressive component and order 2 for the moving-average component,
meaning that information from the last year is used for the autoregressive component and that
information from the last two years is used for the moving-average component. The general form of the model is $y_t \sim \text{Normal}(\mu_t, \tau)$ and $\mu_t = c + \sum_i \theta_i y_{(t-i)} + \sum_j \phi_j \varepsilon_{(t-j)}$, where $y_{(t-i)}$ is the data value in year $i$ years in the past, $\varepsilon_{(t-1)}$ is the error term from year $j$ years, and the sums are over the order of the autoregressive and moving-average components of the model. We assumed a weakly stationary model and thus applied the re-parameterization and Beta prior scheme proposed by (56, 57). We placed a normal prior on the grand mean, $c \sim \text{Normal}(\text{mean} = 0, \text{precision} = 0.01)$, and gamma prior on the precision for the sampling distribution, $\tau \sim \text{gamma}(0.01, 0.001)$.

Each model was fit for each data set and the best model was selected based on deviance information criterion (DIC; the model with the lowest DIC was chosen). When the null model was best, the next best model was used for downstream analyses (the null model would not provide meaningful results for cross-validation or forecasting as the expectation would be the same for each year). Two estimates of DIC were obtained for each model (to verify consistency), each based on 10 Markov chain Monte Carlo (MCMC) chains each with 100,000 iterations, a 50,000 step burn-in and a thinning interval of 50. MCMC analyses were conducted using the *rjags* JAGS interface.

We then quantified the predictability of each evolutionary time series using the best ARMA model. We used two complementary approaches: leave-one-out cross-validation and forecasting. For leave-one-out cross-validation, we fit the relevant ARMA model for each data set, but with one year of the data set removed (this was done with each year in turn). The missing year's data value (evolutionary change) was then predicted from the ARMA model. We used these estimates to assess the relationship (based on a simple linear model) between the true and predicted evolutionary change.

For forecasting, we dropped the most recent *n* years of data, where *n* was (3, 4, ..., 9, 10), and fit the relevant ARMA model to predict the data values for the dropped data. We then calculated the Pearson correlation coefficient and coefficient of determination between the observed and predicted (forecast) change for the dropped years for each value of *n*. This is conceptually analogous to predicting/forecasting future (as of yet unobserved) evolutionary change. Cross-validation and forecasting results were also based on average of results from two independent MCMC model fits, each comprising 10 chains with 100,000 iterations, a 50,000 iteration burn-in and a thinning interval of 50.

The data analyzed include evolutionary time series for discrete trait frequencies, and in the case of Darwin's finches, quantitative traits (mean value). In both cases, we first obtained point estimates of the value (mean or frequency) for each generation and then converted these into evolutionary change data sets (i.e., the data point for year *i* was the value [mean or frequency] in year *i*+1 minus the value in year *i*). The nature and source of each data set are described below. Results are provided in the main text, Database S1, and Figures S3-S6.

Long-term field studies in *T. cristinae*

We compiled data on morph frequencies in *T. cristinae* using samples collected in the spring using sweep nets between 1990 and 2017. All individual were scored as ‘striped’, ‘unstriped’, or ‘melanistic’, or occasionally when it was difficult to distinguish between the first
1 two categories as 'intermediate-striped'. These classifications have been found to be highly
2 repeatable in past work (26, 35, 36, 58). Samples from 1990 to 1999 were taken and scored by
3 Cristina Sandoval, who then trained PN in 2000. PN collected and scored most samples from
4 2000 to 2017.

5
6
7 GPS coordinates of all localities were taken at and then used to estimate elevation using 'point sampling tool' on QGIS 2.16.2 (59). The elevation values were extracted from 1/3 arc-sec
8 Digital Elevation Models (DEM) at the location of the populations' coordinates. All DEMs were
9 obtained from United States Geological Survey Dataset (USGS), available at National Map
10 Viewer (https://viewer.nationalmap.gov/). Host-plant collected on (Ceanothus or Adenostoma)
11 was recorded for all individuals. We estimated the proportion of individuals in a sample that
12 were striped (% striped) using all striped and unstriped individuals (excluding melanistic
13 individuals). We estimated the proportion of individuals in a sample that were melanistic (%
14 melanistic) using all individuals. Detailed information on these localities (i.e., GPS coordinates
15 and elevations), morph frequencies, sample sizes, etc. is provided in Database S1.
16

17 We observed consistent year-to-year increases and then decreases in the frequency of
18 striped morphs at HV (see main text). We thus computed the binomial probability of the
19 observed stripe time series alternating between an increase and decrease in stripe frequency
20 every other year. Specifically, conditional on the first year, we calculated the probability that
21 every other year showed a reversal in the direction of evolution as 0.5¹⁷ = 7.6e⁻⁶ (the full time
22 series includes 18 years, the null probability that evolution reverses direction was assumed to be
23 0.5, and thus the probability of not changing direction was also 0.5).
24

25
26 Climatic data and analyses

27
28 We collated data on mean springtime statewide temperature in California using publicly
29 available records (National Centers for Environmental Information, https://www.ncdc.noaa.gov/cag/time-series/us/4/0/tavg/3/4/1990-
30 2016?base_prd=true&firstbaseyear=1901&lastbaseyear=2000). We focused on temperature
31 averages across March, April, and May as these are the three months that *T. cristinae* is by far
32 most active (most of the rest of the year is spent in egg diapause)(26, 34, 58). Nonetheless, we
33 present results from different combinations of spring months below.

34
35 We fit a hierarchical Bayesian model for the full *T. cristinae* color data set, using data from
36 all populations (i.e., not just HV) collected from 1990 to 2017. We did so to: (i) test for an
37 association between climate and the melanistic morph frequency, and (ii) determine how well
38 climate predicts color morph frequency across space and time.
39

40
41 We assumed a binomial sampling distribution for the observed number of melanistic
42 morphs for a site and year (yᵢⱼ) given the number of *T. cristinae* sampled (nᵢⱼ) and the true
43 melanistic morph frequency (pᵢⱼ). We connected this to a linear model with the logit link
44 function, such that logit(pᵢⱼ) = αᵢ + βᵢ xₜtemperature + θ xₜyear + εᵢⱼ, where αᵢ is a population (site)
45 specific intercept, βᵢ denotes the effect of climate (temperature, see details below) on melanistic
46 morph frequency for population i, θ is an overall effect of year (allowing for a general increase
47 or decrease in melanistic morph frequency), and εᵢⱼ is an error term that accounts for over-
dispersion relative to binomial sampling. We gave the ε values a normal prior with mean of 0 and
precision parameter τ ~ gamma(0.1, 0.01) (we imposed a sum-to-zero constraint on the ε values).
We then defined linear models at the next level of the hierarchy for the population specific α and
β coefficients (for the intercept and effect of temperature, respectively), such that,

$$
\alpha_i = a_1 + b_1 x_{\text{elevation}} + c_1 x_{\text{host}} + d_1 x_{\text{mountain}}
$$

$$
\beta_i = a_2 + b_2 x_{\text{elevation}} + c_2 x_{\text{host}} + d_2 x_{\text{mountain}}
$$

Here, *x*<sub>elevation</sub> is the elevation at a location, *x*<sub>host</sub> is a binary indicator variable for host plant (Adenostoma = 0, Ceanothus = 1), *x*<sub>mountain</sub> is a binary indicator variable denoting the mountain range (0 = Highway 154; 1 = Refugio), and *a*<sub>1</sub>, *a*<sub>2</sub>, *b*<sub>1</sub>, *b*<sub>2</sub>, *c*<sub>1</sub>, *c*<sub>2</sub>, *d*<sub>1</sub> and *d*<sub>2</sub> are regression coefficients (all given Normal priors with mean 0 and precision 0.0001).

We fit this model with three different temperature variables: (i) mean temperature for
March, April and May (when *T. cristinae* are most active), (ii) mean temperature for February,
March, April and May, and (iii) mean temperature for February, March and April. We used the
*rjags* interface with JAGS to obtain Markov chain Monte Carlo (MCMC) parameter estimates
for the model parameters. In each case, we ran three chains, each with a 10,000 iteration burn-in,
25,000 post burn-in iterations and a thinning interval of 10. We used four-fold cross-validation to
determine the predictive power of the models. Specifically, we split the data set into four random
subsets (only considering cases where the sample size was 25 or greater) and used three subsets
to fit the model and validated the model by predicting morph frequencies for the other subset
(MCMC options identical to those for the main models were used).

Temperature was generally associated with a higher frequency of melanistic *T. cristinae* (*a*<sub>2</sub> was positive), but less so at *Ceanothus* sites (*c*<sub>2</sub> was negative) (Table S2; estimates of the effect of temperature for each site and year are shown in the main text). Melanistic morphs were less common at higher elevations and on Refugio independent of temperature. Cross-validation results showed that the models had significant but modest predictive power. For example, with the March, April, May temperature model, the Pearson correlation between observed and predicted melanistic morph frequencies was *r* = 0.16 (95% CIs = 0.040-0.28, *P* = 0.0102, *r*<sup>2</sup> from a linear model = 0.027). The other temperature variables gave similar results: February, March, April, May temperature, *r* = 0.15 (95% CIs = 0.025-0.27, *P* = 0.0188, *r*<sup>2</sup> from a linear model = 0.022); February, March, April temperature, *r* = 0.19 (95% CIs = 0.069-0.31, *P* = 0.0024, *r*<sup>2</sup> from a linear model = 0.037).

Thermoregulatory experiments

We conducted lab thermoregulatory experiments testing the desiccation / heat tolerance of green versus melanistic *T. cristinae*. Heat stemmed from a desk lamp (K-mart model ksn: 0-02546202-9), raised 4.5 inches above two petri dishes that were stacked on top of each other and pushed to touch the base of the lamp. A third petri dish containing an individual *T. cristinae* was placed on top of the other two. The bulb used a Sylvania A19 halogen 100-watt replacement that used 72-watts. A total of four such lamp set-ups were used, allowing simultaneous assays of four *T. cristinae* (always two green and two melanistic, assigned randomly to one of the four lamps at
1 the initiation of an assay, and then randomly re-assigned to one of the four after each weighing
2 census, see below). Details of the procedure were as follows. Each individual was weighed. Each
3 lamp was then turned on for ten minutes. Placing test animals underneath the lamps then started
4 the trials. Every twenty minutes all four individuals were removed simultaneously and weighed
5 in a random order, and scored as dead or alive. They were then assigned randomly back to one of
6 the four test lamps. This procedure was repeated until 180 minutes had passed. A total of eight
7 sets of such trials were run (total n = 32).

8

9 We fit a Cox proportional hazards model to the survival data to test for an effect of morph
10 (green versus melanistic) on survival (60). For this, we used the survival package in R (61). We
11 used the exact partial likelihood method, which is advantageous relative to the more common
12 Efron method when time is measured in discrete intervals and tied times of death are thus more
13 likely. We detected a significant effect of morph on survival time (exp(B) = 3.57, 95% CIs =
14 1.34-9.51, P = 0.0111). Note that exp(B) > 1 indicates melanistic morphs died from desiccation
15 more rapidly than green morphs.

16

Estimating genotype-specific fitness using genomic data

18

19 We estimated selection coefficients/relative fitnesses for different genotypes at the *Mel-*
20 *Stripe* locus based on the within-generation release-recapture experiment and based on patterns
21 of evolutionary change between the 2013 and 2011 FHA samples. Similar to (23) we used PCA
22 and k-means clustering to assign individuals one of six *Mel-Stripe* genotypes: homozygous for
23 the stripe haplotype/allele (s/s), homozygous for the green unstriped haplotype (u/u),
24 homozygous for the melanistic haplotype (m/m), or one of the three possible heterozygotes (s/u,
25 s/m or u/m)(Fig. S2). We conducted a PCA on the individual genotype matrix for each of the two
26 data sets. This was done for all individuals and the 780 SNPs comprising the *Mel-Stripe* locus.
27 We then clustered *T. cristinae* based on the first two genetic PCs with k-means clustering; this
28 was done with the R *kmeans* function with six centers, 100 starts and a maximum of 200
29 iterations. An initial round of clustering was performed to define cluster centers. For this round
30 an equal number of green, striped and melanistic individuals were used (42 of each, which was
31 the number of green individuals). We then used those centers to cluster all individuals with a
32 second round of k-means clustering (this included individuals with no phenotypic data).
33 Assignments from k-means clustering corresponded well with groups of individuals with the
34 same color/pattern (i.e., stripe) phenotype, and were the basis for designating genotypes.
35

36

37 For the within-generation experimental data, we fit a Bayesian beta-binomial model to infer
38 fitness values. Here, we inferred the survival probability of individuals with each genotype using
39 a binomial sampling distribution for the number of recaptures given the probability of survival
40 and recapture (p<sub>genotype</sub>) and the number of individuals released with that genotype (n<sub>genotype</sub>). We
41 assigned an uninformative (Jeffery's) beta prior for each survival probability. Posterior samples
42 (N = 5000 each) were obtained from the closed form solution for the posterior using R (62), and
43 were then used to calculate the relative fitness of each genotype by dividing the survival
44 probability by the survival probability with the highest fitness (based on the point estimate; s/s).
45

46

An alternative model was required for the FHA data, which was based on change over two
generations (2011 versus 2013). During this time haplotype frequencies went from *m* = 0.316, *s*
1 = 0.602, and $u = 0.082$ to $m = 0.360$, $s = 0.570$, and $u = 0.071$. Perhaps more importantly, in both
2 years, we detected an excess of the $s/m$ heterozygotes (0.514 in 2011 and 0.502 in 2013) relative
3 to Hardy-Weinberg expectations (0.380 and 0.410, respectively). For this analysis, we assumed
4 the following relative fitness values: $s/m = 1$ (based on observed patterns of change this genotype
5 appeared to have the highest fitness), $m/m = 1 + s1$, $s/s = 1 + s2$, $u/u = 1 + s2 + s3$, $s/u = 1 + s2 +
6 s3 * s4$, and $m/u = 1 + s1 + s3 * s4$. Thus, $s1$ and $s2$ define the fitness value of the $m/m$ and $s/s$
7 homozygote in a way that allows for any form of dominance. In turn, $s3$ defines the fitness of $u/u$
8 relative to $s/s$ (i.e., after adding $s2$). The $s/u$ heterozygote is $1 + s2 + s3 * s4$, thus $s4$ is the
9 heterozygous effect. This is similar for $m/u$. We took an approximate Bayesian computation
10 (ABC) approach to estimating the selection coefficients. We first sampled selection coefficients
11 from their priors, $U(-0.5, 0.5)$ for $s1$, $s2$, and $s3$, and $U(0,1)$ for $s4$. We then simulated evolution
12 forward in time for two generations according to a Wright-Fisher model with the observed
13 starting genotype frequencies, and dynamics governed by drift and the sampled the selection
14 coefficients (assuming viability selection). We assumed a variance effective population size of
15 110.3, which was inferred from patterns of change across 178,141 SNPs (following general
16 procedures outlined in (63)). We ran 1,000,000 ABC simulations. We then used the ridge
17 regression adjustment method in the R `abc` package to obtain samples from the posterior
18 distribution from the simulation output. We retained the top 0.5% of simulations with the
19 smallest distance between the simulated and observed genotype frequencies in the 2013 sample.
20 We then converted the estimates of selection coefficients to relative fitnesses.
21
22 Field experiment testing for NFDS
23
24
25 We implemented a field transplant experiment testing for NFDS. A total of 1000 individuals
26 were transplanted, collected from March 21-24, 2017 from populations PRNC (latitude 34.53,
27 longitude -119.85), OUTA (latitude 34.53, longitude -119.84), HVC (latitude 34.49, longitude -
28 119.79), and HVA (latitude 34.49, longitude -119.79). Numbers were as follows: green-unstriped
29 morphs, PRCN 220, OUTA 50, HVC 140, HVA 90; green-striped morphs, PRCN 30, OUTA
30 100, HVA 280, HVC 90. Individuals were kept in groups of 10 and each group was randomly
31 assigned to one of two treatments: striped individuals common (40 striped and 10 unstriped
32 individuals) versus striped individuals rare (10 striped and 40 unstriped individuals). Each of
33 these groups of 50 individuals was then randomly assigned to one of 20 experimental bushes (in
34 the general area of latitude 34.51 and longitude 119.80). Each bush was cleared of existing `T`.
35 `cristinae` (the only `Timema` species occurring in this area) by sampling it each day March 21-24.
36 Past work demonstrates that this clears bushes of the overwhelming majority of `Timema` (25, 26,
37 Nonetheless, as an additional measure for ensuring accurate identification of experimental
38 animals, each transplanted individual was marked with fine tip sharpie on the underbelly. This
39 marks are not visible when `Timema` are resting on leaves. Individuals were released on March
40 26th between 9am and 3pm. Each individual was released with tweezers onto an experimental
41 plant and checked to cling well to their transplanted host. Individuals were recaptured using
42 visual surveys and sweep nets on March 31st, as in past work (25, 26, 30, 35, 36, 58), and scored
43 as striped or unstriped.
44
45
46 We fit a Bayesian beta-binomial model to assess the effect of initial stripe frequency on the
recapture stripe frequency. We assumed that the recapture stripe count for bush $i$ was $y_i \sim
1 binomial($p_i$, $n_i$), where $p_i$ is the true stripe recapture rate for a given initial release stripe

2 frequency. We then placed independent, uninformative beta priors on $p_i$ for each treatment.

3 MCMC (via *rjags*) was then used to draw samples from posterior distribution. Stripe frequencies

4 clearly increased when stripe was initially rare (recapture frequency = 0.46, 95% CIs - 0.37-0.55;

5 change in stripe frequency = 0.26, 95% CIs = 0.17-0.35; posterior probability that stripe

6 increased > 0.99). In contrast, we found no clear, consistent pattern of change when stripe was

7 initially common (change in stripe frequency = -0.006, 95% CIs = -0.093-0.063; posterior

8 probability that stripe increased > 0.43). We inferred selection coefficients for each treatment

9 (20% vs. 80% initial stripe frequency) based on the estimated posterior distribution for the true
10 stripe recapture rate. We defined relative fitnesses for striped and green stick insects as $w_{stripe} = 1$

11 and $w_{green} = 1 - s$, respectively. Here $s$ is the selection coefficient. We then estimated $w_{green}$ based

12 on the difference between release and recapture frequencies of the striped morph, such that $p_i =$

13 $(p_0 w_{stripe})/(p_0 * w_{stripe} + (1-p_0) * w_{green})$, which can be rearranged as $w_{green} = (p_0 * p_i - p_0)/((p_0-1) *$

14 $p_i)$. Here $p_0$ is the stripe release frequency (0.2 or 0.8).

15

16 Estimation of differences between hosts

17

18 We fit a hierarchical Bayesian model to quantify the overall difference in stripe frequency
19 between hosts across years. A key aspect of this model was that it allowed us to account for the
20 heterogeneity in sampling, including the fact that a subset of sites was sampled each year. We
21 used all samples from the main mountain, Highway 154. This included 21,067 data points (T.

22 *cristinae* scored as striped versus unstriped, we excluded melanistic morphs) from 274

23 collections (site by year combinations; 29 sites with a mean of 9.4 visits per site) spanning 27
24 years (1990 to 2017).

25

26 We specified generalized linear models for the stripe frequency at each location (site) for
27 each year (nearby or inter-digitated samples from different hosts were considered different sites).

28 We included effects for site and year, and modeled each of these hierarchically by placing a
29 normal prior on them with parameter values estimated from the data (except the means for the
30 year effects, which were fixed at 0 to ensure the model parameters were identifiable). We placed

31 uninformative priors on the site means, normal with mean 0 and precision $1e^{-6}$, and on the
32 precision parameters, gamma(0.01, 0.001). We used Markov chain Monte Carlo to generate

33 samples from the posterior distribution and used these samples to compute several key derived
34 parameters: the yearly mean stripe frequency for each host and the yearly mean difference in

35 stripe frequency between hosts. Inferences were based on three MCMC chains, each with a
36 10,000 iteration burn-in, 20,000 sampling iterations and a thinning interval of 5 (MCMC
37 analyses were conducted with *rjags*). Point estimates (posterior medians) for the difference

38 between hosts (stripe frequency on *Adenostoma* minus *Ceanothus*) ranged from 0.30 to 0.64

39 (mean = 0.56), and for all but one year (2011) the 95% CIs for the difference in stripe frequency
40 excluded 0 (i.e., they were significantly positive).

41

42 Estimating predictability in finches and moths

43

44 The data analyzed were obtained as follows. We obtained data on *Geospiza fortis* and

45 *Geospiza scandens* body size and beak size from (21). The data are from Daphne Major from

46 1973 to 2012. Three measurements were included: principal components (PC) 1 body size, PC1
1 beak size and PC2 beak size. We obtained data on *Panaxia dominula medionigra* allele

2 frequency from (40). We used the data from 1940-1978, as there were no gaps in sampling

3 during this time interval. We obtained data on *Biston betularia* peppered moth morph frequency

4 from (41). We used the data from Leeds, which was most complete, and restricted analysis to
5 years 1967 to 1995 because there were several years after 1995 with very low sample sizes.

6 ARMA Models were fit to the data as described for *T. cristinae* above.

7

8 We then asked whether and to what extent including rainfall data on Daphne Major (also

9 from 1973 to 2012) improved the fit of the *Geospiza* time series data sets. We focused on rainfall
10 as it is thought to be a strong determinant of seed size, which is a key source of selection on
11 these finches (1, 21). We obtained the rainfall data from (21). We fit Bayesian ARMA models of
12 order 0, 1, or 2 with respect to the AR and MA components (as described previously) that also
13 included rainfall (MCMC details were identical to those described above). We placed an
14 uninformative prior, Normal(mean = 0, precision = 1e-5), on the coefficient for rainfall. We then
15 used the best ARMA model that included rainfall (based on DIC) for predictive cross-validation
16 and forecasting as described above for the pure ARMA models (without rainfall). We then
17 compared the predictive performance of the best ARMA models with and without rainfall.
1

melanistic genome

<figure><img src="image_5.png" /></figure>

2

green genome

3

Fig. S1. Schematic illustrating the delimitation of the Mel-Stripe locus using two reference genomes. See text of supplementary materials for details.

4

5
<figure><img src="image_6.png" /><figcaption>Fig. S2. Principal Components Analysis (PCA) ordination of 1102 *T. cristinae* from FHA based on genetic data from the Mel-Stripe locus. Points (left panel) and numbers (right panel) denote individuals, and are colored based on color and pattern phenotypes (we did not have phenotypic data for some individuals). In the right panel, numbers denote cluster/group assignments from k-means clustering with k=6. Cluster assignments were used to assign genotypes when estimating selection.</figcaption></figure>

1
2
3
4
5
6
7
8
9
<figure><img src="image_7.png" /><figcaption>Fig. S3. Evolutionary time series for *Geospiza fortis* body size (a), beak PC1 (b), beak PC2 (c), *G. scandens* body size (d), beak PC1 (e), beak PC2 (f), *Panaxia dominula medionigra* frequency (g), and *Biston betularia* "peppered" frequency.</figcaption></figure>

1

2

3

4

5
<figure><img src="image_8.png" /><figcaption>Fig. S4. Change in mean trait values or morph/allele frequency for *Geospiza fortis* body size (a), beak PC1 (b), beak PC2 (c), *G. scandens* body size (d), beak PC1 (e), beak PC2 (f), *Panaxia dominula medionigra* frequency (g), and *Biston betularia* "peppered" frequency. Data points for each year denote that change observed from that year to the next year.</figcaption></figure>

1

2

3

4

5

6

7
<figure><img src="image_9.png" /><figcaption>Fig. S5. Predictive r² from ARMA forecasting models for evolutionary time series in *Geospiza fortis* body size (a), beak PC1 (b), beak PC2 (c), *G. scandens* body size (d), beak PC1 (e), beak PC2 (f), *Panaxia dominula medionigra* frequency (g), and *Biston betularia* "peppered" frequency. r² between the observed and predicted values of change are shown from models dropping (and predicting) the last three to 10 years (r² was computed from a simple linear model).</figcaption></figure>

1
2
3
4
5
6
7
8
<figure><img src="image_10.png" /><figcaption>Fig. S6. Predictive correlations from ARMA forecasting models for evolutionary time series in *Geospiza fortis* body size (a), beak PC1 (b), beak PC2 (c), *G. scandens* body size (d), beak PC1 (e), beak PC2 (f), *Panaxia dominula medionigra* frequency (g), and *Biston betularia* "peppered" frequency. Pearson correlations (solid line and points) and 95% confidence intervals (shaded polygons) between the observed and predicted values of change are shown from models dropping (and predicting) the last three to 10 years.</figcaption></figure>

1
2
3
4
5
6
7
8
<table><caption>Table S1. Summary of cross-validation and forecasting results (values for forecasting are medians from estimates based on 3 to 10 year forecasts). Bold font denotes cases where the ARMA model was preferred over a null model with a constant expectation.</caption><thead><tr><th>Data set</th><th>Best model</th><th>cross-validation intercept</th><th>cross-validation slope</th><th>cross-validation r<sup>2</sup></th><th>forecasting r</th><th>forecasting r<sup>2</sup></th></tr></thead><tbody><tr><td>Timema stripe</td><td>ARMA(1,2)</td><td>-0.005463</td><td>0.938310</td><td>0.6974</td><td>0.9282905</td><td>0.8618326</td></tr><tr><td>Timema color</td><td>ARMA(1,2)</td><td>0.03373</td><td>-1.24295</td><td>0.1019</td><td>-0.2959806</td><td>0.1388707</td></tr><tr><td>G. fortis body size</td><td>ARMA(2,2)</td><td>-0.07142</td><td>-1.89589</td><td>0.2581</td><td>0.2593157</td><td>0.2565920</td></tr><tr><td>G. fortis beak size (PC1)</td><td>ARMA(0,1)</td><td>-0.1872</td><td>-6.0475</td><td>0.2769</td><td>-0.0460291</td><td>0.1405218</td></tr><tr><td>G. fortis beak size (PC2)</td><td>ARMA(0,1)</td><td>0.002860</td><td>0.462132</td><td>0.03286</td><td>0.06829066</td><td>0.05675793</td></tr><tr><td>G. scandens body size</td><td>ARMA(1,2)</td><td>-0.005377</td><td>-0.159569</td><td>0.05488</td><td>0.6263869</td><td>0.3951535</td></tr><tr><td>G. scandens beak size (PC1)</td><td>ARMA(1,2)</td><td>0.02175</td><td>-0.22193</td><td>0.05206</td><td>0.2741161</td><td>0.1395220</td></tr><tr><td>G. scandens beak size (PC2)</td><td>ARMA(1,2)</td><td>-0.002938</td><td>0.546308</td><td>0.05978</td><td>0.08118622</td><td>0.18602100</td></tr><tr><td>P. dominula medionigra</td><td>ARMA(1,1)</td><td>-0.0005844</td><td>0.4419179</td><td>0.01698</td><td>0.000594893</td><td>0.180800958</td></tr><tr><td>B. betularia peppered</td><td>ARMA(1,0)</td><td>-0.05174</td><td>-1.80925</td><td>0.6584</td><td>0.15945025</td><td>0.02692756</td></tr></tbody></table>
<table><caption>Table S2. Posterior median and 95% credible intervals for key model parameters from the March, April, May melanistic morph model. All continuous covariates were standardized.</caption><thead><tr><td>Parameter</td><td>Median</td><td>Lower bound 95% CI</td><td>Upper bound 95% CI</td></tr></thead><tbody><tr><td>a<sub>1</sub></td><td>-2.31</td><td>-2.44</td><td>-2.20</td></tr><tr><td>a<sub>2</sub></td><td>0.187</td><td>0.063</td><td>0.309</td></tr><tr><td>b<sub>1</sub></td><td>-0.163</td><td>-0.260</td><td>-0.061</td></tr><tr><td>b<sub>2</sub></td><td>-0.0060</td><td>-0.0151</td><td>0.0274</td></tr><tr><td>c<sub>1</sub></td><td>0.164</td><td>-0.001</td><td>0.341</td></tr><tr><td>c<sub>2</sub></td><td>-0.197</td><td>-0.362</td><td>-0.249</td></tr><tr><td>d<sub>1</sub></td><td>-0.500</td><td>-0.749</td><td>-0.249</td></tr><tr><td>d<sub>2</sub></td><td>-0.050</td><td>-0.323</td><td>0.219</td></tr></tbody></table>

3
4
Table S3. Summary of model fit for the *Geospiza* data when rainfall is included in the model (based on rainfall and trait measurements from 1973-2012). We report the r² (mean across 3-10 years) for forecasting for the best ARMA model with rainfall, as well as the change in forecasting r², r (unsquared), and the lower and upper bounds on of the 95% confidence interval on r (lb and ub, respectively)(all of these values are averages across 3-10 year forecasts) obtained by including rainfall (positive values mean that rainfall improved the predictive forecast).

<table><thead><tr><th>Data_set</th><th>Model</th><th>r<sup>2</sup></th><th>Change in r<sup>2</sup></th><th>Change in r</th><th>lb</th><th>ub</th></tr></thead><tbody><tr><td><em>G. fortis</em></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>body size</td><td>ARMA(2,2)</td><td>0.434</td><td>0.178</td><td>0.238</td><td>0.192</td><td>0.102</td></tr><tr><td>beak PC1</td><td>ARMA(0,1)</td><td>0.174</td><td>0.034</td><td>0.365</td><td>0.027</td><td>0.116</td></tr><tr><td>beak PC2</td><td>ARMA(0,1)</td><td>0.080</td><td>0.023</td><td>0.207</td><td>0.078</td><td>0.047</td></tr><tr><td><em>G. scandens</em></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>body size</td><td>ARMA(2,1)</td><td>0.059</td><td>-0.337</td><td>-0.632</td><td>-0.421</td><td>-0.206</td></tr><tr><td>beak PC1</td><td>ARMA(1,2)</td><td>0.249</td><td>0.109</td><td>-0.476</td><td>0.014</td><td>-0.210</td></tr><tr><td>beak PC2</td><td>ARMA(1,2)</td><td>0.121</td><td>-0.065</td><td>-0.054</td><td>0.002</td><td>0.100</td></tr></tbody></table>

8
9
1 Database S1. Raw population data. See attached .csv sheet. Variable names are as follows:

2 location = population/locality, year = year collected, latitude = latitude, longitude =

3 longitude, elevation = elevation in meters, host = host plant collected on (A = *Adenostoma*,

4 C = *Ceanothus*), melanistic = number of melanistic individuals collected, striped = number

5 of striped individuals collected, unstriped = number of unstriped individuals collected,

6 intermediate = number of intermediately striped individuals collected, total = total number

7 of individuals collected, proportion_melanistic = proportion of the sample that was

8 melanistic, proportion Striped_no_mel = proportion of the sample that was striped

9 (excluding melanistics), mean_FebMarApr_temp = mean temperature in Fahrenheit for

10 February, March, and April, mean_MarAprMay_temp = mean temperature in Fahrenheit

11 for March, April, and May, mean_FebMarAprMay_temp = mean temperature in

12 Fahrenheit for February, March, April, and May, refugio yn = Mountain collected on (1 =

13 Refugio, 0 = Highway 154).

14

15

16

17

18

19

20

21
The three morphs of *Timema cristinae*

<figure><img src="image_11.png" /><figcaption>Unstriped</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Striped</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Melanistic</figcaption></figure>
A. GWAS mapping of color variation

<figure><img src="image_14.png" /></figure>

B. Allele frequency change through time
(FHA 2011 vs. 2013)

<figure><img src="image_15.png" /></figure>

C. Allele frequency change through time
(within-generation experiment)

<figure><img src="image_16.png" /></figure>

D. Allele frequency change through time
(between-generation experiment)

<figure><img src="image_17.png" /></figure>
A. Predicting evolution via multi-year forecasting

<figure><img src="image_18.png" /></figure>

B. Color morph frequencies through time

<figure><img src="image_19.png" /></figure>

C. Pattern morph frequencies through time

<figure><img src="image_20.png" /></figure>

D. Change in color morph frequencies

<figure><img src="image_21.png" /></figure>

E. Change in pattern morph frequencies

<figure><img src="image_22.png" /></figure>

F. Predicting change in color morph frequencies (r²)

<figure><img src="image_23.png" /></figure>

G. Predicting change in pattern morph frequencies (r²)

<figure><img src="image_24.png" /></figure>

H. Predicting change in color morph frequencies (r²)

<figure><img src="image_25.png" /></figure>

I. Predicting change in pattern morph frequencies (r²)

<figure><img src="image_26.png" /></figure>
A. Yearly temperature and melanistic frequencies

<figure><img src="image_27.png" /></figure>

B. Heat tolerance in lab experiments

<figure><img src="image_28.png" /></figure>

C. Selection in the within-generation experiment

<figure><img src="image_29.png" /></figure>

D. Selection in FHA between 2011 and 2013

<figure><img src="image_30.png" /></figure>
A. Posterior probability of selection coefficient per treatment

<figure><img src="image_31.png" /></figure>

B. Frequency change (% stripe recaptured minus % striped released) on each experimental bush

<figure><img src="image_32.png" /></figure>

C. Posterior probability of frequency change per treatment

<figure><img src="image_33.png" /></figure>

D. Striped morph frequency across the 25-year study period

<figure><img src="image_34.png" /></figure>
A. Predicting evolution in studies ($r^2$)

<figure><img src="image_35.png" /></figure>

B. Predicting evolution in studies ($r$)

<figure><img src="image_36.png" /></figure>