We get by with a little help from our friends: diversity begets diversity through shared adaptive genetic variation

22 Adaptive radiations involve astounding bursts of phenotypic, ecological, and species diversity. 23 However, the microevolutionary processes that underlie the origins of these bursts are still poorly 24 understood. We report the discovery of a cryptic ‘wide-mouth’ intermediate scale-eating 25 ecomorph in a sympatric radiation of Cyprinodon pupfishes which provides crucial information 26 about the evolutionary and ecological transition from a widespread algae-eating generalist to a 27 novel microendemic scale-eating specialist. We first show that this ecomorph occurs in sympatry 28 with generalist C. variegatus and scale-eating specialist C. desquamator across several 29 hypersaline lakes on San Salvador Island, Bahamas, but is genetically differentiated, 30 morphologically distinct when reared in a common garden, and sometimes consumes scales. We 31 then compared the timing of selective sweeps on shared and unique adaptive variants in both 32 scale-eating species to characterize the evolutionary path to scale-eating. We predicted that 33 adaptation to the intermediate C. sp. ‘wide-mouth’ scale-eating niche aided in the rapid 34 divergence of the more specialized scale-eater C. desquamator. Therefore, selection for shared 35 adaptive variants should occur first in ‘wide-mouth’. Contrary to our prediction, four of the six 36 sets of shared adaptive alleles in both scale-eating species swept significantly earlier in C. 37 desquamator. Adaptive introgression from desquamator into the ‘wide-mouth’ ancestor may 38 have resulted in parallel evolution of their dietary niche. Conversely, no adaptive alleles for 39 scale-eating were reused in a third sympatric molluscivore specialist C. brontotheriodes, despite 40 sharing 9% of hard selective sweeps. Our work provides a microevolutionary framework for 41 investigating how ‘diversity begets diversity’ during adaptive radiation. 42

Here we first investigate the position of the 'wide-mouth' on the ecological spectrum 137 from generalist to scale-eating specialist using a combination of morphological, behavioral, 138 dietary, and genomic data. We then estimated the demographic history of the 'wide-mouth' and 139 explore the timing of selection on shared and unique genetic variation involved in adaptation to 140 scale-eating to better understand this ecological transition. Additionally, we explore the extent to 141 which all trophic specialist pupfishes on this island share adaptive genetic variation, despite 142 occupying divergent ecological niches. Our results show that diversification of one species in a 143 radiation can promote further diversification by sharing adaptive genetic variation useful for 144 specialization. 145 146 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint shown). Shapes of points represent different lakes. PC1 is mainly described by lower jaw length 175 and PC2 by adductor mandibulae insertion height, buccal width, and neurocranium width. B) 176 Depictions of the three external measurements that best distinguished C. sp. 'wide-mouth' from 177 both C. desquamator and C. variegatus, measured using digital calipers. C-E) The relationship 178 between standard length (mm) of individuals and their C) lower jaw length, D) buccal cavity 179 width, and E) adductor mandibulae insertion height (AM insertion) across individuals of the 180 three species in Osprey Lake. 95% confidence bands for linear models in gray. 181

