Local adaptation does not lead to genome‐wide differentiation in lava flow lizards

Abstract Adaptation can occur with or without genome‐wide differentiation. If adaptive loci are linked to traits involved in reproductive isolation, genome‐wide divergence is likely, and speciation is possible. However, adaptation can also lead to phenotypic differentiation without genome‐wide divergence if levels of ongoing gene flow are high. Here, we use the replicated occurrence of melanism in lava flow lizards to assess the relationship between local adaptation and genome‐wide differentiation. We compare patterns of phenotypic and genomic divergence among lava flow and nonlava populations for three lizard species and three lava flows in the Chihuahuan Desert. We find that local phenotypic adaptation (melanism) is not typically accompanied by genome‐wide differentiation. Specifically, lava populations do not generally exhibit greater divergence from nonlava populations than expected by geography alone, regardless of whether the lava formation is 5,000 or 760,000 years old. We also infer that gene flow between lava and nonlava populations is ongoing in all lava populations surveyed. Recent work in the isolation by environment and ecological speciation literature suggests that environmentally driven genome‐wide differentiation is common in nature. However, local adaptation may often simply be local adaptation rather than an early stage of ecological speciation.

Of course, adaptation does not always lead to genome-wide differentiation, especially when locally adaptive phenotypes have a simple genetic basis and no consequences for reproductive isolation.
While examples of IBA and IBE are increasingly common in the literature (Funk et al., 2011;Orsini, Vanoverbeke, Swillen, Mergeay, & Meester, 2013;Seehausen et al., 2008;Sexton et al., 2014;Shafer & Wolf, 2013), some authors suggest that there may be a publication bias against studies that do not find a correlation between genetic differentiation and environmental differences (Feder, Egan, & Nosil, 2012;Shafer & Wolf, 2013). Thus, it is possible that cases where patterns of genomic variation across environmental gradients do not follow expectations of IBE or IBA are underreported in the literature. Thus, replicated studies across species and evolutionary timescales are needed to understand the frequency with which-and the conditions under which-local adaptation does or does not lead to genome-wide differentiation in the wild.
The lava flows of the Chihuahuan Desert in New Mexico represent an ideal place to study whether patterns of IBA and IBE are shared across evolutionary and environmental replicates because multiple species are found on multiple independent lava flows. The dark basalt rocks of the lava flows contrast with the surrounding adobe colored soils, and many vertebrates exhibit darker-melanistic (Majerus, 1998)-color morphs on the lava flows (Best, James, & Best, 1983;Lewis, 1949Lewis, , 1951Rosenblum, 2005). Given that numerous lava flows are within approximately 100 miles of each other in the Chihuahuan Desert, they draw from a similar pool of colonists and represent multiple replicates of colonization, adaptation, and community assembly. The lava flows vary in age and substrate color, representing snapshots of evolution since initial colonization. The youngest lava flow is only 5,000 years old and has more homogeneously black substrate (Dunbar, 1999). The more ancient lava flows are orders of magnitude older (106,000 and 760,00 years old) and have more sand accumulation creating a more heterogeneously colored background (Bachman & Mehnert, 1978;Hoffer, 1975).
A number of reptile species have been reported to be significantly darker on lava flows in the Chihuahuan Desert than in the surrounding desert (Krohn & Rosenblum, 2016;Lewis, 1951;Rosenblum, 2005). Melanistic phenotypes in lava flow lizards are typically considered to be an adaptation for crypsis because visually hunting diurnal predators are known to prey on color-mismatched lizards (Luke, 1989). Melanistic lava flow lizards of the focal species typically match their background extremely well, often better than lizards from surrounding nonlava flow populations (Krohn & Rosenblum, 2016;Rosenblum, 2005). Dark coloration in lava flow lizards also appears to be heritable as the degree of melanism does not change ontogenetically, is not influenced by the juvenile's substrate color, and correlates strongly with parental coloration (Corl, Lancaster, & Sinervo, 2012;Micheletti, Parra, & Routman, 2012;Rosenblum, 2005). Additionally, while lizards can change color over short time periods, physiological color change is insufficient to explain the observed differences in coloration in these species (Krohn & Rosenblum, 2016;Micheletti et al., 2012;Rosenblum, 2005).
Although melanism appears to be heritable and associated with environmental variation in many reptile species, it is not clear whether local adaptation in coloration is expected to lead to genomewide differentiation. On one hand, patterns may simply reflect migration-selection balance, where selection on relatively few genes is balanced by ongoing gene flow and genome-wide divergence is unlikely (Corl et al., 2018;Rosenblum, Hickerson, & Moritz, 2007). On the other hand, color variation in reptiles can have consequences for both natural and sexual selection as many lizard species use visual cues including coloration for mate choice (Hardwick, Robertson, & Rosenblum, 2013;Robertson & Rosenblum, 2009). In this case, local adaptation could lead to genome-wide divergence and ultimately reproductive isolation. Thus, empirical tests are necessary to determine whether local adaptation typically leads to genomic differentiation between lava flow populations and nonlava flow lizards.
Here, we focus on understanding patterns of local adaptation and genomic differentiation in three lizard species distributed across three lava flows in the Chihuahuan Desert. The three focal species-Crotaphytus collaris, Sceloporus cowlesi, and Urosaurus ornatus-have all been reported to have melanic populations on local lava flows (Krohn & Rosenblum, 2016;Lewis, 1951;Rosenblum, 2005). The three lava flows differ in age-from 5,000 to 106,000 to 760,000 years old-providing an opportunity to evaluate whether divergence accumulates (or decays) over time. For populations on and off the lava flows, we quantified color variation, measured genetic variation using restriction-associated DNA sequencing (Peterson, Weber, Kay, Fisher, & Hoekstra, 2012), and estimated migration rates using demographic simulations. Our sampling approach leveraged population trios (where we compared each focal lava flow population with two nearby nonlava flow populations) to disentangle the effects of IBD and IBE.

