Parallel genetic changes underlie integrated craniofacial traits in an adaptive radiation of trophic specialist pupfishes

Many factors such as divergence time, shared standing genetic variation, frequency of introgression, and mutation rates can influence the likelihood of whether populations adapt to similar environments via parallel or non-parallel genetic changes. However, the frequency of parallel vs non-parallel genetic changes resulting in phenotypic convergence is still rarely estimated. In this study, we used a QTL mapping approach to investigate the genetic basis of highly divergent craniofacial traits between scale- and snail-eating trophic specialist species across similar hypersaline lake environments in an adaptive radiation of pupfishes endemic to San Salvador Island, Bahamas. We raised F2 intercrosses of scale- and snail-eaters from two different lake populations of sympatric specialists, estimated linkage maps, and scanned for significant QTL for 30 skeletal craniofacial and body traits. We then compared the location of QTL in each lake for each trait to determine whether convergence in this system is due primarily to parallel or non-parallel genetic changes. We found strong support for parallel genetic changes in both lakes for five traits in which we detected a significant QTL in at least one lake. However, many of these shared QTL affected different, but highly correlated craniofacial traits in each lake, suggesting that pleiotropy and trait integration should not be neglected when estimating rates of parallel evolution.We also found that shared QTL regions between lakes were six times more likely to contain at least one highly differenciated variant under selection in wild pupfish populations than the rest of the genome. We further observed a 23-52% increase in adaptive introgression within shared QTL, suggesting that introgression may be important for parallel evolution. Overall, our results suggest that the same genomic regions contribute to the convergent phenotypes of snail- and scale-eating across lakes, and that those genomic regions contain more introgressed adaptive alleles than expected by chance. We also highlight the need for more expansive searches for shared QTL when testing for parallel evolution.


