Connecting micro and macroevolution using genetic incompatibilities and natural selection on 1 additive genetic variance 2

Evolutionary biologists have long sought to identify the links between micro and macroevolution to better understand how biodiversity is created. Despite the pursuit, it remains a challenge to understand how allele frequency changes correlate with the evolution of morphological diversity, and the build-up of reproductive isolation. To connect micro and macroevolution, we tested the adaptive importance of alleles underlying genetic incompatibilities, and the consequences for predicting evolutionary trajectories. Using a quantitative genetics crossing design, we produced an F4 Advanced Recombinant Form (ARF) between four contrasting ecotypes, which we phenotyped in the glasshouse (N=770) and transplanted into the four natural habitats (N=14,265 seeds), alongside the parental ecotypes. F2 hybrid breakdown was associated with the loss of extreme phenotypes and environment-specific genetic variation in field performance. We found evidence of genetic trade-offs among environments, but only in axes describing smaller amounts of genetic variance for fitness. Habitats that showed stronger patterns of adaptive divergence for native versus foreign ecotypes, also showed lower genetic variance in fitness of the ARF. Integrating data from the field and glasshouse predicted patterns of selection on morphological traits in a similar direction to the parental ecotypes. Overall, our results provide strong empirical evidence linking ecotype specific alleles with fitness trade-offs, phenotypic divergence and the rise in genetic incompatibilities among recently derived ecotypes. Our data connects microevolutionary change with macroevolution through adaptive radiation, where increases in environment specific alleles create changes in the distribution of genetic variance, ameliorating genetic constraints to adaptation as adaptive divergence proceeds.


Introduction 40
Despite the fundamental importance of divergent natural selection for adaptive evolution, research lacks a 41 clear empirical relationship between micro and macroevolution, largely due to the difficulty of estimating 42 natural selection in the wild (Pujol et al. 2018), and connecting it to the accumulation of reproductive 43 isolation (Baack et al. 2015). There are some examples where divergent natural selection acting on genes 44 of major effect (e.g., Eda gene in sticklebacks) describe patterns of adaptive divergence (Barrett et al. 45 2008), but few studies have connected divergence in genetic variation underlying quantitative traits with 46 observed evolutionary trajectories and the accumulation of reproductive isolation. Consequently, a more 47 detailed understanding of natural selection is required to understand if alleles underlying adaptation are 48 also those contributing to reproductive isolation, and how this creates adaptive radiation. 49 Alleles can confer an adaptive advantage in one environment, but with neutral or deleterious effects in 50 other environments, leading to fitness tradeoffs (Anderson et al. 2011;Anderson et al. 2013). If such 51 environment specific alleles evolve with reduced or restricted gene flow, they will be novel in relation to 52 genotypes from alternative environments and fail when tested in alternative genetic backgrounds, creating 53 reproductive isolation (Coyne and Orr 2004). Under this scenario, environment specific alleles will lead 54 to Bateson-Dobzhansky-Muller genetic incompatibilities when incompatible with alternative genetic 55 backgrounds (Bateson 1909;Dobzhansky 1937;Muller 1942), meaning alleles underlying adaptive traits 56 in one ecotype may lead to hybrid breakdown when they are introgressed into an alternative ecotype 57 (Gavrilets 2003;Kondrashov 2003;Navarro and Barton 2003). If alleles underlying adaptation to 58 contrasting environments concomitantly create fitness trade-offs and reproductive isolation, we can use 59 changes in environment specific allele frequencies to connect micro and macroevolution. 60 Natural selection is unlikely to affect single traits in isolation, favoring beneficial combinations of traits 61 and the evolution of multivariate phenotypes (Lande 1979;Cheverud 1982). Adaptation will be 62 constrained when traits share genetic variance, and genetic architecture rather than natural selection 63 determines evolutionary trajectories (Lande and Arnold 1983;Arnold 1992; Arnold et al. 2008). In this 64 way, adaptive alleles will likely increase in frequency, but only if selection on genetically correlated traits 65 We created the ARF ensuring each ecotype contributed equally and ensuring that at each generation (see 120 Figure 1C), all full-sibling families (hereafter, 'families') contributed equally to the next generation. First, 121 we grew plants for the base population from seeds sampled from the natural populations and performed 122 crosses among the ecotypes (n = 41-60 individuals/ecotype) to create all combinations of F1 hybrids (n = 123 12 crossing combinations; n = 20-25 families/cross type). We then mated among all combinations of 124 crosses in the F1 generation such that all F2 families (n = 24 crossing combinations; n = 17-22 families 125 /cross type) possessed a grandparent from each of the original parental ecotypes (e.g., F1 Dune,Headland × 126 F1 Tableland,Woodland ). Given strong reductions in intrinsic fitness was observed in a previous Dune x 127 Headland F2 hybrid ), we maximized the number of F1 crosses to produce 458 F2 128 families in total. We grew one individual from each family. Reductions in fitness were observed as F2 129 hybrid sterility (42% of F2 individuals were successfully mated compared to >90% in F1 hybrids) and 130 reduced fertility (49% reduction in seed set compared to F1 hybrids) . Consequently, 131 we divided the F2 individuals that produced flowers into three replicate crossing lines to maintain 132 replicates of the construction of the ARF. We then randomly mated among all F2 individuals within each 133 line (n = 4-12 families/F2 cross type; total F2 families crossed N = 202) to produce the F3 generation (N 134 = 259 families), ensuring that each family contributed equally. We then produced the F4 generation by 135 first growing one individual from each F3 family and randomly designating each individual as a sire or 136 dam. We then mated 115 sires to 114 dams in a full-sibling, half-sibling crossing design to produce 198 137 families for the F4 generation. The numbers of families and individuals used to create each generation of 138 the ARF are listed in Table S1. 139 In the following analyses we examine results from two experiments using the ARF. In experiment 1, we 140 grew the ARF in the glasshouse to estimate genetic variance underlying morphological traits. In 141 experiment 2, we transplanted seeds of the ARF into the four habitats to compare the fitness of the ARF 142 with the parental ecotypes. We then used the field fitness of the ARF (experiment 2) to quantify the 143 genetic covariance in performance among transplant habitats and identify genotype-by-environment 144 interactions that would indicate genetic trade-offs among habitats. Finally, to quantify differences in 145 natural selection among habitats, we combined the data from both experiments and estimated the genetic covariance between morphological traits in the glasshouse (experiment 1) and field fitness in each habitat 147 (experiment 2). 148

