Microgeographic morphological variation across larval wood frog populations associated with environment despite gene flow

Abstract Gene flow has historically been thought to constrain local adaptation; yet, recent research suggests that populations can diverge despite exchanging genes. Here I use a common garden experiment to assess the combined effects of gene flow and natural selection on morphological variation of 16 wood frog (Rana sylvatica) populations, a species known to experience divergent selection pressures in open‐ and closed‐canopy ponds across relatively small geographic scales. Wood frog tadpoles from different ponds showed significant morphological variation associated with canopy type with a trade‐off between tail length and body depth consistent with previous research. In contrast, neutral genetic differentiation of nine microsatellite loci as measured by Jost's D was not associated with canopy type, indicating no pattern of isolation by environment. Genetic structure analyses indicated some substructure across the 16 ponds (K = 4); however, three out of four assigned clusters included both open‐ and closed‐canopy ponds. Together, these results suggest that morphological divergence among these wood frog populations is occurring despite gene flow and that selection within these environments is strong. Furthermore, morphological variation among ponds differed across two sampling periods during larval development, demonstrating the importance of evaluating phenotypic divergence over multiple time periods and at a time relevant to the processes being studied.


| 2505
ZELLMER to which it is affected by the interplay between gene flow and selection. Genomes are heterogeneous (Wu, 2001), with the strength of selection and rate of introgression varying across different genes. As a result, genes coding for phenotypic traits under intense selection pressures may diverge between populations, whereas nonselected genes become homogenized (Nosil, Egan, & Funk, 2008). Understanding how this interplay between gene flow and selection affects adaptation is crucial, as the extent to which populations can become locally adapted and phenotypically divergent is important for understanding the processes of divergence and speciation. Further, knowing the relative roles of gene flow and selection is necessary for predicting the potential for populations to adapt and cope with human-induced environmental changes following reduced gene flow due to habitat loss and fragmentation.
In this study, I sought to assess the contribution of gene flow and selection to local adaptation of the wood frog, Rana sylvatica, a North American amphibian that inhabits a relatively broad environmental gradient, with larvae occupying both open-and closed-canopy ponds (Werner & Glennemeier, 1999). Open-canopy ponds on average have greater resource availability, have higher dissolved oxygen levels, have longer hydroperiods, are warmer than closed-canopy ponds (Werner & Glennemeier, 1999), and also harbor more predators (Relyea, 2002b), whereas closed-canopy ponds-due to resource scarcity-have higher levels of intraspecific competition (Werner, Skelly, Relyea, & Yurewicz, 2007). As a result, selection pressures in open-versus closed-canopy ponds are strongly divergent. Selection by predators favors individuals that invest in antipredator defenses, including reduced activity and increased tail development, at the expense of decreased growth rates (Relyea, 2000(Relyea, , 2001aRelyea & Werner, 2000;Van Buskirk, 2002;Van Buskirk, McCollum, & Werner, 1997;Van Buskirk & Relyea, 1998). In contrast, intense competition (or low resource levels) favors individuals that maximize growth rates as opposed to antipredator defenses (Relyea, 2002a). In fact, previous research has demonstrated distinct morphological differences among ponds in common garden and reciprocal transplant experiments, suggesting that these selection pressures promote the adaptation of wood frog populations to local environmental conditions (Relyea, 2002b;Skelly, 2004).
While these population-level differences were initially attributed to isolation of ponds due to limited dispersal (Relyea, 2002b), more recent research suggests that amphibians have greater dispersal capabilities than previously thought (Smith & Green, 2005). Since the alternative habitats are interspersed across the landscape, with open-and closed-canopy ponds well within the dispersal capabilities of wood frogs (Berven & Grudzien, 1990), there is potential for gene flow to occur among opposing selective regimes. Thus, the phenotypic differences among ponds may be due to either environmental barriers preventing gene flow among populations (Wang & Bradburd, 2014) or strong selection on maladapted individuals (Relyea, 2002b). Although local adaptation and phenotypic variation of wood frogs in response to environmental variation have long been studied through both experimental research and natural history observations (Relyea, 2000(Relyea, , 2001aRelyea & Werner, 2000;Van Buskirk, 2002;Van Buskirk & Relyea, 1998;Van Buskirk et al., 1997), the genetics of this variation and contribution of gene flow to these phenotypic differences has yet to be assessed.
To determine the relative roles of gene flow and local environmental conditions on phenotypic variation, I conducted a common garden experiment paired with an evaluation of population genetic structure of wood frog larvae from 16 ponds with divergent environments. Since environmental differences are thought to promote local adaptation, I expect morphological differences between open-and closed-canopy ponds, consistent with previous research (Relyea, 2002b). However, this variation may also be impacted by the amount of gene flow occurring among ponds. I thus investigated the extent to which gene flow occurs among ponds and among canopy types and tested for an association between gene flow and the extent of morphological differences among ponds. If morphological variation is the result of strong selective differences among ponds despite gene flow, then there should be a significant difference in morphology among openand closed-canopy ponds with evidence of gene flow. If on the other hand selection influences gene flow, then I expect a significant pattern of isolation by environment, with neutral genetic divergence increasing with environmental differences. I evaluated these hypotheses at two separate time periods during larval development, while larvae are more or less susceptible to predation by gape-limited predators, to determine whether the effects were lasting beyond the selection period. The traits evaluated include mass and five morphological measurements: tail length and body depth, which are heritable and under selection by predators (Relyea, 2005;Van Buskirk & Relyea, 1998), tail depth and body length, which have heritable plasticity (Relyea, 2005), and muscle depth, which is heritable (Relyea, 2005)

| METHODS
To test hypotheses about the role of gene flow and environment on morphological divergence, I paired a common garden mesocosm experiment with an analysis of environmental and genetic variation of the original populations. I conducted a common garden experiment to quantify the morphological variation among ponds independent of environmental influences. To quantify environmental differences, I assessed pond canopy cover of each of the populations. To assess genetic divergence, I used microsatellite data to quantify pairwise genetic differences as well as population genetic structure. Finally, I compared morphological differences among ponds in relation to canopy and genetic cluster.

| Study species
The wood frog (Rana sylvatica) is a pond-breeding amphibian that inhabits much of eastern North America, with a range extending up to Alaska. The wood frog is an explosive breeder that lays eggs during a 1-2 week period in the spring. Females usually lay their eggs in a single egg mass consisting of on average 711 eggs (Benard, 2015), and most females from a single breeding chorus lay their egg masses next to one another. Breeding populations are therefore usually discrete units. Eggs take approximate 1-2 weeks to hatch. Tadpoles grow in the ponds for about 6 weeks and then after metamorphosis, enter forested habitat.

| Common garden experiment
To quantify morphological differences among wood frog populations, I conducted a common garden experiment raising individuals from 16 ponds that were classified as either open-or closed-canopy in mesocosms with similar conditions. Ponds were selected to ensure that dispersal and gene flow were possible among populations. Each of the populations used in this study were within the average wood frog dispersal distance (approximately 1.2 km: Berven & Grudzien, 1990) from at least one other pond and were either within or (in one case) just beyond the maximum-recorded dispersal distance (2.53 km: Berven & Grudzien, 1990) from at least one other pond of the opposite canopy type ( Figure 1). Moreover, some ponds were within less than 165 m from a pond of the opposite habitat type.
During the 2008 breeding season, 16 ponds were visited routinely to determine breeding chorus locations. Approximately 100 eggs were collected from each of 10 egg masses from each pond to ensure a broad sampling of the population. For one pond (Cassidy 1), egg masses could not initially be located, and instead approximately 15 amplectant pairs were caught, kept separate, and returned to the laboratory to breed. Eggs from 10 masses (laid within 24 hours of collection) were kept for this study, and the adults were returned to their pond of origin. Egg masses were later located to confirm that breeding did occur in this pond. All eggs were kept until hatching in outdoor wading pools covered by shade cloth. Individuals were raised in common garden mesocosms (1000 L polyethylene cattle watering tanks). Each population was raised in separate tanks replicated four times across four spatial blocks, for a total of 64 tanks (16 populations × 4 blocks). The mesocosms were set up 16-21 April, 2008. Each tank was filled with aged well water, inoculated with zooplankton and approximately 6 L of filtered pond water to initiate phytoplankton growth, supplemented with approximately 300 g of leaves to serve as a substrate for phytoplankton, and cov- then preserved in 10% formalin for morphological measurements. The experiment was initiated on 23 April, with 100 hatchlings added to each tank.

| Morphological measurements
On days 18 and 37, ten tadpoles were removed from each tank and preserved in 10% formalin, for a total of 40 individuals per population.
The sampling dates were chosen for comparison to previous studies (Relyea, 2002b) and to evaluate changes over time. At day 18, tadpoles were still small enough to be susceptible to gape-limited predation, whereas at day 37, tadpoles should be large enough to be outside the limits of many predators. It is important to note that they may still be susceptible to gape-unlimited predators. All tadpoles collected were photographed and weighed. Morphological measurements were made on the photographs using ImageJ software. Since the effects of selection and gene flow can vary across the genome, phenotypic measurements were made on a number of traits that are heritable (Relyea, 2005) and under selection by aquatic predators (Van Buskirk & Relyea, 1998). Five morphological measurements were made on each individual, including body length, tail length, body depth, tail depth, and muscle depth (see Relyea, 2000). While geometric morphometric methods provide additional information over linear measurements (Rohlf & Marcus, 1993), linear measurements were used because previous studies have identified adaptive differences based on linear measurements in wood frog larvae. Consistency with these previous studies allows for direct tests of hypotheses generated from previous work. All phenotypic measurements were conducted without knowledge of the source population to reduce the possibility for bias in the measurements.
The morphological measurements were ln-transformed for linearity. All statistical analyses were conducted in R (R Core team 2015).
After averaging across individuals within tanks for each time period, I conducted a nested, repeated-measures ANCOVA to test for significant variation for each trait across ponds and over time taking into account any differences among spatial blocks. I included ln-transformed mass to control for differences due to body size. P-values were corrected using a sequential Bonferroni correction (Rice, 1989). To facilitate pairwise comparisons, morphological measurements were also averaged across tanks to a single value for each pond at each sampling period. These pond averages were used for all subsequent analyses.
To quantify a combined measure of morphological differences among ponds, I created a morphological distance matrix using the R "vegan" package (Oksanen, Blanchet, Kindt, Legendre, & O'Hara, 2016) for each sampling period.

| Genetic divergence
To quantify the extent of gene flow among ponds, I evaluated genetic population structure among the 16 ponds by using a clustering analysis to detect genetic clusters and by calculating pairwise genetic differentiation. Both measures of genetic differentiation were based on data from nine microsatellite loci from approximately 20 individuals per population published in a previous study (Zellmer & Knowles, 2009 iterations to assure alpha converged in each run. For each set of analyses, K was set from 1-17 and was replicated three times for each K. Following the STRUCTURE analyses, I evaluated the data using Structure Harvester (Earl & vonHoldt, 2012) and CLUMPP (Jakobsson & Rosenberg, 2007) via CLUMPAK (Kopelman, Mayzel, Jakobsson, & Rosenberg, 2015). The highest likelihood and the max ΔK value (Evanno, Regnaut, & Goudet, 2005) in addition to histograms of population assignment values were each used to identify the number of clusters with the best support.
Genetic differentiation was calculated as D est (Hedrick, 1999;Jost, 2008) using the R "diveRsity" package (Keenan, Mcginnity, Cross, Crozier, & Prodöhl, 2013), since D est is less susceptible to gene variation resulting in a better estimator of allelic differentiation among populations as compared to G ST values (Jost, 2008). Populations were considered diverged if the 95% confidence intervals for Jost's D did not overlap 0.
To determine whether gene flow is limited by environmental difference, I conducted an Isolation by Environment (IBE) analysis (Wang, 2013 Geographic distances were calculated as Haversine distances using the R "geosphere" package (Karney, 2013). MMRR was conducted in R (Wang, 2013).

| Predictors of morphological divergence
To quantify the relative role of environmental differences and gene flow in phenotypic differentiation among populations, I conducted a PERMANOVA to test whether canopy or gene flow predict pairwise differences in morphology across ponds using the R "vegan" package (Oksanen et al., 2016). The model included canopy as a factor and genetic cluster (as determined by STRUCTURE) as predictor variables in addition to ln-transformed mass as a covariate. The analysis was repeated using the combined morphological distance matrix for each time period (days 18 and 37) and also for each individual trait. p-values were corrected with sequential Bonferroni.
I also conducted a principal components analysis (PCA) to evaluate differences in morphology across canopy types and genetic clusters.
Pond averages for each of the five ln-transformed morphological measurements in addition to ln-transformed mass were used for the

| Predictors of Morphological Variation
Across all morphological traits combined, the differences in morphology across ponds were significantly different across canopy types (PERMANOVA: p < .023) at day 18 only but significantly differed across clusters (p < .001), and was significantly associated with ln-transformed mass (p < .001) on both days (Table 2). There were differences among morphological traits in the extent of association with canopy and cluster (Figures 5-6; Tables S3-S4). At day 18, tail length, body length, and body depth differed across canopy (PERMANOVA: p < .002, p < .001, p < .001, respectively; Table S3).
At day 37, only body length significantly differed across canopy types (p < .001), although tail depth also showed a trend for differences across canopy type (p < .036; Table S4). All morphological variables aside from muscle depth (p > .05) differed across clusters and were significantly associated with mass for both time periods (Tables   S3-S4).
Combining morphological traits in a principal components analysis elucidated overall differences across ponds.  Figure 7). PC2 showed no variation in association with canopy or genetic cluster (p > .05; Table 3). PC3 and PC4 showed a trend toward significant differences across canopy type at day 18 (p = .077, p = .024), although was not significant following sequential Bonferroni correction (Table 3; Figure 7b).

| DISCUSSION
The presence of fine-scale phenotypic variation in wood frog populations has been hypothesized to be due to either population F I G U R E 2 Likelihood and Δ K scores for increasing number of genetic clusters.
Likelihood (a) and ΔK scores (b) for increasing number of clusters (K) using either no prior (open circles, dotted line), individual ponds as a prior (closed circles, solid line), pond environment (closed squares, dashed line) as a prior isolation or due to very strong selection within ponds (Relyea, 2002b).
Determining the relative contribution of selection and isolation to population divergence is important for understanding fine-scale evolutionary processes and ultimately speciation (Ohmer, Robertson, & Zamudio, 2009). Using a common garden experiment, I demonstrated that larval wood frog populations exhibit significant morphological variation across an environmental gradient of canopy cover ( Figure 5; Table 2), consistent with previous research (Relyea, 2002b(Relyea, , 2005. Wood frog tadpoles from different ponds showed significant morphological variation associated with canopy type with a trade-off between tail length and body depth ( Figure 5). However, contrary to expectations, there was no evidence that gene flow is limited across the environmental gradient (Figures 3-4). The morphological differences occur over surprisingly small geographic scales, with gene flow between neighboring ponds of opposite canopy type (Figure 1).
While there was some genetic population structure across the entire study area, there was evidence of high levels of admixture among ponds ( Figure 3). Genetic clusters corresponded to location rather than canopy types (Figures 2-3), with each genetic cluster containing both open-and closed-canopy ponds (excluding the cluster that contained only a single pond, Sullivan 3). Genetic analysis confirms that dispersal and gene flow likely occurs across nearby ponds and even across canopy types (Figures 2-3). Further, I found no evidence that dispersal or gene flow was limited by environmental differences among ponds. Although genetic divergence increased with geographic distance, it was not associated with canopy type (Figure 4; Table S2).
As a result, there is discordance among genetic and morphological variation among larval wood frog populations. This discordance is strongest during the selective period, when larvae are still small enough to be vulnerable to many aquatic predators (day 18), but dissipates when larvae are larger (day 37), with morphological variation becoming more closely associated with genetic divergence instead of environmental differences (Table 2). Taken together, the genetic and morphological results support the hypothesis that for some traits, selection is sufficiently strong within local environments so as to counteract gene flow that occurs.

| Variation across traits
When considering individual traits, there was variation in the extent to which each trait diverged across pond types or with gene flow. There are a number of differences among the traits that may explain the varying patterns of phenotypic divergence, including differences in selection on each trait, the type of trait, the degree of plasticity, the amount of heritability or heritable plasticity, and correlations among traits. Interestingly, the results fit the a priori expectations derived from previous research on heritabilities of and selection on each trait (Relyea, 2005;Van Buskirk & Relyea, 1998 Relyea, 1998;Relyea, 2005) showed variation due to environment (Tables 2-3; Table S3), whereas traits that have been reported to have heritable plasticity, including tail depth (Relyea, 2005), or heritable (Relyea, 2005) (Tables 2-3; Table S3). While the effects of trait differences were not explicitly assessed in this study, the results presented here in combination with results from previous studies (e.g., Relyea, 2001aRelyea, , 2005Van Buskirk & Relyea, 1998)  For species in which plasticity of some traits is under selection, such as the wood frog (Relyea, 2005), these processes may be even more difficult to disentangle. Gene flow may even facilitate the evolution of heritable plasticity by providing variation upon which selection can act (Crispo, 2008;Rasanen & Hendry, 2008). Further theoretical models and empirical examples are needed to fully understand how different phenotypic traits respond to gene flow and selection.

| Temporal patterns
Although both selection and gene flow were associated with some of the variation across traits during both time periods, these patterns were not always constant over time. At day 18, there was significant variation in morphological differences associated with environmental differences among ponds; however, this difference dissipated by day 37 (Tables S3-S4). In comparison, there was significant variation in phenotype associated with gene flow during both time periods. The only trait that continued to be associated with canopy across both time periods was body length, while the associations with both tail length and body depth dissipated (Tables   S3-S4). What these results suggest is that the effects of both selection and gene flow may vary through larval development, demonstrating the importance of measuring phenotypic divergence over multiple time points and at times that are relevant to the processes | 2511 ZELLMER being evaluated. Many studies assessing the effects of selection and gene flow on local adaptation have primarily focused on measuring traits at a single time point (e.g., Smith et al., 1997;Storfer, Cross, Rush, & Caruso, 1999;Nosil & Crespi, 2004). By incorporating temporal samples into these types of studies, we will have a better opportunity to understand the significance of phenotypic variation in the face of gene flow and may be able to uncover the occurrence of additional processes that otherwise would have been missed. For example, for some traits, the effects of selection may be compensated for over time (e.g., body depth), whereas for others (e.g., tail length) the effects of selection may be longer lasting. Such longer lasting effects may indicate potential for carry-over effects, with environmental conditions in the larval stage impacting traits in adults (Denver & Maher, 2010).

| Spatial patterns
Interestingly, in addition to differences across canopy types, there was also some unexpected evidence that morphological variation is associated with geographic genetic structure across the four geographic regions identified by the clustering analysis. Specifically, differences in overall morphology (Tables 2-3)  to different predators (Relyea, 2001a). In addition to variation in F I G U R E 3 Individual genetic cluster assignments from STRUCTURE analysis. Proportion assigned to each cluster for each individual using pond locations as a prior and assuming K = 4. Pond names above and canopy classification, open (O) and closed (C), below  correlated with the open-and closed-canopy gradient, such as resource availability, dissolved oxygen, hydroperiod, warmth (Werner & Glennemeier, 1999), or competition (Werner et al., 2007). Further research should investigate variation in selective pressures due to environmental differences that may occur at these larger spatial scales. The variation in morphology across space illustrates the importance in taking into account both environmental variation and gene flow, as differences among open-and closed-canopy ponds could be obscured by this additional variation.

| Alternative hypotheses
Two alternative hypotheses could explain the observed phenotypic differences among canopy types, including early environmental cues or maternal effects; however, there is little evidence to support either of these hypotheses. First, the phenotypic differences could be due to exposure to cues (e.g., predator chemical cues) in the ponds during the approximately 24 hours before eggs were collected. However, recent research suggests that for larval wood frogs, such cues must be associated with actual costs (e.g., chemical cues from depredated conspecifics) in order for predator-related morphologies and behaviors to be induced (Ferrari & Chivers, 2009 Rasanen & Hendry, 2008 and references therein). There is increasing evidence for gene flow among populations experiencing opposing selection pressures (Crispo, Bentzen, Reznick, Kinnison, & Hendry, 2006;Crispo & Chapman, 2008;Richter-Boix, Teplitsky, Rogell, & Laurila, 2010). Theory predicts this pattern when increasing immigration of maladapted individuals into a population increases the strength of selection within populations due to a "migration load" (Bolnick & Nosil, 2007). As a result, there may be no net change in trait frequencies across time despite immigration (Bolnick & Nosil, 2007). This mechanism has been proposed to explain trait means within isolated and connected Timema walking-stick populations (Bolnick & Nosil, 2007) and could additionally explain why there is little effect of gene flow on divergence of these wood frog populations across open-and closed-canopy ponds. Future research assessing differences in selection differentials and fitness within more isolated and connected populations will be necessary to determine if this mechanism is responsible for the observed patterns of divergence.

| CONCLUSION
Previous research on the effects of gene flow on local adaptation has provided mixed results, with support for gene flow as a constraining force (Bridle & Vines, 2007;Kawecki & Ebert, 2004;Lenormand, 2002;Slatkin, 1985), a facilitating force (Bridle & Vines, 2007), or alternatively having little influence on divergence of populations (Emelianov et al., 2004;Niemiller et al., 2008;Smith et al., 1997). The association between morphological variation in larval wood frog populations and environment within genetic clusters suggests that local adaptation may be occurring in the face of gene flow. This study adds to the mounting evidence that even over small spatial scales selection may be strong enough to overpower the homogenizing effects of gene flow (Emelianov et al., 2004;Niemiller et al., 2008;Smith et al., 1997 Genetic clusters shown in colors based on best-supported STRUCTURE run, K = 4, as shown in Figure 3. PC1 significantly differed across genetic clusters at both time periods (ANOVA: p = .007, p = .002, respectively). PC3 and PC4 showed a trend for variation across canopy (ANOVA: p = .077, p = .024, respectively) at day 18 only