Introduction
Convergent evolution describes the independent evolution of similar phenotypes in response to similar selective pressure and can provide support for ecological adaptation (Losos 2009;Schluter 2000). Historically, it was assumed that convergence was brought about via different genetic mechanisms (hereafter referred to as non-parallel genetic changes); however, recent work now provides empirical examples of the same genetic mechanisms (hereafter referred to as parallel genetic changes) leading to convergent phenotypes (Besnard et al., 2009;Bricelj et al., 2005;Chan et al., 2010;Christin, Salamin, Savolainen, Duvall, & Besnard, 2007;Cresko et al., 2004;Feldman, Brodie, Brodie, & Pfrender, 2009Jost et al., 2008). Despite this substantial attention, the frequency and likelihood of convergence via parallel or non-parallel changes is still relatively unknown (Rosenblum, Parent, & Brandt, 2014;Stern, 2013;Stern & Orgogozo, 2008).
Many factors influence whether phenotypic convergence is produced by parallel or nonparallel genetic mechanisms. First, recently diverged species exhibit increased probabilities of genetic parallelism when adapting to similar environments. Recently diverged taxa may inhabit environments that are similar more frequently or they may have similar genetic architectures, similar genetic variance-covariance matrices, or similar genetic backgrounds that produce similar epistatic interactions (Conte et al. 2012;Rosenblum et al. 2014). Second, any mechanism that allows the use of the same adaptive genetic mechanism should increase the likelihood of convergence via parallelism, including the availability of shared standing genetic variation and introgression (Rosenblum et al., 2014). For example, threespine sticklebacks colonized freshwater thousands of times and converged on similar phenotypes largely due to selection on an ancient shared pool of marine standing genetic variation (Jones et al. 2012;Feulner et al. 2013;Nelson and Cresko 2018;Haenel et al. 2019; but see: Chan et al. 2010;Stuart et al. 2017).
Finally, de novo mutations are more likely to result in convergent phenotypes via non-parallel genetic pathways, particularly when there is a large mutational target size or polygenic adaptive phenotype (Wittkopp et al. 2003;Kowalko et al. 2013;Bolnick et al. 2018, but see: Colosimo et al. 2004;Chan et al. 2010;Xie et al. 2019).
Quantitative trait locus (QTL) mapping is often used to infer whether parallel or nonparallel genetic changes underlie convergent phenotypes. However, many QTL studies only investigate a limited number of traits that are controlled by large effect loci, which may bias the literature towards supporting genetic parallelism (Conte, Arnegard, Peichel, & Schluter, 2012).
This bias may be exacerbated by the fact that in many QTL studies the genomic regions associated with a convergent trait are large, contain many genes, and their effects on phenotypic variance are overestimated in under-powered studies (Beavis, 1998). These methodological and experimental limitations reduce confidence in the specific genomic regions associated with a convergent phenotype and, by extension, reduce confidence in whether convergence was brought about by parallel or non-parallel genetic changes. One possible solution is to compare the genomic regions associated with many different phenotypes across populations (Erickson et al., 2016). In this scenario, shared genomic regions across populations provide strong support for genetic parallelism, except in the likely rare circumstance of independent de novo mutations within the same region (O'Brown et al. 2015;Xie et al. 2019;Chan et al. 2010).
The San Salvador Island, Bahamas (SSI) pupfish radiation is an excellent system for investigating the genetic underpinnings of parallel ecomorph phenotypes because novel trophic specialists occur in sympatry across multiple hypersaline lake populations on the island. The radiation includes three pupfish species: a generalist pupfish (Cyprinodon variegatus), a scaleeating pupfish (C. desquamator), and a snail-eating (durophagous) pupfish (C. brontotheroides; Martin and Wainwright 2013). The snail-and scale-eating pupfish are endemic to SSI, and occur in sympatry with one another and the generalist pupfish.
Among lakes, specialists have converged on multivariate phenotypes that are adaptive for their given ecological niche. For example, scale-eaters across all lakes exhibit increased oral jaw size (Martin & Wainwright, 2013;Hernandez et al. 2018) and reduced lower jaw angles during scale-eating strikes which may play a critical role in scale-biting performance during high-speed strikes on their prey (St. John, Holzman, & Martin, 2020). Similarly, the snail-eating pupfish exhibits a novel nasal protrusion which may improve snail-eating performance or result from sexual selection (Martin & Wainwright, 2013;St. John, Dixon, & Martin, 2020). Furthermore, the nasal protrusion of the snail-eating speices varies substantially among lake populations (Hernandez et al., 2018;Martin & Feinstein, 2014). Despite the importance of these species characteristics, we still do not understand how their genetic architecture varies across populations.
There is some evidence to suggest that parallel genetic changes underlie specialist phenotypes on SSI. First, the SSI radiation is very young, diverging only about 10 kya (Hagey & Mylroie, 1995). Second, previous genomic analyses suggest that many of the alleles associated with trophic specialization arrived on SSI from Caribbean-wide standing genetic variation within generalist pupfish populations, but there also appear to be a few de novo adaptive mutations associated with scale-eating (Richards et al., 2021). Scale-eaters form a monophyletic group, suggesting a shared genetic component to the scale-eating phenotype across lakes (Richards & Martin, 2017). In contrast, snail-eaters and generalists often genetically cluster together by lake instead of by species-suggesting that non-parallel genetic changes could underlie convergent snail-eater phenotypes (Martin and Feinstein 2014;Richards and Martin 2017). Furthermore, previous studies have documented strong genetic divergence between scale-eaters from Crescent Pond and all other populations of scale-eater (Richards & Martin, 2017;Richards et al., 2021).
In this study we mapped the genetic basis of 30 craniofacial and body traits associated with snail-and scale-eating using lab-reared F2 intercrosses from Crescent Pond and Little Lake.
We called variants, estimated linkage maps, and performed QTL analyses independently for each F2 population. We found that only one trait-cranial height-mapped to the same genomic region in both Crescent Pond and Little Lake, but 4 of the 5 remaining significant QTL detected in one lake mapped to the same genomic region as a highly correlated craniofacial trait in the second lake. Ultimately, we conclude that parallel evolution through reuse of introgressed adaptive alleles is acting to produce snail-and scale-eating phenotypes across lake populations on SSI.
For individuals from Crescent Pond, we recorded their sex using their breeding coloration. Male pupfish develop a blue iridescent coloration along their anteriodorsal surface and a black marginal band along their caudal fin (Echelle & Echelle, 2020).
Once F2 hybrids reached sexual maturity, we performed mating assays using a subset of the hybrid females from Crescent Pond to estimate mate preferences for snail-or scale-eating mates (N=74). Prior to the mating assays, female fish were isolated for at least twelve hours and conditioned on frozen bloodworms with a 12:12 light:dark cycle. Mating assays occurred in three 1.1 m diameter kiddie pools (5-10ppts salinity). Pools were covered with gravel substrate and divided in half. In each half, we placed three clear plastic 7.5-L Kritter Keepers in a row containing three conspecific males housed individually to avoid aggression. Size-matched scaleeater males were placed on one side of each arena and snail-eating males on the other. Once the males were placed in their clear boxes, a female F2 hybrid from Crescent Pond was placed into the center of one of the three the kiddie pools, chosen at random. We considered females acclimated to the pool once they had visited both rows of males, after which time we started the seven-minute trial period. During each trial we recorded the amount of time a female spent within one body-length of each species. Each female was tested consecutively in all three pools, and we used the mean of her association time (scale-eater association time / total association time during each 7-minute trial) across all three pools for QTL analysis. Size-matched males were periodically rotated through kiddie pools during testing.

Morphological Traits
To measure skeletal phenotypes in our F2 intercrosses, we cleared and double-stained each specimen with alizarin red and alcian blue. Before clearing and staining, each fish was skinned and fixed in 95% ethanol. We then fixed specimens in 10% buffered formalin for at least one week and stained batches of individually labeled specimens following Dingerkus and Uhler's (1977) protocol. We suspended cleared and stained specimens in glycerin, and photographed their left lateral side using a Canon EOS 60D digital SLR camera with a 60 mm macro lens. For each individual, we took two types of photographs: first, we took a whole-body photograph to calculate fin and body measurements and second, a lateral skull image to calculate craniofacial measurements (Error! Reference source not found.). We used the DLTdv8 software to digitize 11 landmarks on each whole body image and 19 landmarks on each lateral skull image following the morphometric methods laid out in Martin et al. (2017). For individuals from Crescent Pond, we also weighed the adductor mandibulae muscle mass. Each image included a standardized grid background which we used to calibrate and transform our measurements from pixels into millimeters. In total, we measured 354 individuals from Crescent Pond and 287 individuals from Little Lake. We used R to convert the 30 landmarks into linear distances. To reduce measurement error due to the lateral positioning of the specimens, we took the mean distances from the two clearest skull and whole-body photographs for each individual when possible. If an individual did not have two clear photographs for each position, we measured the single clearest photograph. Finally, we size-corrected each trait by using the residuals from a linear model including the log-transformed measurement of each trait as the response variable and log-transformed standard length as the predictor variable. We investigated whether sizecorrected traits varied between the two populations using a PCA and a MANOVA test, but found no appreciable difference between them (Error! Reference source not found., num df = 28, approximate F-value= 0.34, P = 1)

Genotyping
We genotyped individuals using three different methods: First, we used whole genome resequencing for the wild-caught F0 parental generation of our Crescent Pond and Little Lake intercrosses. We used DNeasy Blood and Tissue Kits (Qiagen, Inc.) to extract DNA from the muscle tissue of each fish and quantified it on a Qubit 3.0 fluoromether (Thermofisher Scientific, Inc.). Genomic libraries were then prepared at the Vincent J. Coates Genomic Sequencing Center (QB3) using the automated Apollo 324 system (WaterGen Biosystems, Inc.). Samples were fragmented using Covaris sonication and barcoded with Illumina indices. A quality check was also performed on all samples using a Fragment Analyzer (Advanced Analytical Technologies, Inc.). We used 150 paired-end sequencing on an Illumina Hiseq4000 for these four parental samples along with an additional 38 samples that were included in a previous study (Richards & Martin, 2017).
Second, in addition to the 190 previously sequenced individuals from Crescent Pond used for a QTL mapping study , we included an additional 164 F2 individuals from Crescent Pond sequenced using double-digest restriction site associated sequencing (ddRADseq) following similar library prep and sequencing methods described in Martin et al. (2015Martin et al. ( , 2016Martin et al. ( , 2017. Briefly, we prepared four indexed libraries each containing 96 barcoded individuals. We sequenced these using 100 single-end high-output mode on two lanes of Illumina Hiseq4000 at the Vincent J. Coates Genomic Sequencing Center (QB3).
Finally, we sequenced all F2 individuals from Little Lake and a subset of previously sequenced, but low-coverage Crescent Pond F2's (N=84), using Nextera-tagmented reductivelyamplified DNA (NextRad) sequencing (Russello, Waterhouse, Etter, & Johnson, 2015). We followed the above methods for DNA extraction and sent samples to SNPsaurus (SNPsaurus, LLC) for quality checking, NextRad library preparation, and 150 single-end sequencing on two lanes of Illumina Hiseq4000 at the University of Oregon sequencing core.

Calling Variants
We used the following methods to call variants separately for: 1) the Crescent Pond intercross (2 parents and 354 F2 hybrids), and 2) the Little Lake intercross (2 parents and 285 F2 hybrids): First, we inspected raw read quality using FastQC (Babraham Bioinformatics Institute,v0.11.7) and trimmed reads to their appropriate length (100bp for samples sequenced with ddRAD, and 150bp for samples sequenced with NextRAD) using TrimGalore! (v0.6.4). For samples that were sequenced using both ddRAD and NextRad methods, we concatenated trimmed raw reads into a single file. We next used bwa-mem to map reads from all individuals in an intercross, both parents and offspring, to the Cyprinodon brontotheroides reference genome (v 1.0; total sequence length = 1,162,855,435 bp; number of scaffolds = 15,698, scaffold N50 = 32 Mbp; (Richards et al. 2021)). We identified duplicate reads using MarkDuplicates and created BAM indices using BuildBamIndex in the Picard package (http://picard.sourceforge.net(v.2.0.1)).
Following the best practices guide from the Genome Analysis Toolkit (v 3.5; (Depristo et al., 2011)), we called and refined our single nucleotide polymorphism (SNP) variant data set using the program HaplotypeCaller. Pupfish lack high-quality known alleles because they are a nonmodel organism; we therefore used the recommended hard filter criteria (QD < 2.0; FS <; 60; MQRankSum < -12.5; ReadPosRankSum < -8; (Depristo et al., 2011;Marsden et al., 2014)) to filter our SNP variant dataset. Ultimately, we detected 13.7 million variants in our Crescent Pond dataset and 14.4 million variants in our Little Lake dataset.
We used the program STACKS to further filter our dataset and convert our vcf files into phenotype and genotype comma-separated values files that could be imported into the Rqtl program. Specifically, we used the populations program to filter out variants that were not present in both the parental and F2 populations, and to filter out variants found in 10% or less of the population. From this filtering step we retained 36,318 variants with 46.5 mean mappable progeny per site in Crescent Pond and 87,579 variants with 85.984 mean mappable progeny per site in Little Lake.
We continued to filter our datasets using the Rqtl (v1.46-2), and ASMap (v1.0-4) packages (K. W. Broman, Wu, Sen, & Churchill, 2003;J. Taylor & Butler, 2017). We started filtering by removing individuals that did not contain any filtered variants and any duplicate individuals. This reduced our Crescent Pond data set to 227 individuals, and our Little Lake data set to 281 individuals. Next, we filtered markers that had >0.98 or <0.1 heterozygosity (Crescent Pond: markers =15,247, Little Lake: markers=14,661). This step also filtered out 13 individuals from Crescent Pond which only contained markers with >0.98 or <0.1 heterozygosity. Before constructing our genetic maps, we set aside markers that appeared to suffer from segregation distortion. We used the pullCross() function from the ASmap package to set aside markers in both data sets that were missing in >75% of individuals, departed from Mendelian ratios (1:2:1), or any co-located markers for the initial construction of the linkage maps. This filtering retained more than twice the number of markers for Crescent Pond than Little Lake. We therefore used a stricter filtering threshold for missing data (i.e., removing markers with >72% missing data) for our Crescent Pond dataset to construct linkage maps of comparable sizes for downstream comparative analyses. At the end of this filtering process the Crescent Pond dataset contained 214 individuals and 657 markers and the Little Lake dataset contained 281 individuals with 490 markers.

Linkage Map Construction
We used the mstmap.cross() function to form initial linkage groups and order markers, using the kosambi method for calculating genetic distances and a clustering threshold of P=1 x10 -14 for Little Lake and P=1x10 -20 for Crescent Pond. After forming these initial linkage groups, we used the pushCross() function from the ASmap package to integrate previously set aside markers back into our map. We pushed markers back based on a segregation ratio of 3:4:3, and we pushed back any markers that had previously been designated as co-located. This increased our map sizes to 817 markers for Crescent Pond and 580 markers for Little Lake. With these additional markers, we re-estimated our linkage map using the est.rf() and formLinkageGroups() functions from the Rqtl package. We used a max recombination fraction of 0.35 and a minimum LOD threshold of 5 to estimate linkage groups for both data sets. We used the droponemarker() command from Rqtl with an error probability of 0.01 to identify and drop problematic markers from the genetic maps, including dropping linkage groups with 3 or fewer markers. Finally, we visually inspected our linkage groups using plotRF() from the Rqtl package, and merged linkage groups which had been incorrectly split up using the mergeCross() function from the ASmap package. Ultimately our final genetic maps included: 1) Crescent Pond: 214 individuals, 743 markers, 24 linkage groups and 2) Little Lake: 281 individuals, 540 markers, and 24 linkage groups (Figure 2).

QTL Analyses
We mapped QTL for 29 skeletal traits for both populations, and additional morphological (adductor mandibulae muscle mass) and behavioral traits (mate preference) for Crescent Pond.
We used the Rqtl2 package (v0.22-11) to calculate genotype probabilities with a multipoint hidden Markov model using an error probability of 0.0001 and a Kosambi map function. We also calculated overall kinship and kinship using the leave-one-chromosome-out method (LOCO), which we used to reduce bias while scanning for QTLs. We used the scan1() function to perform three separate genome scans using a single-qtl model by: 1) Haley-Knott regression, 2) a linear mixed model using the overall kinship matrix, and 3) a linear mixed model using the LOCO kinship matrix. For our Crescent Pond data set we also included sex as an additive covariate. We assessed the significance of all three models using two significance thresholds P < 0.1 and 0.05 based on 1000 permutations each, using the scan1perm() function. As noted above the scan1() function can use several different methods to determine if a region is significantly associated with a given phenotype (Broman et al., 2019;Haley & Knott, 1992;Yang, Zaitlen, Goddard, Visscher, & Price, 2014;Yu et al., 2006), however, it is clear from previous theoretical work that many of these methods may suffer from type II error depending on the size of an organism's genome, the density of markers in a linkage map, or the complexity of the phenotypic traits being measured (E. S. Lander & Botstein, 1989;Risch, 1990). We therefore relaxed the P-value cut off from 0.05 to 0.1 to capture potentially important genomic regions. This relaxation is further supported by the LOD scores associated with regions significant at the P<0.1 level because they all exceed the traditional threshold of 3 (Nyholt, 2000), the more conservative threshold of ~3.3 (E. Lander & Kruglyak, 1995;Nyholt, 2000), the suggestive threshold of 1.86 (E. Lander & Kruglyak, 1995), and are in line with estimates of significant LOD thresholds in previous studies (Erickson et al., 2016). All three of these methods detected similar QTLs and moving forward we only used the Haley-Knott regression method.
For each trait, we calculated the location of the maximum LOD score, and used the fit1() function to re-fit a single-QTL model at this location. We used the newly calculated LOD score to estimate the proportion of variance explained by the QTL and to calculate a P-value associated with each significant QTL ‫ݔ(‬ 2 -test). We also used the location of the maximum LOD score to calculate 95% Bayes credible intervals using the bayes_int() function from the Rqtl2 package. We note that the maximum LOD score associated with every trait across both ponds exceeded the suggestive threshold of 1.86 (Lander and Kruglyak 1995). Finally, we used the find.markerpos() function from Rqtl to determine where markers in each linkage map fell within the reference genome. With this information we were able to determine the scaffolds/positions from the reference genome that fell within the 95% credible intervals for each putative QTL.