Experiment 1: Glasshouse phenotypes 149
To estimate genetic variance underlying morphological traits we grew four individuals from each full 150 sibling family of the ARF (n = 198 full-sibling families, total N = 770 individuals) in 30 cell growth trays 151 containing the same potting media described above. Alongside the ARF we grew four individuals from 152 ~2 5 full sibling families for each of the parental ecotypes (N = 366 individuals). Plants were grown in a 153 25°C controlled temperature room with a 12h:12h day:night photoperiod. After eight weeks of growth we 154 measured plant height and sampled one fully mature leaf for each plant. We used the software 'Lamina' 155 to analyse the scanned leaf and quantify six variables relating to leaf size and leaf shape (Bylesjo et al. 156 2008). Using the outputs of Lamina, we quantified leaf morphology using leaf area, leaf area 2 / leaf 157 perimeter 2 as a measure of leaf complexity, leaf circularity, number of indents standardized by leaf 158 perimeter, leaf indent width and leaf indent depth. 159

Experiment 2: Field transplant 160
Seeds from the F4 generation of the ARF were transplanted into each of the habitats. At each transplant 161 habitat, we planted 18 seeds from each full-sibling family (n = 198) divided equally amongst six 162 experimental blocks (habitat n ≈ 3,500 seeds, total N = 14,265 seeds). Alongside the ARF we 163 replicate natural germination conditions we suspended shadecloth (50%) 15cm above each experimental 170 block and watered them daily for three weeks. During the initial three-week period we measured 171 emergence and mortality daily. Following the initial three weeks we measured survival and development 172 at weeks 4, 5, 7 and 9, and then monthly until 20 months at which time there were fewer than 20% of 173 germinated plants remained, and we ceased the experiment. The measures of fitness we recorded were: 174 whether each seedling emerged, whether each seedling reached 10 leaves (as a measure of seedling 175 establishment) and produced a bud (reached maturity). All measures of fitness were collected as binary 176 data. 177