| Sampling and color quantification
We sampled five populations of three species of lizard across three lava flows (Figure 1). We caught lizards from the Carrizozo Lava Flow (CZLF; ~5,200 years old (Dunbar, 1999)), Aden Afton Lava Flow (~106,000 years old, (Hoffer, 1975)), and Pedro Armendariz Lava Flow (PALF; ~760,000 years old (Bachman & Mehnert, 1978)) in south-central New Mexico. Our sampling scheme focused on population trios: for each lava flow lizard population, we also caught lizards from two nonlava flow populations (Figure 1). The nonlava flow populations were as geographically close to the lava flows as possible. For clarity, each population trio will be referred to by an abbreviated version of the lava flow name and a species. All sampling occurred between June and August, 2013 and 2015.
At CZLF, we quantified the dorsal coloration of 18 C. collaris (CZLF-Crotaphytus) and 26 S. cowlesi (CZLF-Sceloporus). We included additional phenotypic data of S. cowlesi at CZLF from a previous study (Rosenblum, 2006). Surrounding the CZLF (non-  ; Table S1).  For all populations' trios except the CZLF-Sceloporus population trio, we used photographs to quantify dorsal coloration. We followed the procedures outlined in Krohn and Rosenblum (2016) to take standardized photographs of the lizards, using the same camera and camera settings as in Krohn and Rosenblum (2016). Briefly, we photographed lizards at their mean operating body temperature (C. collaris at 36°C, S. cowlesi at 35°C and U. ornatus at 35°C) according to Brattstrom (1965). We photographed all lizards on white paper in an upside-down bucket in direct sunlight to allow for a diffuse, constant light source for all photographs. All photographs included a white, gray, and black standard (Adorama QPcard 101) for color standardization. We imported photographs in RAW format, linearized the photographs using the gray and white standards in Photoshop, and then exported the photographs in the lossless TIF format (Krohn & Rosenblum, 2016). To quantify overall dorsal coloration, we took the average R, G, and B values across a 75 × 75 pixel square in C. collaris (Krohn & Rosenblum, 2016), across the dorsal midline stripe in S. cowlesi (Rosenblum, 2006), and across the entire dorsal body surface of U. ornatus in ImageJ using the plugin RGB measure (v1.47; National Institutes of Health). To convert the RGB values to a color space meaningful for most vertebrate visual systems, we used properties of color and brightness opponency (Endler, 2012;Endler & Mielke, 2005;Osorio, Vorobyev, & Jones, 1999) to calculate values of hue, chroma, and luminance from the RGB values as in Krohn and Rosenblum (2016). Given that luminance explains significant differences in coloration between lava flow and nonlava flow lizards (Krohn & Rosenblum, 2016;Rosenblum, 2006), we used luminance to compare lizard coloration in all our subsequent analyses.
We quantified coloration of the CZLF-Sceloporus population trio using data collected for a prior study using a spectrophotometer. For CZLF-Sceloporus and RI, we integrated the raw reflectance data from Rosenblum (2006). Because lizards from TU were collected after the Rosenblum (2006) study, we collected and analyzed reflectance data identically, but used an Ocean Optics Jaz spectrophotometer.
As above, and in Rosenblum (2006)