Identifying adaptive alleles within QTL regions
For each scaffold that fell within a QTL's credible interval we calculated the minimum and maximum position for that scaffold (that was identified in the putative QTL region) and searched the C. brontotheroides reference genome for annotated genes within the region. We then compared this list to a previously published list of genes that 1) contained or were adjacent to (within 20Kbp) fixed or nearly fixed (Fst > 0.95) SNPs between specialist species on SSI, and 2) showed significant evidence of a hard selective sweep in both the site frequency spectrum-based method SweeD (Pavlidis, Živković, Stamatakis, & Alachiotis, 2013) and the linkagedisequilibrium-based method OmegaPlus (Alachiotis et al. 2012). We hereafter refer to these loci as adaptive alleles. We also noted whether adaptive alleles within QTL regions were classified as de novo, introgressed, or as standing genetic variation on SSI (Richards et al., 2021). We used a bootstrapping resampling method to determine whether the observed proportions of adaptive alleles originating from de novo, introgression, or standing genetic variation found within QTL 95% credible interval were different than the proportions expected when drawn from the genome at random. We used the boot (v. 1.3-25) package (Buckland, Davison, & Hinkley, 1998;Canty & Ripley, 2021) to resample our entire adaptive allele dataset (with replacement) 10,000 times.
We then used the boot.ci() command from the boot package to calculated the 95% credible intervals around expected proportions of de novo, introgressed, and standing adaptive alleles. We performed these calculations separately for scale-eater and snail-eater adaptive alleles.