Implementation of Bayesian models 178
In the subsequent analyses we implemented Bayesian models to 1) compare field performance 179 (experiment 2) of the ARF with the parental ecotypes, 2) identify whether genotype-by-environment 180 interactions create genetic trade-offs among transplant habitats (experiment 2), 3) estimate genetic 181 (co)variance of morphological traits for the ARF (experiment 1), and 4) estimate the genetic covariance 182 between morphological traits (experiment 1) and field performance (experiment 2), to identify differences 183 in natural selection on morphological traits, among habitats. For the analyses estimating genetic variance, comparing estimates of genetic variance with zero provides 194 an uninformative test of significance because estimates are restricted to be greater than zero (positive-195 definite). To create an informative significance test, we re-implemented each model with randomized 196 data, created by shuffling the parental information. For each model implemented on the observed data, we 197 re-implemented the same model on 1,000 randomizations of the data, and extracted the posterior mean for 198 each randomization. We then compared the distribution of means from models conducted on the 199 randomizations, to the mean of the observed posterior distribution. If the mean of the observed 200 distribution occurred outside the 95% Highest Posterior Density (HPD) interval for the random 201 distribution, we took this as evidence that we captured biologically important information for the 202 comparison of interest. As we were only interested in estimating the posterior mean of models 203 implemented on each randomization of the data, we could reduce computing time by reducing the total 204 number of sampling iterations. To do so, we maintained the same burn-in period and sampling interval to 205 ensure an identical mixing of MCMC chains, reducing only the total number of sampling iterations to the 206 number required to obtain a stable estimate of the mean. We calculated the number of sampling iterations 207 required using the models implemented on the observed data, which was different for each of the analyses 208 outlined below (Table S2). 209

Comparing ARF and ecotype morphology 210
To compare differences in multivariate phenotype between the ARF and parental ecotypes, we 211 implemented a multivariate analysis of variance (MANOVA) on the seven morphological traits measured 212 in experiment 1. We first standardized all seven morphological traits to a mean of zero, and standard 213 deviation of one before including them as a multivariate response variable. To test whether the ARF was 214 phenotypically different to each ecotype we conducted a separate MANOVA for each pairwise 215 comparison between the parental ecotypes, and the ARF. We used a bonferroni corrected α-value of 216 0.0125 (α = 0.05 / n, where n represents the number of tests). To visualize differences among all ecotypes 217 and the ARF we estimate D, the variance-covariance matrix representing multivariate phenotypic 218 divergence. To do so, we first conducted another MANOVA that included all ecotypes (but not the ARF). 219 From this, we extracted the sums of squares and cross-product matrices for the ecotypes (SSCP H ) and 220 error terms (SSCP E ) to calculate their mean-square matrices by dividing by the appropriate degrees of multivariate mean phenotype, among the parental ecotypes, after removing the residual phenotypic variation. To visualize the phenotypic space occupied by the ARF relative to the parental ecotypes, we 226 decomposed D into orthogonal axes (eigenvectors) and calculated the phenotype scores for the first two 227 eigenvectors for all ecotypes, and the ARF. 228