'Wide-mouth' occupies a distinct intermediate scale-eating niche 183
We found that 'wide-mouth' ingest scales, but at a significantly lower frequency than C. 184 desquamator (Wilcoxon Rank Sum test, P = 0.004; Fig. 3A). We did not detect any scales in C. 185 variegatus guts (Fig 3A). A previous gut content analysis of more individuals also found nearly 186 zero occurrence of scales in the C. variegatus diet [67]. All three sympatric species from Osprey 187 Lake collected on the same day from the same site differed in δ15N levels (ANOVA, P >  Table 1. Support for the five best-fitting demographic models for the evolution of the C. sp. 244 'wide-mouth' from the site frequency spectrum. The likelihood and AIC scores for the top five 245 best-fitting models estimated in fastsimcoal2 for the Osprey populations of C. variegatus (var), 246 C. brontotheroides (bro), C. desquamator (des), and C. sp. 'wide-mouth' (wid) are presented 247 here with a complete list of all models tested reported in Table S2. Change in likelihood (∆LnL) 248 represents the difference in likelihood from the expected simulated SFS in the demographic 249 model tested. Change in AIC (∆AIC) represents the difference in AIC scores from the model 250 with the smallest ∆AIC. All models presented here represent different divergence scenarios with 251 recent gene flow allowed, which were better supported than models with no gene flow or early 252 gene flow. Visual representations of the top five models are depicted in Figure S5. 255 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint

Shared signatures of selection across the three specialists in the radiation 270
Next, we investigated whether the evolution of new trophic specialists on San Salvador Island 271 may have been aided by reuse of the same underlying genetic variation. We characterized 272 regions of the genome under selection in the Osprey Lake populations of the four species using 273 the sliding window SFS-based approach SweeD (Pavlidis et al. 2013) and neutral simulation 274 thresholds based on demographic changes in effective population size over time (Fig S4; Table  275 S3). A higher percentage of the genome was under positive selection in all specialist populations 276 compared to the C. variegatus population (Fig 5A). This pattern might be expected given the 277 divergent selection pressures the specialists face to adapt to different trophic niches. Of all the 278 specialists, C. brontotheriodes exhibited both the most selective sweeps and the longest selective 279 sweeps ( Fig 5A). This suggests that C. brontotheriodes have undergone selection most recently 280 among all the species in Osprey Lake and is supported by the recent divergence time between the 281 C. brontotheriodes and C. variegatus populations (Fig 4D). 282 Additionally, we assessed patterns of genetic divergence across the four species in 283 Osprey Lake. Allowing for some level of gene flow between species, we calculated the number 284 of fixed and nearly-fixed SNPs (Fst ≥ 0.95) between all pairwise combinations of populations in 285 Osprey Lake. All specialists, including the two scale-eating species, were more genetically 286 diverged from each other than to C. variegatus species based on the number of SNPs fixed or 287 nearly-fixed between them (Fig 5B). This pattern of stronger genetic divergence between 288 specialists also held across fixed SNPs, the top 1% of Fst, and genome-wide average Fst (Table  289 S4). Next, we identified a set of candidate adaptive alleles for each specialist by filtering for 290 SNPs that were fixed or nearly fixed (Fst ≥ 0.95) relative to C. variegatus and found within a 291 hard selective sweep in the focal specialist. In this set of candidate adaptive alleles, C. 292 brontotheriodes had the fewest adaptive alleles while C. desquamator had the most (Fig 5B). was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made Next, we looked for shared patterns of genetic divergence and selection across the 310 specialist genomes. Despite divergent trophic specializations and morphology, we found 311 evidence of 44 shared selective sweeps across all three specialist populations that were not 312 shared with C. variegatus populations (Fig 5C). These shared regions under selection in all 313 specialists were significantly enriched for genes annotated for metabolic processes (Fig 5C), 314 suggesting shared selection for a metabolism associated with a more protein-rich diet across 315 trophic specialists (also see [70]). We also found evidence of shared selection across all pairwise 316 combinations of two specialists. The shared selective sweeps were shorter than any of the 317 selective sweeps unique to each of the specialists, suggesting that selection for these shared 318 regions was not the most recent or strongest in any of the specialists (Fig 5C). We also found no 319 shared fixed or nearly fixed (Fst ≥ 0.95) adaptive alleles across all three specialists, although 320 shared alleles did occur at lower frequencies (Fst < 0.62; Fig S6). 321 322 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint

Shared and unique adaptive alleles in 'wide-mouth' and C. desquamator 324
We discovered that the two scale-eating populations shared a set of the same adaptive alleles not 325 found in C. brontotheriodes or C. variegatus (Fig 5D). This consisted of 15 alleles in 6 shared 326 hard selective sweeps in C. desquamator and 'wide-mouth': 10 SNPs were in unannotated 327 regions, two were in the introns of the gene daam2, and three were in regulatory regions of the 328 genes usp50, atp8a1, and znf214 (one variant each). binding factor, and the ten unannotated variants were not associated with lower jaw size 340 variation in a previous genome-wide association study of the radiation (Richards et al. 2021). 341 Despite the divergent craniofacial features of the 'wide-mouth', none of the adaptive 342 alleles unique to the 'wide-mouth' appear to be in or near genes annotated for craniofacial 343 phenotypes in model organisms. In C. desquamator, three of the 13 sets of unique adaptive 344 alleles are in or near genes annotated for craniofacial phenotypes: a de novo non-synonymous 345 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. first step and the adaptive alleles involved may have set the stage for the more specialized scale-371 eater C. desquamator to evolve. To test this, we estimated the age at which the adaptive alleles 372 first started sweeping in each population using the coalescent-based approach starTMRCA [78]. 373 Across the 30 sets of adaptive alleles with selective sweep ages, we found an overall 374 difference in timing of selection between shared and unique adaptive alleles in the two scale-375 eater populations (ANOVA P-value = 0.00478). In desquamator, shared adaptive alleles were 376 under selection before any unique adaptive alleles (Tukey HSD P-value = 0.003217; Fig. 6). We 377 found a similar trend in 'wide-mouth': shared adaptive alleles were generally under selection 378 before those unique to the species (Fig. 6). However, this difference was not significant due to 379 one set of unique adaptive alleles near the gene slitrk5 which was the oldest estimated sweep 380 (ANOVA, Tukey HSD; P = 0.8367). 381 In contrast to our predictions, selection on 4 of the 6 shared sets of adaptive alleles 382 occurred significantly earlier in desquamator than 'wide-mouth'. Only a single set of these 383 adaptive alleles had an older median age estimate in 'wide-mouth' than desquamator, although 384 the 95% HPD intervals surrounding these point estimates overlapped between the two 385 populations (Fig. 6). The other adaptive allele set that had overlapping 95% HPD intervals 386 between the two populations were a set of unannotated de novo mutations. 387 388 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint

Shared adaptive variation may have resulted from introgression between desquamator and 399 'wide-mouth' scale-eaters 400
Finally, we explored several selection scenarios surrounding the shared adaptive alleles between 401 'wide-mouth' and desquamator. The most parsimonious answer for why these two populations 402 share adaptive alleles is that these alleles swept in their common ancestor and were retained in 403 both post divergence. However, the timing of selection on these alleles is much more recent in 404 both populations than the inferred divergence time (Fig 6), suggesting that these alleles were not 405 selected upon in their most recent common ancestor. Second, the timing of selection on most of 406 these shared adaptive alleles is significantly different between the populations (Fig 6), suggesting 407 that these alleles were not selected upon in a single event in the shared ancestor. This left us to 408 explore two other scenarios that could explain the shared adaptive alleles: 1) these alleles were 409 segregating in both populations after they initially diverged from each other and independently 410 swept to fixation or 2) these adaptive alleles introgressed from one population to another after 411 being selected for in one population. 412 To distinguish between these two scenarios, we assessed the regions surrounding these 413 shared adaptive alleles for signatures of introgression. Introgression between recently diverged 414 sister species is often more challenging to detect in comparison to introgression with a more 415 distantly related outgroup because of the high genetic similarity between sister species. 416 However, introgression between sister species can be detected from population level 417 polymorphism data by looking for regions of the genome with lower absolute genetic divergence 418 than expected given the genome-wide average [79-82]. We thus calculated Dxy in sliding 419 windows across the genomes of the Osprey Lake populations and compared the distribution of 420 Dxy values in the regions surrounding the shared adaptive alleles to the distribution of Dxy values 421 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint across the genome. The average Dxy of the shared adaptive alleles was two times lower than 422 average genome-wide values between desquamator and 'wide-mouth' (Fig 7A). It was also 423 slightly lower (1.2X) than the average Dxy for all selective sweep regions in the 'wide-mouth' 424 genomes, suggesting it is not due to selection on these regions in general (Table S5). This lower 425 Dxy was predominantly due to four out of the six regions containing shared adaptive alleles 426 having ~2-5 times lower Dxy values than the genome-wide average. The other two regions had 427 comparable (znf214) and higher (unannotated region on scaffold 43) Dxy values than the genome-428 wide average ( Fig 7B). The lower genetic distance between C. desquamator and the 'wide-429 mouth' in these adaptive regions suggests they may have introgressed between the two groups 430 rather than independent selection on ancestral polymorphisms since they diverged approximately 431 11,000 years ago. 432 433 434 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint and C. sp. 'wide-mouth', did we find shared nearly-fixed or fixed adaptive alleles. These species was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint We found that all of the shared adaptive alleles between the scale-eating populations 493 occurred at low frequency in generalist populations on other Caribbean islands in a previous 494 study ( [68]; Table S7). For example, three copies of the shared usp50 adaptive allele were found 495 outside of San Salvador Island (Table S7). This indicates that initial bouts of selection occurred 496 on available standing genetic variation and that staggered timing of hard selective sweeps in each 497 trophic specialist most likely reflects mutation-order processes in which selection on a beneficial 498 allele was contingent on prior fixation of other adaptive alleles in each specialists' genetic 499 background. However, several adaptive alleles originated from introgression or de novo 500 mutations found only on San Salvador Island (Table S7) Although the two scale-eating species shared a set of adaptive alleles and a recent common 507 ancestor, 5 of 6 shared adaptive regions swept at significantly different times between the two 508 species. This difference in timing may result from several different scenarios: 1) independent 509 adaptive walks to the same scale-eating niche, 2) independent adaptive walks to different scale-510 eating niches or 3) the adaptive walk of one population depends on the adaptive walk of the other 511 population. We explore each of these scenarios in turn. 512 First, the difference in timing of selection on the same shared adaptive alleles could 513 indicate independent adaptive walks to the same scale-eating peak that occurred at different 514 times and/or by slightly different routes. C. desquamator and 'wide-mouth' populations have 515 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint predominantly abutting ranges, with only a small amount of geographic overlap on San Salvador 516 Island in four lakes, Osprey, Oyster, Little Lake, and Mermaid's Pond (Fig 1). If this current 517 distribution is representative of their historical ranges and the two lineages began diverging 518 about 11,000 years ago, it is possible that the adaptive walks took place in different lakes and 519 were largely independent of one another. The differences in timing on the shared adaptive alleles occurred on the shared adaptive variants between the two scale-eater populations (Fig 5). alleles and selective sweeps, the majority of selective sweeps are unique to each species (Fig 5 &  534   6), including neurogenesis, brain, and nervous system development (slitrk5, sema4g, and 535 smarce1), whereas unique adaptive alleles in desquamator include craniofacial development 536 annotations (olfm1, gnaq, twist1) as well. This difference in gene annotations can also be seen at 537 the broader level of regions of the genome under selection (Fig S7). Therefore, the differences in 538 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint timing at shared adaptive alleles might reflect differences in relative strength of selection due to 539 different selective regimes experienced by the two scale-eating populations. While starTMRCA 540 is fairly robust in its timing estimates to different strengths of the selection [78], we cannot rule 541 out that differences in relative timing of selection might partially reflect differences in strength of 542 selection occurring in the two scale-eating populations. 543 Third, another scenario that could have led to the older timing of selection across shared 544 adaptive alleles in desquamator than 'wide mouth' is one in which the adaptive walk of one 545 population depended on adaptive alleles from another. The significantly older timing of 546 selection in desquamator on shared adaptive alleles and the low absolute genetic divergence 547 between these two scale-eating population in these regions (Fig 6) suggests that these alleles may 548 have introgressed between the two populations. Given that these shared adaptive alleles appear 549 as standing genetic variation across Caribbean populations, albeit at very low frequency in our 550 sampling, it is likely that the shared alleles were present in the ancestor and were segregating in was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint microbiomes in the radiation [73]. 586 We also see pairwise shared selective sweeps between all combinations of two 587 specialists. In contrast to the sweeps shared between all specialists, we see stronger genetic 588 divergence among species in these regions and a proportion of the selective sweeps shared 589 between two specialists do contain fixed or nearly-fixed alleles. However, the specialists that was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint suggesting that selection may have occurred most recently in the C. brontotheriodes. This 608 matches with the most recent divergence event being between generalist and C. brontotheriodes 609 in Osprey Lake. An intriguing implication of this result is the potentially long wait time for 610 speciation of C. brontotheriodes species compared C. desquamator despite the shallower fitness 611 valley isolating brontotheroides from generalists. However, we cannot yet rule out that the recent 612 divergence time estimated in Osprey Lake reflects a more recent arrival of C. brontotheriodes to 613 that lake after initially diverging long ago in another lake not incorporated in our demographic 614 model. 615