Linkage Map Construction
We identified 24 linkage groups from 743 markers for Crescent Pond, and 24 linkage groups from 540 markers for Little Lake (Figure 2). Previous karyotypes of Cyprinodon species estimated 24 diploid chromosomes, matching the linkage groups in this study (Liu & Echelle, 2013;Stevenson, 1981). The total map length for Crescent Pond was 7335 cM and the total map length for Little Lake was 5330; the largest linkage groups for each map were 740 cM and 380 cM, respectively, and inter-marker map distance did not exceed 20cM in either map. To compare our maps and to determine if the same genomic regions were being reused across lakes, we identified where each marker was located in our reference genome. Overall, we found 324 markers in both maps that were within 10 Kbp of one another, indicating that 60% of the Little Lake map was also present in the Crescent Pond map and 44% of the Crescent Pond map was present in the Little Lake map (Figure 3).

Craniofacial QTL
We detected three significant QTL in Crescent Pond and five QTL in Little Lake ( Table 1,   Table 2). In Crescent Pond, we identified QTL associated with the depth of the dentigerous arm of the premaxilla, cranial height, and adductor mandibulae muscle mass. Cranial height in Crescent Pond mapped to linkage group (LG) 10. Dentigerous arm depth and adductor mandibulae muscle mass both mapped to LG 13, which also contained the max LOD scores for two additional jaw traits (jaw opening in-lever and maxillary head height; Table 2). The 95% credible intervals for all these traits overlapped, suggesting that linkage group 13 may contain a single pleiotropic locus or many loci that affect all four traits.
In Little Lake, we detected significant QTL associated with jaw closing in-lever (i.e. height of the coronoid process on the articular: LG9), width and depth of the dentigerous arm of the premaxilla (LG3 and LG6), maxillary head protrusion (LG10), and cranial height (LG1; Table 1, Table 2). The 95% credible interval for dentigerous arm width on LG3 also contained the max LOD score for lower jaw length, suggesting that either a single pleiotropic locus or a cluster of loci in this region may be controlling both traits.