Comparing ARF and ecotype field performance 229
We estimated fitness at early life history stages for the ARF and parental ecotypes transplanted into all 230 four habitats. To do so, we created a dummy variable that represented the ARF and native versus foreign 231 ecotypes in each habitat. We then used MCMCglmm to implement the model, 232 (1) 233 where transplant habitat ( ! ), ARF/ecotype ( ! ) and their interaction ( ! × ) were included as fixed 234 effects. Blocks within transplant habitat ( ! ! ) and replicate genetic lines within the ARF ( ! ! ) were 235 included as random effects, and !(!"#$) represented the model error. We implemented equation 1 with 236 emergence, seedling establishment and maturity as a multivariate response variable ( !"#$% ). As such, for 237 all ecotypes and the ARF, equation 1 calculated the probability of reaching maturity, conditional on the 238 previous life history stages. 239

Quantifying divergent natural selection 240
We used the ARF to investigate differences in natural selection among contrasting natural habitats. To do 241 so, we conducted two further analyses. First, to identify whether natural selection created genetic trade-242 offs among habitats, indicated by a negative genetic covariance among habitats, we analysed field 243 performance using a genotype-by-environment covariance framework described below. Next, we 244 examined whether genetic selection on traits occurred in the direction of the native ecotypes. To do so, 245 we combined the morphology data from experiment 1, with field performance data in experiment 2 and 246 used the Robertson-Price Identity to estimate the genetic covariance between morphological traits and 247 field performance for each transplant habitat. We predicted that if natural selection on morphology 248 occurred in the direction of the original ecotypes, differences in selection gradients (among habitats) 249 would align with divergence in phenotype mean of the parental ecotypes. 250

Genetic trade-offs among contrasting habitats 251
We investigated genotype-by-environment (G×E) interactions in the ARF using a character state 252 approach, where different environments represent different traits (Robinson and Beckerman 2013). To do 253 so, we used the field performance of the ARF and implemented 254 (2) 255 where replicate genetic line of the ARF ( ! ) and transplant habitat ( ! ) were included as fixed effects. 256 We included sire ( ! !" ), dam ( !(!") ) and block within habitat ( ! ! ) as random effects, with 257 !(!"#$%) representing the residual error variance. For each term in the random component, we estimated 258 random intercepts for each habitat and the covariance among habitats. As such, for the sire and dam 259 components we estimated a 4×4 covariance matrix representing variance in each habitat, and covariance 260 among habitats. Information for estimating covariance among habitats is taken from individuals of the 261 same full-sibling families transplanted in each habitat. Consequently, we implemented equation 2 with a 262 heterogeneous residual covariance matrix. This allowed for different variances in each habitat, but fixed 263 residual covariances at zero because individuals (seeds) could not be planted in two habitats 264 simultaneously. We used three separate implementations of equation 2 for emergence, seedling 265 establishment and maturity included as binary univariate response variables (y ijklmn ). 266 From equation 2, S l(im) represents one quarter of the additive genetic variance in each habitat, and one 267 quarter of the additive genetic covariance between habitats. We multiplied the sire variance component 268 (S l(im) ) by four and used the posterior mean as our observed estimate of additive genetic (co)variance for 269 field performance among the four habitats. The diagonal of the resulting G-matrices represents additive 270 genetic variance within a transplant habitat, with the off diagonal representing the genetic covariance 271 between habitats. 272

Divergent natural selection on morphological traits 273
By linking genetic variance underlying morphological traits in the laboratory, with genetic variance 274 underlying field performance, we sought to quantify differences in natural selection among the transplant 275 habitats using the Robertson-Price Identity. A requirement for natural selection is genetic variance in both 276 morphological traits and field performance. The analysis of genetic variance underlying field 277 performance (described in the previous section) identified significant genetic variance for the ability to 278 reach maturity, in all four transplant habitats (see results). To identify the morphological traits with 279 genetic variation we used the morphology data from experiment 1 and implemented 280 where replicate genetic line of the ARF ( ! ) was included as a fixed effect and sire ( ! !" ), dam ( !(!") ) 282 and their interaction ( ! !" × !(!") ) were random effects, with !(!"#) as the residual variance. We 283 implemented equation 3 with plant height and six leaf morphology traits as a multivariate response 284 variable (y ijkl ). To prevent traits on different scales affecting the analysis, we centered all traits to a mean 285 of zero and standardized to a standard deviation of one prior to analysis. We then calculated the additive 286 genetic (co)variance matrix as four times the sire variance component. As traits were standardized prior to 287 analysis, genetic variances represent heritabilities. We found only four traits with heritabilities greater 288 than 0.1 (plant height, leaf area, leaf perimeter 2 / area 2 and leaf indent width; see Table S3), which we 289 then combined with field performance to study natural selection in the subsequent analyses, described 290

below. 291
We estimated the genetic covariance between the morphological traits and field performance by 292 implementing the Robertson-Price Identity 293 where the response to selection (R) is analogous to the selection differential ( s g), calculated as the genetic 295 covariance between a trait ( z ) and fitness ( w ). Equation 4 then generalizes to multivariate form by 296 including more phenotypic traits and estimating a genetic variance-covariance matrix (G), with fitness as 297 the final trait. In this framework, s g generalizes to the vector of selection gradients (s g ) representing the 298 multivariate response to selection. Estimating the response to selection in this way includes both direct 299 and indirect selection. To isolate the effect of direct selection on phenotypic responses, we can calculate 300 the genetic selection gradient by combining G and s g with 301 where ! now represents a vector of genetic selection gradients (Lande and Arnold 1983;Rausher 1992), 303 after removing the effect of genetic correlations among traits. 304 To estimate the predicted response to selection (s g ) in the ARF we estimated the (co)variance between the 305 four morphology traits and field performance by implementing 306 307 where replicate genetic line ( ! ) was the only fixed effect. Sire ( ! !" ), dam ( !(!") ) and their interaction 308 ( ! !" × !(!") ) were included as random effects along with block within habitat ( ! ). The multivariate 309 response variable ( !"#$% ) included four phenotypic traits as well as ability to reach maturity in each 310 habitat. Fitness and morphology was measured on separate individuals (field versus glasshouse 311 experiments), and so similar to equation 2, we estimated a heterogeneous residual covariance matrix. vector of covariances between morphological traits and field performance (rows one to four of the fifth 320 column). 321 To identify whether we captured biologically meaningful differences in selection among habitats, we 322 conducted two analyses. First, s g and ! being vectors, we calculated the dot product (representing vector 323 length) of the observed and random matrices. If the observed length was greater than the length calculated 324 from the random distribution, we took this as evidence we detected biologically meaningful estimates of 325 selection for a given habitat (Stinchcombe et al. 2014;Walsh and Lynch 2018). Second, to quantify the 326 differences in s g among the four transplant habitats we estimated variance is selection gradients using 327 suggested we captured greater among-habitat differences in selection than expected by random sampling. 334

Comparing ARF and ecotype morphology 336
Ecotypes showed strong differences in leaf morphology ( Figure 1A  The first eigenvector of D (d max ) described 84% of phenotypic divergence mostly created by phenotypic 344 differences between the Tableland and Headland ecotypes ( Figure 1D). The second eigenvector (d 2 ) 345 described 14% of variation in multivariate phenotype, describing differences between the Woodland and 346 Tableland ecotypes ( Figure 1D). The ARF occupied an area in phenotypic space close to the Dune 347 ecotype, and intermediate between the Headland, Tableland and Woodland ecotypes. However, the ARF 348 mean was not similar to that of the overall mean of all ecotypes, but exhibited high phenotypic variance 349 that appeared to be missing some of the extreme phenotypes, especially from the Tableland and 350 Woodland ecotypes ( Figure 1D).

Genetic trade-offs among contrasting habitats 368
If habitat-specific natural selection creates genetic trade-offs between contrasting habitats, we expected 369 the ARF to show genetic variance for field fitness, and negative genetic covariance between contrasting 370 habitats. However, given the ARF lacked the extreme phenotypes of the Headland, Tableland and 371 Woodland ecotypes (Figure 1), and exhibited relatively high field performance (Figure 2), we might 372 expect low genetic variance associated with either zero genetic covariance or a positive genetic 373 covariance between habitats. We found that across the four habitats, additive genetic variance increased 374 as life history stages progressed (Figure 3). Observed estimates of genetic variance in field performance 375 were within the random distribution at emergence, but were greater than the random distribution for 376 maturity (Figure 3). In the headland and tableland habitats we detected lower genetic variance than 377 expected by chance for seedling establishment.

383
Decomposing the genetic covariance matrix described orthogonal axes of genetic variation underlying 384 field fitness. Decomposing the matrix for each life history stage, we found the first three eigenvectors for 385 maturity described more genetic variance than expected under random sampling (Figure 4). Interpreting 386 the loadings of each eigenvector reveal how each habitat contributes to describing the genetic variance in 387 fitness quantified by that eigenvector. Habitats with loadings of the same sign describe shared genetic 388 variance for fitness, whereas loadings of different signs describe differences in genetic variance and 389 provide evidence of fitness trade-offs. We found all habitats contributed equally to describing genetic 390 variance underlying the first eigenvector, suggesting it described heterosis or shared genetic variation 391 needed to function in stressful environments (Table 1). However, eigenvectors two and three provided 392 evidence of genetic tradeoffs, describing genetic variance in fitness that differed between the woodland 393 and dune ecotypes (e 2 ), and between the tableland, and the dune and woodland transplant habitats (e 3 ; 394 Table 1). Eigenvector 4 did not describe biologically meaningful genetic variance (Figure 4), but 395 described differences in genetic variance between the headland, and dune and tableland habitats. The

406
Overall, our results showed strong patterns of adaptive divergence (Figure 2), and although there appears 407 to be a common genetic basis to fitness in all environments (e 1 ; Table 1) we also detected genetic trade-408 offs for fitness among certain habitats (Table 1). Despite strong adaptive divergence in Figure 2 Eigenvector Genetic variance adaptive genetic variation that was lost in the ARF, reducing genetic variance for field performance in 413 certain environments and producing weaker genetic trade-offs than expected. To test this, for each habitat 414 we compared the strength of adaptive divergence (Figure 2; native ecotype performance -foreign ecotype 415 performance) against the level of genetic variance exhibited by the ARF. As predicted, we found a strong 416 negative association for seedling establishment and a weaker negative association for maturity ( Figure 5), 417 suggesting alleles associated with strong adaptive divergence were also responsible for genetic 418 incompatibilities.

Natural selection on morphological traits 427
To quantify selection in each habitat we calculated s g as the genetic covariance between morphological 428 traits measured in the glasshouse, and field performance measured in each of the four transplant habitats. 429 We then isolated direct selection by calculating ! , the genetic selection gradient for each habitat. 430 Comparing the length of observed and random s g and ! suggested we captured biologically meaningful 431 selection within each habitat ( Figure S5A). To quantify differences in selection among habitats we 432 estimated B and Z as the among-habitat (co)variance in selection vectors. Comparison of observed and 433 random eigenvalues showed that both selection vectors exhibited greater differences among habitats than 434  (Figure S5B), suggesting differences in our observed selection vectors 435 described biologically meaningful differences in natural selection among transplant habitats. 436 If differences in natural selection among the four habitats occurred in the direction of adaptive evolution, 437 we would expect differences in s g , but not ! , to align with divergence in mean phenotype of the parental 438 ecotypes. Eigenanalyis of B and Z quantifies the axes that describe differences among the original 439 selection vectors, with the first axis for each matrix representing 83% (HPD 56-98%) and 81% (HPD 55-440 98%) of the total variance, respectively. We tested whether the first axis from each selection vector 441 aligned with d max , the axis describing the greatest difference in multivariate phenotype mean. To do so, 442 we calculated the angle between the first eigenvector of B and Z, and d max . We found the alignment 443 between Z and D, but not B, was closer than expected with random sampling ( Figure 6A). To 444 complement this analysis, we conducted a more extensive analysis using a covariance tensor approach, 445 which is provided in supplementary material. Results obtained from both analyses matched closely, 446 suggesting the response to selection, but not the direction of selection, aligned with divergence in parental 447 ecotype morphology.  The ARF exhibited a multivariate phenotype intermediate to the four parental ecotypes, and lacking in 455 some of the more extreme phenotypes. Field performance of the ARF was similar, or higher, than the 456 native ecotypes, suggesting heterosis has persisted. Despite this, we found trade-offs in field performance 457 among habitats, but only in axes describing smaller amounts of genetic variance underlying field fitness. The ARF also showed reduced genetic variance for fitness when transplanted into habitats associated with 459 stronger adaptive divergence, providing evidence that divergent natural selection shaped genetic variation 460 in traits associated with fitness. Among-habitat differences in natural selection on morphological traits 461 aligned with the direction of morphological divergence of the original ecotypes, providing evidence for 462 divergent natural selection towards the original ecotypes, despite only one generation of selection. this system, it is likely the evolution of coadapted gene complexes, rather than speciation genes (Nosil 485 and Schluter 2011), were likely responsible for the rise of intrinsic reproductive isolation during the early 486 stages of speciation (Corbett-Detig et al. 2013). We can then view the evolution of these ecotypes from a 487 perspective that combines Fisher's geometric model (Fisher 1930) and Wright's focus on epistasis 488 (Wright 1931), where selection acts upon additive genetic variation by increasing allele frequencies at 489 independent loci, but limited recombination due to small population size or strong selection creates 490 coadapted gene complexes, making them functionally epistatic. The strength of coadaptation within a 491 population will then determine how genetic incompatibilities arise among populations and lead to 492 speciation. 493 The strength of divergence (and hence, reproductive isolation) among coadapted gene complexes will be 494 population and environment specific, and depend on the interaction between mutation, migration, drift these contrasting environments are repeatedly selected in the same environment, they may form 499 coadapted gene complexes within each population, with drift or local adaptation causing differences 500 among localities (Goodnight 2000). Whether locally adapted coadapted gene complexes between 501 locations of the same species will give rise to reproductive isolation remains unexplored, but could 502 provide important insights into the relationship between adaptation and divergence across a 503 heterogeneous landscape. 504 Genetic variance for life history and fitness traits is often low (e.g., McFarlane et al. 2014), and often 505 decreases with ontogeny (e.g., Aguirre et al. 2014). In contrast, our results suggested increases in genetic 506 variance with development, which has been found previously in laboratory conditions (Styga et al. 2018). 507 Changes in genetic variance underlying fitness have profound implications for understanding adaptation 508 and responses to environmental change (Sgrò and Hoffmann 2004). If genetic correlations among traits 509 under selection change during ontogeny, the effects of selection will not be linear over time and will 510 depend on changes in the combination of genetic variation and selection pressures over time. As different 511 trait combinations will be available to selection at different developmental points, patterns of adaptation 512 will be determined by the combination of traits visible when selection is strong (Bourret et al. 2017;513 Styga et al. 2018). Consequently, it will be important to consider the relationship between changes in 514 genetic correlations and changes in natural selection, as development proceeds. 515 The alignment between differences in the response to selection (s g ), but not the genetic selection gradients 516 (β g ), and phenotypic divergence (D), suggests that constraints to adaptation would exist if the ARF was 517 left to evolve in the natural environments. This is because after one generation of selection, the mean 518 phenotype was expected to follow divergence towards the parental ecotypes, but selection in the absence 519 of genetic correlations among traits was in a direction different to phenotypic divergence. We must be 520 which aligned with the direction of morphological evolution. This suggested that genetic constraints have 525 limited capacity to constrain adaptation during adaptive radiation, or genetic variance can evolve to 526 reduce genetic constraints as evolution proceeds. Given we observed divergence in ecotype G, but also 527 genetic constraints in the ARF after ecotype-specific adaptive alleles were lost, we believe the loss of 528 ecotype-specific adaptive alleles has re-created the constraints present during the very early stages of 529 adaptive divergence. Thus, adaptive radiation can occur when environment-specific alleles increase in 530 frequency, causing changes in the distribution of genetic variance and ameliorating genetic constraints as 531 adaptive divergence proceeds. Therefore, during adaptive radiation, the early stages of adaptive 532 , but longer-term evolution is determined by the long-term correlated response to selection (Zeng 534 1988). 535 In this work, we identified patterns of phenotypic divergence and adaptive divergence among recently 536 derived ecotypes as a result of the accumulation of environment-specific alleles in response to natural 537 selection. We show that these alleles likely created fitness trade-offs between habitats, lead to outlier 538 phenotypes absent in alternative ecotypes and were likely the cause of genetic incompatibilities. Through 539 these experiments we identify the connection between microevolutionary genetic changes and 540 macroevolutionary diversification in the context of an adaptive radiation. 541