Methods 617
Sampling 618 The C. sp. 'wide-mouth' individuals were collected on San Salvador Island in the Bahamas 619 using hand and seine nets between 2017 and 2018. 'Wide-mouth' individuals were collected 620 from 5 lakes in which taxa occurs in sympatry with just C. variegatus (Great Lake, and Stout's 621 Lake) or with generalist and C. desquamator (Oyster Pond, Osprey Lake, and Mermaid's Pond,; 622 Multiple individuals from 3 lake populations (Osprey Lake, Great Lake, and Oyster Pond) in 634 which we had measurable specimens of the species (n=84) were measured for 9 traits using dial 635 calipers. Traits were selected for specific connections to foraging performance [103,104] and that 636 differed across the three groups in a previous study when externally measured [10]. Measurements on the full dataset were not taken until repeatability was above 85% per trait (r 2 in 650 a linear regression of both sets of measurements). To remove the effects of body size on trait 651 size, we used the residuals from a linear regression of the log-transformed trait on log 652 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint transformed standard length. We also compared these morphological measurements between lab-653 reared F2 and wild caught specimens to rule out strong plasticity in traits of interest. 654  Duplicate reads were identified using MarkDuplicates and BAM indices were created using 675 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint pooling all entries with less than 20 SNPs). For parameter estimates, we used a wide search 722 range with log-uniform distributions with the range of some priors informed by previous 723 estimates of effective population size from MSMC and the age of the last glacial maximum 724 (~20kya) prior to which the lakes would have been dry on SSI. These ranges were not upper 725 bounded by these specified priors and so the simulations were free to explore parameter space 726 that exceeded the priors.