Cranial height
Cranial height was the only trait with statistically significant or marginally significant QTL in both lakes (P < 0.1). While the QTL occurred on different linkage groups between maps, we found a high degree of synteny between these linkage groups indicating that the QTL is located in the same genomic region in both lakes (Table 2, Figure 3). We found 44 genes within scaffold 33 that fell partially or fully within the 95% credible intervals of the QTL in both lakes (Table 1, Error! Reference source not found.). Only three of these genes contained adaptive alleles within 20 kb: wdr31, bri3bp, and gnaq (Table 3). Interestingly, gnaq is well known to be associated with craniofacial development (Hall, Cadle, Morrill-Cornelius, & Bay, 2007;Shirley et al., 2013) and is differentially expressed between our specialist species in developing larvae (McGirr & Martin, 2020).

Dentigerous Arm Width
We found that regions on scaffolds 58 and 24 were associated with a significant QTL for dentigerous arm width in Little Lake and max LOD scores for maxillary head protrusion and female mate preferences in Crescent Pond ( Table 1, Table 2). We found 161 genes which fell partially or completely within these shared regions, but only 2 genes, dysf and cyp26b1, which contained adaptive alleles within 20Kbp ( Table 3). The dysf gene provides instructions for making a protein called dysferlin, which is found in the sarcolemma membrane that surrounds muscle fibers (J. Liu et al., 1998). This could indicate a role for muscle development in affecting skeletal development of the maxilla and premaxilla.