| Sampling and RAD sequencing
We

| RAD sequence processing and filtering
Before analyzing our data, we filtered our reads to ensure only highquality loci. We followed the RADToolKit pipeline (https ://github. com/CGRL-QB3-UCBer keley/ RAD) to demultiplex individuals, remove low-quality reads, create individual contigs, and to create a de novo reference assembly. We demultiplexed individuals while allowing for one mismatch per barcode. After demultiplexing, we separated individuals by species and continued the pipeline on each species independently. We filtered out reads that were low quality (Phred score <20), were contaminated by bacteria, had anomalous runs of a single base (more than 50% of the read), or did not contain the restriction enzyme cut site (allowing a mismatch of 1 bp). After filtering, we removed individuals that did not have enough reads to continue (fewer than 250,000 reads per individual; n = 14 C. collaris, n = 9 S. cowlesi, n = 2 U. ornatus). Next, we masked repetitive elements known from vertebrate metazoans using RepeatMasker (version 4.0, 2013-2015) and clustered reads into contigs for each individual using CD-hit (Fu, Niu, Zhu, Wu, & Li, 2012;Li & Godzik, 2006), with a sequence identity threshold between 0.8 and 0.95 (optimized for each species using the program "cluster test" from RADToolKit

| Data analysis
First, to understand basic population structure, we visualized genetic data for each species using a principal components analysis (PCA). We used the program ngsCovar to create a covariance matrix for each individual's genotypes (Fumagalli et al.., 2013). ngsCovar accounts for the uncertainty of each genotype by incorporating genotype posterior probabilities to generate the covariance matrix.
For each trio of populations, we input genotype posterior probabilities from ANGSD. We called genotype posterior probabilities (not including genotypes with a posterior <0.95) at sites that, in addition to being present after the filtering described above, had a minimum coverage of 5X, were polymorphic (using a likelihood ratio test with a minimum p-value of 1 × 10 −2 ), had a minor allele frequency >0.05, and were present in at least 80% of individuals. We plotted the eigenvalues of the covariance matrix in R to visualize population structure in each of the population trios.
Second, to understand the relationship between genetic variation and geographic distances, we calculated a pairwise genetic distance matrix for all individuals in each of the population trios. We used ngsDist to calculate genetic distance between all individuals (Vieira, Lassalle, Korneliussen, & Fumagalli, 2016). ngsDist takes genotype uncertainty into consideration by calculating genetic distances from genotype likelihoods or posterior probabilities instead of directly from the genotype calls. Using the same filters as above,   yses, but using -doGeno 2 in ANGSD to output allele calls. We converted the allele calls into the ∂a∂i input format using custom Python scripts (available at https ://github.com/alexk rohn/LavaF lowLi zards/ ).
Using ∂a∂i, we created joint site frequency spectra for each pair of populations. We ran four models of demographic scenarios to model how alleles might change with the formation of the lava flow and ongoing migration. We ran a neutral model with no divergence, a model with divergence but no migration, a model with divergence and symmetric present-day migration, and a model with divergence and asymmetric present-day migration (as in Portik et al., 2017). For each population comparison and each model, we followed the procedure of Portik et al. (2017) to optimize parameters. We optimized our original parameters over three rounds. To get the best final estimates of our parameters, each round of optimization had an increasing amount of iterations of parameter optimization, an increasing amount of total replicates from modified starting parameters, and a decreasing degree to which parameters were perturbed (Portik et al., 2017). From each round, we took the parameters with the lowest AIC score to continue to the next round of optimization (Ese et al., 2016). The first round had 20 iterations over 50 replicates, the second round had 30 iterations over 50 replicates, and the final round had 50 iterations over 100 replicates. We determined that a model best described the data when its AIC score was at least 10 less than the next best model (Ese et al., 2016). Scripts of models used and optimization protocols are available online (https ://github.com/dport ik/dadi_pipel ine/tree/maste r/Two_Popul ations).

| RE SULTS
Although we observed strong geographic color variation among   Figure 1. While there is at least some population-level clustering in all three species, lava populations are not more distinct than expected based on background divergence. Figure 2 and  Figure 3).
Crotaphytus was the only species in which genetic distance was not correlated with geographic distance.
Partial Mantel tests confirmed that genetic differentiation could not be explained by environmental differences, even when correcting for geographic distances (Table 1). Only CZLF-Crotaphytus showed an increase in genetic distance with environmental distance when con-  (Table 2). Demographic analyses using ∂a∂i revealed that migration onto and off of the lava flows is likely relatively common. All species showed evidence of migration between at least one lava flow and nonlava flow population (Table S2). There was no consistent pattern observed in terms of patterns of migration and lava flow age or background color homogeneity. In U. ornatus, the best model predicted asymmetric migration with stronger migration from nonlava flow populations to the lava flow populations. In C. collaris, the best model indicated symmetrical migration between the PALF-Crotaphytus and one nonlava flow population. For the other lava flow in C. collaris, the best model suggested relatively high migration between CZLF-Crotaphytus and one nonlava flow population (OS), but AIC values were too close (>10) to determine whether migration was symmetric or asymmetric.
Finally, in S. cowlesi, the best model predicted migration between PALF-Sceloporus and one nonlava flow population (PALF -MC) but could not distinguish whether migration is symmetric or asymmetric.
Thus, regardless of levels of phenotypic divergence between lava flow and nonlava flow populations, migration is ongoing across populations within species.

| D ISCUSS I ON
We demonstrated that genome-wide differentiation does not accompany local adaptation in three species of lava flow lizards in the Chihuahua Desert. Consistent with previous literature (Lewis, 1949(Lewis, , 1951, we found that C. collaris and S. cowlesi from the CZLF (~5,200 years old) and U. ornatus from the Aden Afton Lava Flow (~106,000 years old) were melanistic, indicating strong phenotypic divergence from nonlava flow populations. Despite what was reported in earlier literature (Best et al., 1983), C. collaris and S. cowlesi from the PALF (~760,000 years old) were not significantly darker than surrounding nonlava flow populations. The phenotypic patterns themselves are intriguing: the only lava flow lizards that were not melanistic were those from the oldest lava flow, where the substrate is less homogenously dark. However, our most important finding was that, irrespective of phenotypic differentiation, the lava populations typically did not exhibit more divergence from nonlava conspecifics than expected based on geography. Most comparisons also revealed that migration occurs between lava flow and nonlava flow populations. Thus, we find some evidence of local adaptation, but little evidence of the accompanying genetic isolation often reported in the IBE and IBA literature.
In our replicated comparisons, we observed a lack of strong environmentally driven genome-wide differentiation, regardless of the geological age of the lava flow. The three lava flows differ in age, substrate color, and habitat uniformity, but populations typically showed patterns more consistent with IBD rather than IBE. If there was habitat-specific genetic differentiation, we would expect greater genetic divergence between lava and nonlava populations (compared with background nonlava to nonlava comparisons). In most cases, overall divergence between lava and nonlava populations was low (F ST < 0.15; Hartle & Clark, 1997;Frankham, Ballou, & Briscoe, 2002 nonlava populations was comparable to background differentiation. Moreover, we did not uncover a signal of genetic divergence accumulating or decaying with lava flow age. For example, in Crotaphytus, F ST between lava and nonlava populations was comparably low regardless of whether the lava flow was 5,200 or 760,000 years old. Using partial Mantel tests, only one comparison (CZLF-C. collaris) showed significant genetic differentiation due to the environment, and even in this comparison, the magnitude of genetic differentiation between habitats was small ( Figure 3 and Table 2). Thus, when we detected any genetic differentiation, it appears, as found in some other recent studies (Aguillon et al., 2017;Heath, Schrey, Ashton, Mushinsky, & McCoy, 2012;Wang, Glor, & Losos, 2013), that geographic factors such as distance are more important in shaping genetic patterns than environmental factors or local adaptation (Figure 3).
One explanation for lack of genome-wide differentiation is ongoing migration between lava flow and nonlava flow populations.
Shallow divergence times and ongoing migration in geologically novel habitats (like recently formed lava flows) can make it difficult to distinguish between demographic scenarios, even with thousands of loci (Pinho & Hey, 2010). Yet, for three of the four population trios tested, we were able to detect recent migration between at Melanism has a genetic basis in many animal species (van't Hof, 2011;  phenotypic plasticity contributes to melanism in these lava flow lizards (Corl et al., 2018;Luke, 1989). Previously, we demonstrated that physiological color change plays a role in the coloration of C. collaris at PALF (Krohn & Rosenblum, 2016). However, we have no evidence for developmental plasticity, and other studies demonstrate that melanism is heritable in other lava flow lizards (Corl et al., 2012;Micheletti et al., 2012;Rosenblum, 2005). If melanism has a genetic basis in the species studied here, it is likely few genes control coloration, such that relatively little of the genome is under divergent selection (Gagnaire, Pavey, Normandeau, & Bernatchez, 2013;Nachman et al., 2003;Nosil, Funk, & Ortiz-Barrientos, 2009;Rosenblum, Rompler, Schoneberg, & Hoekstra, 2010). Additionally, if adaptive loci are not linked to traits involved in reproductive isolation, genome-wide differentiation is less likely to occur (Feder et al., 2012;MacNair & Christie, 1983;Presgraves, 2013;Slatkin, 1987). Interestingly, a previous study in this system found no preference for local mates in either lava flow or nonlava S. cowlesi (Hardwick et al., 2013). Thus, there is no evidence for preferential mate choice or reproductive isolation across lava and nonlava S. cowlesi populations, despite phenotypic differentiation. Without a link between the genetic mechanisms of melanism and reproductive isolation, gene flow could easily persist even with selection for background-matched lizards.
Our replicated results show that regardless of the species studied or the age of the lava flow, locally adapted populations show ongoing gene flow and no more genome-wide differentiation than expected by geographic distance. Clearly, adaptation does not always lead to genome-wide divergence, or speciation. Given the similarity of our results across populations, it is possible that selection-migration balance without genome-wide differentiation may be very common (Hendry, 2009). However, given the claim that incipient ecological speciation is widespread, that IBE is common, and that studies that fail to find large effects of the environment on genetic differentiation are seldom published (Feder et al., 2012;Hendry, 2009;Sexton et al., 2014;Shafer & Wolf, 2013), it is important to also highlight natural systems where local adaptation occurs without genome-wide differentiation. Animal Care and Use protocol R347-0314.

CO N FLI C T O F I NTE R E S T S
None declared.

AUTH O R CO NTR I B UTI O N S
A.R. Krohn designed the experiment, conducted the field and laboratory work, analyzed the data, and wrote the manuscript. E.T.
Diepeveen contributed to fieldwork and laboratory work. K. Bi contributed to data analysis. E.B. Rosenblum contributed to the experimental design and manuscript preparation.

DATA AVA I L A B I L I T Y
Raw sequence data are archived in the Short Read Archive