Patterns of directional selection 739
Various demographic histories can shift the distribution of low-and high-frequency derived 740 variants to falsely resemble signatures of hard selective sweeps. In a previous study, we used 741 MSMC analyses to infer histories of population size changes in C. variegatus, C. desquamator 742 and C. brontotheriodes (Richards et al. 2021). In this study, we repeated the same analysis for 743 the 'wide-mouth'. In order to account for demography of the 'wide-mouth' population in 744 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint downstream analyses, we used the MSMC (v. 1.0.1;[112] to infer historical effective population 745 size (N e ) changes in the 'wide-mouth'. We ran MSMC on unphased GATK-called genotypes 746 separately for a single individual 'wide-mouth' from Osprey Lake with 16x mean coverage 747 across its genome (Fig S4). As recommended in the MSMC documentation, we masked sites 748 with less than half or more than double the mean coverage for that individual or with a genotype 749 quality below 20. We also excluded sites with <10 reads as recommended by Nadachowska- were at least 100-kb in length (99 scaffolds; 85.6% of the genome). We also calculated CLRs 760 across 100,000 scaffolds consisting of neutrally evolving sequences simulated with ms-move 761 [116], controlling for the impact of the inferred population size decreases over time for each 762 population from MSMC runs mentioned above (Fig. S4; Table S3). The CLR ratios for the 763 simulated datasets were then used to assess outlier CLR ratios from the empirical dataset. 764 Regions with CLR ratios above the 95 th percentile value of CLR from the neutral simulated 765 dataset were considered candidate hard selective sweep regions (Table S3). We compared 766 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint selective sweeps across C. variegatus, C. desquamator, C. brontotheriodes, and 'wide-mouth' to 767 look for shared and unique selective sweeps among the four groups. adaptive allele and all runs were checked for convergence between and within runs. Most runs 809 rapidly converged within the initial 6000 steps, but 5 runs did not converge after an additional 810 4000 steps and were discarded from further analysis. 811 812 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. relationship between genotype or phenotype and fitness ( 1 ). The landscape is a concept , 869 metaphor , and an empirical measurement that exerts substantial influence over all 870 evolutionary dynamics ( 2 -6 ). Fitness landscapes were originally depicted as high-871 dimensional networks spanning genotypic space in which each genotype is associated 872 with fitness ( 1 ). Simpson ( 7 ) later described the phenotypic evolution of populations 873 through time on a rugged landscape , in which isolated clusters of fitness peaks represent ' 874 adaptive zones ' to which populations evolve from adjacent regions of low fitness ( 8 ).

875
Lande and Arnold formalized the analysis of selection and estimation of phenotypic 876 fitness landscapes ( 9 -11 ), leading to empirical studies of fitness landscapes in 877 numerous empirical systems ( 12 -18 ). Fitness surfaces also provide a central component 878 of speciation models and theory ( 19 -21 ). Populations may speciate through a process of 879 divergent selection ( 22 ), either on static landscapes or negative frequency-dependent 880 disruptive selection ( 23 -28 ). A fundamental concept in fitness landscape theory is that 881 not all genotypic trajectories are accessible ( 5 , 29 , 30 ); only those paths that 882 monotonically increase in fitness at each mutational step are evolutionarily accessible , an 883 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; https://doi.org/10.1101/2021.07.01.450755 doi: bioRxiv preprint Figure S1. Relationships between three divergent traits and standard length across ponds. The 1198 relationship between standard length (mm) of individuals and their D) lower jaw length, E) 1199 buccal cavity width, and F) adductor mandibulae muscle insertion height across individuals the 1200 three species C. desquamator (teal), C. variegatus (gold), and C. sp. 'wide-mouth' (red-orange) 1201 across different ponds (Osprey, Oyster, and Great Lake). Colored lines represent linear model of 1202 these relationships for each species with their 95% confidence bands in gray. 1203 1204 1205 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. populations for every SNP found in regions under shared selection across two or more specialist 1256 species on San Salvador Island: C. brontotheroides (bro); C. desquamator (des), C. sp. 'wide-1257 mouth' (wid). Regions shared across all three specialists (bro+des+wid) do not contain highly 1258 divergent SNPs between any of the species in the radiation. Highly divergent SNPs (Fst > 0.75) 1259 are only found in regions shared between one or two specialists. 1260 1261 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. categories that were significantly enriched for relevant terms corresponding to craniofacial 1267 development, the major axes of morphological divergence in this radiation, are highlighted in 1268 red. All terms included were significant at an FDR < 0.05 and full list of terms in Data S1-3. . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted July 28, 2021. ; Table S2. Support for the 28 demographic models for the evolution of the 'wide-mouths' from