Dentigerous Arm Depth
The QTL for dentigerous arm depth in Little Lake was associated with LG 6, which corresponds to LG 7 in Crescent Pond, however, no traits from Crescent Pond mapped to this linkage group (Table 2, Figure 3). Instead, dentigerous arm depth in Crescent Pond was associated with LG 13, and did not share any similar genomic regions with those associated with dentigerous arm depth in Little Lake. We found 80 genes completely or partially within the 95% credible region for this QTL in Little Lake but none contained adaptive alleles (Error! Reference source not found.).
In fact, only a single adaptive allele was found in this QTL region, but it was in an unannotated region of the genome (Table 3).

Maxillary Head Protrusion
Maxillary head protrusion in Little Lake mapped to a QTL region on LG10 which corresponds to the max LOD scores for both lower jaw length and caudal peduncle height in Crescent Pond (Table 2, Figure 3). Across lakes, all three traits were associated with scaffolds 53, 2336, and 6275. We found 528 genes partially or fully within these shared regions, but only 21 of these genes contained adaptive alleles within 20 Kbp ( Table 3). One of these genes, twist1, contains a non-synonymous substitution fixed in scale-eating pupfish on San Salvador Island, Bahamas (Richards et al. 2021). Twist1 is a transcription factor and oncogene associated with palate development and oral jaw size in model organisms (Parsons et al., 2014;Teng et al., 2018).

Jaw closing In-Lever
The QTL for jaw closing in-lever was associated with LG 9 in Little Lake, which corresponds to the max LOD scores for orbit diameter and anterior body depth in Crescent Pond (Table 2, Figure 3). Scaffolds 8 and 8020 were associated with all three of these traits. We found 13 genes which partially or completely fell within these shared regions, and only two genes, map2k6 and galr2, which contained adaptive alleles within 20Kbp (Table 3). Galr2 was also previously detected within a QTL region for lower jaw length in pupfish .

Dentigerous Arm Depth and Adductor Mandibulae Muscle Mass
Finally, in Crescent Pond the QTL for dentigerous arm depth and adductor mandibulae muscle mass mapped to the exact same location on LG 13 (95% CI dentigerous arm depth (0, 250), adductor mandibulae muscle mass (0,70). This linkage group corresponds to LG14 in Little Lake, which contains the max LOD scores for both palatine height and suspensorium length (Table 3). We found 52 genes that overlapped between these regions, 18 of which contained adaptive alleles. Furthermore, three genes of the genes-ube2w, ncoa2, and prlh-contain adaptive alleles that introgressed from Lake Bavaro in the Dominican Republic to snail-eating pupfish (ube2w), from Lake Cunningham, New Providence Island, to scale-eating pupfish (ncoa2), or from North Carolina, USA, to scale-eating pupfish (prlh). We also found four genes that contained adaptive alleles within 20Kbp that arose from de novo mutations: cd226, cmbl, slc51a, and zfhx, however, only one adaptive allele in slc51a is found within a coding region of the gene.

Parallel genetic changes underlie 5 out of 6 of craniofacial QTL
We found evidence supporting both parallel and non-parallel genetic changes in an adaptive radiation of trophic specialist pupfishes. A single significant QTL was associated with cranial height in both lakes and mapped to the same genomic region, suggesting that parallel genetic changes are responsible for variation in this trait in both lakes. On the other hand, significant QTLs were identified for the premaxilla dentigerous arm depth in each lake, but they mapped to different locations, indicating that this trait is associated with non-parallel genetic changes. We found an additional three traits with significant QTLs (dentigerous arm depth, jaw closing inlever, maxillary head protrusion) in the Little Lake population that were not detected in Crescent Pond. However, all genomic regions associated with these traits in Little Lake also mapped to the max LOD score for this same integrated suite of craniofacial traits in Crescent Pond. Therefore, rather than assume independent QTL for each trait, we conservatively conclude that the same genomic regions are being reused in each lake and affect a highly integrated suite of craniofacial traits. Overall, we found that 5 out of the 6 significant QTLs were reused in some way across lakes suggesting that parallel genetic changes underly adaptive phenotypes in the San Salvador Island pupfish radiation.

High level of QTL reuse across ponds
Overall, we found that about 16% (1 out of 6) of the identified QTL regions corresponded to non-parallel changes and 84% (5 out of 6) corresponded to parallel genetic changes-either affecting the same phenotypic trait or a tightly correlated craniofacial trait-across populations.
The presence of both non-parallel and parallel genetic changes leading to convergent phenotypes across ponds, has also been documented in similar QTL studies. For example, Colosimo et al. (2004) investigated the genetic basis of armor plate morphology in two independent threespine stickleback populations and found a single large effect locus on LG 4 in the two populations.
However, they also noted a potential difference in the dominance relationship of alleles across ponds at this location, and found additional differences in modifier QTLs between populations, suggesting that both parallel and non-parallel genetic changes could lead to the loss of armor plating. Similarly, Erickson et al. (2016) found evidence for both parallel (43% of QTL regions overlapped between at least two populations) and non-parallel (57% of QTL regions were found in only a single population) evolution in a QTL study investigating the genetic basis of 36 skeletal phenotypes in three independent threespine stickleback populations. However, our findings suggest that pupfish exhibit a much higher proportion of parallel evolution than previously documented. In fact, Conte et al. (2012) estimated that the probability of convergence via gene reuse is only 32-55% -which is 1.5 to 2.5 times lower than our current findingalthough this may be underestimated (Stern, 2013).
Pupfish may have a higher rate of parallel evolution than other model fish speciation systems for a few reasons. First, the pupfish radiation is exceptionally young, with specialist species diverging less than 10kya (Hagey & Mylroie, 1995;Martin & Wainwright, 2013), and parallel evolution is predicted to be more likely when populations or species have recently diverged (Rosenblum et al., 2014). This may be because recently diverged species are more likely to experience similar environments, have access to similar pools of genetic variation (either due to standing genetic variation or introgression), or similar genetic constraints. Second, the genomic basis of pupfish skeletal traits may be primarily controlled by cis regulatory elements, which evolve more quickly and have less negative pleiotropy which may make them more likely to undergo parallel evolution (Stern & Orgogozo, 2008). A previous study of allelespecific expression in the pupfish system found strong evidence that two cis-regulatory alleles were associated with skeletal development, but trans-acting elements predominated overall .
In part, the increased proportion of parallel evolution estimated in this study results from our relaxed thresholds for detecting and categorizing regions as "parallel". Previous QTL studies have typically searched for evidence of parallel evolution by only looking for one-to-one mapping in which the same genomic regions are associated with the same trait across populations (Colosimo et al., 2004;Conte et al., 2012). While this method provides the most clear-cut examples of parallel evolution, we argue that it vastly underestimates its frequency in nature. For example, this method would not consider reuse of the same genomic regions for integrated morphological traits as parallel evolution, a pattern seen in this study and in Erickson et al. (2016). Furthermore, the strict one-to-one significance method for detecting parallel evolution does not include consideration of the hierarchy and diversity of convergence and parallel evolution, which can span morphological traits, ecotypes, performance, or even fitness (James et al. 2020;Rosenblum et al., 2014;Stern, 2013). Ultimately, we argue that this method of quantifying parallel evolution provides a more wholistic view of the process and better captures the frequency of reuse of adaptive genetic variation in nature.

Few QTL may affect many craniofacial traits
There are several processes that may cause the same genomic regions to be associated with different traits between lakes. First, these genomic regions may be highly pleiotropic and affect several traits simultaneously. For example, Albert et al. (2007) found that that on average a single QTL affected 3.5 phenotypic traits in an analysis of 54 body traits in three-spine stickleback. Wagner et al. (2008) found a similar pattern in QTL analyses of 70 skeletal traits in mice, where a single QTL affected on average 7.8 phenotypic traits (the maximum being 30).
A second possibility is that a single QTL region may contain several tightly linked causative variants that are responsible for variation in many traits. Correlated phenotypic traits are generally assumed to have a shared genetic basis, but it is extremely difficult to determine if this is due to pleiotropy or tight linkage between genomic regions (Lynch and Walsh 1998;Gardner and Latta 2007;Paaby and Rockman 2013;Wright et al. 2010).
Finally, it may be that differences in methodology or sample sizes between lakes enable us to detect significant QTL for some traits in one lake and not the other. For example, our analyses of Little Lake allowed us to detect significant QTL for effect sizes greater than 6.54 PVE at 80% power, but we could only detect significant QTL for effect sizes greater than 8.41 PVE at 80% power in Crescent Pond due to our lower sample size for this cross (Sen, Satagopan, Broman, & Churchill, 2007). However, this level of power is typical in many non-model QTL studies (Ashton, Ritchie, & Wellenreuther, 2017). The ability to detect a significant QTL in one lake but not the other may be further explained by our use of different sequencing methods between populations. However, a critical component of our analyses involved searching for regions within 10Kbp of one another across maps to provide confidence that if we detected a significant QTL in one lake and not the other that it was not simply because that genomic region was not captured by the sequencing. For example, in Little Lake we detected a significant QTL associated with dentigerous arm depth on LG 6 but did not find any traits associated with this region of the genome in Crescent Pond.

QTL are associated with different traits across different lakes
In this study we found an intriguing pattern of different traits mapping to the same region of the genome across lake populations. One potential explanation for this is that there are different relationships between traits in each lake, and we find some evidence of this in our phenotypic data. Error! Reference source not found. depicts correlation matrices between traits in 1) Little Lake and 2) Crescent Pond, and

Identifying causative regions within QTL
Multiple mapping populations across lakes may also be particularly useful for identifying candidate causal alleles. We found that one out of our six unique QTL regions mapped to the same genomic location across lakes and was associated with the same phenotypic trait-cranial height. In Crescent Pond, we found that a region of 110cM was associated with this trait (LG10, position: 204, 95% CI (130,340)), which contained 426 genes. However, when we compared this region to the region independently identified by our Little Lake analysis, we found that the overlapping region was reduced to 20cM (LG1, position: 259, 95% CI (250-270)) and contained only 44 genes-a reduction of more than 80%. We found a similar pattern in the additional four QTL regions that mapped to the same genomic location across maps but were associated with different phenotypic traits and observed an average 56% reduction in region size. As noted above, Erickson et al. (2016) used a similar method of identifying candidate QTL regions across three hybrid populations of stickleback, and found that 43% of identified QTL regions were shared across two or more populations; however, they did not investigate whether these QTL regions completely or partially overlapped.
We also searched for adaptive alleles within QTL region that were identified in a previous study as 1) nearly fixed between species (F st >0.95) and 2) showed significant evidence of a hard selective sweep (Richards et al. 2021). Overall, we found 789 shared genes within shared QTL regions across lakes, and that 45 of those genes contained adaptive variants (5.7%).
This is a six-fold increase from the genome-wide expectation of 0.91% (176 genes associated with at least one adaptive variant / 19304 annotated genic regions), suggesting that these specific regions are important for adaptation to scale-and snail-feeding in wild pupfish. For example, a variant in twist1 was found within the region associated with maxillary head protrusion in Little Lake (which also overlapped with lower jaw length and caudal peduncle height in Crescent Pond). In model organisms, twist1 is associated with palate and jaw development (Parsons et al., 2014;Teng et al., 2018), and previous genome-wide association scans in pupfish showed that a region containing twist1 was significantly associated with oral jaw size in the system (Richards et al., 2021). Similarly, we found that variants associated with galr2 fall within the QTL region associated with jaw closing in-lever in Little Lake (which also overlapped with regions associated with orbit diameter and anterior body depth in Crescent Pond; scaffolds 8 and 8020), and previous QTL mapping studies, gene expression studies, and genome-wide association analyses have all supported galr2's importance for jaw development in pupfish McGirr & Martin, 2016;Richards et al., 2021).

Increased use of introgressed adaptive variants in QTL regions
We found that most genetic variation within candidate parallel regions is standing across the Carribbean (86.04% within scale-eater populations, and 53.32% within snail-eating populations).
Furthermore, we found more introgressed adaptive variants in both scale-eater (observed: 12.13% introgressed, expected 95% CI: (7.96%-9.88%)) and snail-eater populations than expected by chance (observed: 46.67%, expected 95% CI: (32.22%-37.06%)). This supports the prediction that standing genetic variation and/or introgressed variants should underlie parallel genetic changes (Stern, 2013;Thompson, Osmond, & Schluter, 2019). Finally, we found that only 1.83% of adaptive variants associated with traits across both lakes originated from de novo mutations on San Salvador Island. While this percentage did not differ significantly from the expected estimates (expected 95% CI: 1.3%-2.17%) it does not eliminate the possibility that de novo mutations play an important adaptive role in pupfish evolution. For instance, a single de novo mutation fell within the QTL region of dentigerous arm depth and is within an exon of twist1. Future work should functionally test the effects of implicated regions to confirm their phenotypic effects.

Conclusion
In conclusion, we found that a single QTL region was responsible for variation in cranial height in both populations, and an additional four QTL regions were responsible for variation in different craniofacial traits across lakes, suggesting that parallel genetic changes underlie integrated suites of adaptive craniofacial phenotypes on San Salvador Island. Adaptive alleles were more commonly found within these detected QTL regions, and more of these adaptive alleles arrived on SSI via introgression than expected by chance. Finally, we argue that investigating QTL regions across populations in concert with estimation of adaptive alleles in the wild is a powerful tool for identifying potential causative regions of the genome affecting adaptive divergence. Table 1. All phenotypic traits measured in Little Lake and Crescent Pond. A genome scan with a single-QTL model by Haley-Knott regression, was used to identify the position with the highest LOD score, the 95% bayes credible intervals, and the region's significance level for each trait (traits significant at the P<0.1 level are identified with a circle and those significant at the P<0.05 level are identified with an asterisk). Additionally, we report the scaffold numbers of genomic regions that fell within max LOD interval for each trait, the number of individuals phenotyped, percent variance explained (PVE) by a max LOD region, and P-value associated with each max LOD region.   Table 3. Count of adaptive alleles within genomic regions associated with positions of max LOD in both lakes. Adaptive alleles were categorized as either standing genetic variation (SGV), introgression (Intro.), or de novo mutations, and are estimated independently for snail-eaters and scale-eaters. Asterisks represent traits that were significant at least at the P<0.1 level in the genome scan, while crosses show traits that corresponded to the same locations in the alternate lake.

Figure 2
Cyprinodon linkage maps for A) Crescent Pond and B) Little Lake. The Crescent Pond Linkage map was estimated from 743 markers, and the Little Lake linkage map was estimated from 540 markers. Both maps were generated from F2 crosses between a scale-eater (C. desquamator) and snail-eater (C. brontotheroides) from the respective lakes.
d Figure 3 Circos plot depicting the relationship between the Crescent Pond (red) an linkage maps (blue), which share 324 markers within 10Kbp of one another. Numb surrounding each semi-circle represent linkage group numbers in each lake. Marke shared across lakes are connected via the colored lines.
and Little La mbers rkers that are