The influence of breeding phenology on the genetic structure of four pond‐breeding salamanders

Abstract Understanding metapopulation dynamics requires knowledge about local population dynamics and movement in both space and time. Most genetic metapopulation studies use one or two study species across the same landscape to infer population dynamics; however, using multiple co‐occurring species allows for testing of hypotheses related to different life history strategies. We used genetic data to study dispersal, as measured by gene flow, in three ambystomatid salamanders (Ambystoma annulatum, A. maculatum, and A. opacum) and the Central Newt (Notophthalmus viridescens louisianensis) on the same landscape in Missouri, USA. While all four salamander species are forest dependent organisms that require fishless ponds to reproduce, they differ in breeding phenology and spatial distribution on the landscape. We use these differences in life history and distribution to address the following questions: (1) Are there species‐level differences in the observed patterns of genetic diversity and genetic structure? and (2) Is dispersal influenced by landscape resistance? We detected two genetic clusters in A. annulatum and A. opacum on our landscape; both species breed in the fall and larvae overwinter in ponds. In contrast, no structure was evident in A. maculatum and N. v. louisianensis, species that breed during the spring. Tests for isolation by distance were significant for the three ambystomatids but not for N. v. louisianensis. Landscape resistance also contributed to genetic differentiation for all four species. Our results suggest species‐level differences in dispersal ability and breeding phenology are driving observed patterns of genetic differentiation. From an evolutionary standpoint, the observed differences in dispersal distances and genetic structure between fall breeding and spring breeding species may be a result of the trade‐off between larval period length and size at metamorphosis which in turn may influence the long‐term viability of the metapopulation. Thus, it is important to consider life history differences among closely related and ecologically similar species when making management decisions.


| INTRODUCTION
A primary goal in landscape genetic studies is to understand the degree to which dispersal, as measured by gene flow, is facilitated or impeded by environmental factors (Manel & Holderegger, 2013;Wagner & Fortin, 2013). Dispersal, the permanent movement of an individual away from its natal location, is influenced by both intrinsic (e.g., morphology, physiology, behavior, life history) and extrinsic factors (e.g., density dependence, habitat quality), and is the primary mechanism for maintaining gene flow among populations (Clobert, Le Galliard, Cote, Meylan, & Massot, 2009;Einum, Sundt-Hansen, & Nislow, 2006). Successful dispersal is sometimes aided by morphological adaptations, such as the winged versus nonwinged morphologies of crickets (Simmons & Thomas, 2004), or behavioral adaptations, such as the propensity of cane toads to move faster and make more directed movements at the range front (Phillips, Brown, Travis, & Shine, 2008). Although these mechanisms are important drivers of dispersal ability, successful dispersal also relies on the timing and duration of dispersal events.
Life history traits, such as timing of breeding, oviposition strategies, growth rates, and fecundity, can all affect a species' ability to survive and successfully disperse to new habitat patches. Reduced dispersal between suitable habitat patches in turn decreases gene flow among habitat patches, which results in greater genetic differentiation. Across taxa, observed differences in genetic differentiation among closely related, co-occurring species can often be attributed to differences in life history traits (Dawson, Louie, Barlow, Jacobs, & Swift, 2002;Kierepka, Anderson, Swihart, & Rhodes, 2016;Whiteley et al., 2004).
In fragmented landscapes, habitat specialist species may exhibit higher degrees of genetic isolation than generalist species as is the case with eastern chipmunks (Tamias striatus) and white-footed mice (Peromyscus leucopus) in Indiana . Whiteley et al. (2004) found evidence that differences in population size and spawning habitat specificity can lead to higher levels of philopatry in bull trout (Salvelinus confluentus) which results in higher levels of genetic differentiation than the generalist mountain whitefish (Prosopium williamsoni) within the same river network. Similarly, differences in fecundity, larval period length, and population size were hypothesized to drive differential phylogeographic structure in two sympatric marine gobies (Dawson et al., 2002). Therefore, when estimating gene flow for species with complex life cycles (i.e., species that undergo an abrupt change in ontogeny, physiology, and behavior; Wilbur, 1980), it is important to consider the contribution of differences in life history traits.
Organisms with complex life cycles, such as pond-breeding amphibians, parasites, aquatic invertebrates, and butterflies, require breeding habitat that is present for the duration of their larval period to successfully metamorphose without catastrophic reproductive failure (Taylor, Scott, & Gibbons, 2006;Wilbur, 1980). Without successful metamorphosis, there is no recruitment of individuals and thus no contribution of genes to the breeding population, dispersal to new breeding populations, or colonization of new populations (Petranka, 2007;Semlitsch, Scott, Pechmann, & Gibbons, 1996;Werner, Relyea, Yurewicz, Skelly, & Davis, 2009). Furthermore, the timing of metamorphosis may greatly influence the ability of an individual to disperse.
For organisms such as amphibians that are prone to desiccation, individuals that metamorphose later in the summer may encounter a more inhospitable habitat matrix than those that metamorphose earlier in the year when moisture is higher due to spring rains. Additionally, animals that metamorphose at smaller sizes have a higher surface area to volume ratio and are more prone to desiccation which can further reduce the dispersal ability of the organism, especially when encountering an inhospitable landscape (Grover & Ross, 2000;. In this study, we investigated patterns of genetic structure and diversity for four salamander species that co-occur on the same landscape in Missouri. The ringed salamander (Ambystoma annulatum;  (Petranka, 1998). Ambystoma annulatum and A. opacum are fall breeding species whose larvae overwinter in the ponds and metamorphose in April-June (Hocking et al., 2008;Semlitsch, Anderson, Osbourn, & Ousterhout, 2014). Both species breed in late August-November and female A. annulatum oviposit eggs on submerged substrates, whereas A. opacum oviposit on the margin of the wetlands in shallow depressions and eggs hatch upon inundation (Hocking et al., 2008). Adult A. maculatum make an annual breeding migration in February-March, females oviposit on submerged substrates, and larvae metamorphose in late June-August with some metamorphosing as late as October (Hocking et al., 2008;Semlitsch & Anderson, 2016). Adult N. v. louisianensis migrate to ponds in February-March, females lay eggs singly over a period of days to weeks, adults remain present in ponds throughout the summer, and juveniles metamorphose in August-October (Gill, 1978;Hocking et al., 2008). In all four species, juveniles emigrate from ponds during or immediately after pulses of rain (Pittman, Osbourn, & Semlitsch, 2014) to forested habitats where they are highly fossorial, residing in small mammal burrows (Petranka, 1998;Semlitsch, 1981).
Although these four species are ecologically similar in that they are all forest dependent and breed in fishless ponds, differential interactions with landscape features can exist among co-occurring taxa on the same landscape (Goldberg & Waits, 2010;Mims et al., 2015;Peterman et al., 2015;Whiteley et al., 2014). Mims et al. (2015) found that water dependency was the primary driver of genetic diversity of three anuran species in the Madrean Sky Islands of Arizona. The anuran species with longer larval periods were dependent on more permanent sources of water for breeding and larval survival, and this led to decreased metapopulation connectivity as these patches were more limited (Mims et al., 2015). Similarly, comparative studies of two ambystomatid salamanders by Peterman et al. (2015) and Whiteley et al. (2014)  Our study investigates whether the genetic structure and diversity of four co-occurring salamander species is influenced by breeding phenology. We test this by extending the work of Peterman et al. (2015) by collecting additional genetic data for A. annulatum and A. maculatum as well as adding genetic data for a replicate fall breeding (A. opacum) and spring breeding (N. v. louisianensis) salamander species. We first assessed the patterns of genetic diversity and genetic structure for all four species across the same study landscape. We then tested the effect of the interpond landscape matrix on genetic differentiation using a resistance modeling approach.   Figure 2). Tissue samples were collected from all four species between spring 2013 and summer 2014 (Table 1).

| Study area and sample collection
We collected all samples during the same breeding season to minimize the among-year variation in breeding effort as most female ambystomatids do not breed annually (Semlitsch et al., 1996;Titus, Madison, & Green, 2014). When ponds were located within 100 m of each other, we pooled all individuals from those ponds into one group for downstream analyses. We chose the 100-m cutoff because 95% of adult salamanders are assumed to use the terrestrial landscape within 300 m from the edge of their breeding pond (Rittenhouse & Semlitsch, 2007) and ambystomatid adults have been observed utilizing a different breeding pond within 100 m when an annual breeding pond is dry (Gamble, McGarigal, & Compton, 2007;Trenham, Koenig, & Shaffer, 2001).

| Genetic analyses
We extracted DNA from tissue samples using Instagene (Bio-Rad, Hercules, CA, USA) following the protocol outlined in Peterman,

Connette, Spatola, Eggert, and Semlitsch (2012). We genotyped
A. annulatum at 19 microsatellite loci (  (Table   S3; , and N. v. louisianensis at nine microsatellite loci (Table S4;   USA). Samples that amplified at <80% of loci were removed from the dataset. We tested for the presence of full siblings with COLONY v2.0.5.9 (Jones & Wang, 2010). We set male and female mating to polygamous without inbreeding and ran the analysis as a long run with full likelihood, high precision, and no sibship prior. If pairs of samples from the same population were identified as having a >95% posterior probability likelihood of being related at the level of parent-offspring or full-sibling, we haphazardly selected one of the samples to keep for downstream analysis.
We tested for deviations from expected heterozygosity values under Hardy-Weinberg Equilibrium (HWE) and linkage disequlibrium among pairs of loci with Genepop on the Web (Raymond & Rousset, 1995;Rousset, 2008). We conducted both tests using 1,000 dememorization steps and 100 batches with 1,000 interations per batch and assessed the significance of our results following a Bonferroni correction for the number of comparisons (Rice, 1989). We used the "PopGenReport" package (Adamack & Gruber, 2014) in R (R Core Team 2016) to test for the presence of null alleles.
We did not correct for sample sizes across species because the metrics we used for genetic diversity (rarefied allelic richness) and differentiation (F' ST ) account for differences in sample sizes. We calculated rarefied allelic richness in the program HP-Rare (Kalinowski, 2005).
Then, we used GenoDive v2.0 (Meirmans & Van Tienderen, 2004) to calculate observed and expected heterozygosity, the standardized fixation index F' ST , and inbreeding coefficients (F IS ). We tested for isolation by distance (IBD) using the "vegan" package (Oksanen et al., 2016) for R (R Core Team, 2016) and tested for significant differences in the slopes between species with an ANCOVA in R. We assessed population genetic clustering using Bayesian assignment methods implemented in STRUCTURE v2.3.4 (Pritchard, Stephens, & Donnelly, 2000) and BAPS v6.0 (Corander & Marttinen, 2006). In STRUCTURE, we used 100,000 burn-in steps followed by 500,000 MCMC iterations for K = 2-15 under an admixture model with correlated allele frequencies and no location prior. To evaluate the STRUCTURE results, we calculated the rate of change in the log likelihood between K values (∆K; Evanno, Regnaut, & Goudet, 2005) in Structure Harvester (Earl & vonHoldt, 2012). In BAPS, we determined the most likely number of genetic partitions per species with a two-step process. First, we used the spatial clustering of individuals option to assign each sample to its most likely genetic partition. Then, we then ran an admixture analysis to refine our results (Corander & Marttinen, 2006). The admixture approach uses spatial information, Voroni tessellation, and Markov Random fields to determine the maximum number of population clusters (K). As in STRUCTURE, we tested for 2-15 potential clusters using ten replicates for each potential cluster number within each species. In both STRUCTURE and BAPS, we tested for hierarchical substructure within each putative cluster in a separate analysis for K = 2-15 with the same criteria used in the initial analyses.

| Landscape resistance analyses
We generated landscape resistance surfaces using ArcGIS v10.3 (ESRI, Redlands, CA, USA) to test our hypothesis that the species differ in their response to the landscape. We obtained our land cover data (30 × 30 m resolution) and digital elevation model (DEM; 90 × 90 m resolution) from the U.S. Geological Survey National Map server (USGS, http:// view.nationalmap.gov). All other resistance surfaces were derived from 90-m resolution USGS DEM layer in ArcGIS (ESRI) following the methods outlined in . Because the resolution of land cover data was 30 m, we resampled this surface to 90-m resolution to match our DEM-derived resistance surfaces. Resistance surfaces were defined as follows: eastness (sine of aspect; values range from 1 = east to −1 = west), northness (cosine of aspect; values range from 1 = north to −1 = south), streams (binary), percent slope, topographic position index (TPI; Jenness, 2006), topographic wetness index (TWI; Theobold, 2007), and distance from ravines. We used the slope position classification for our TPI resistance surface where the landscape was classified as follows: 1-hilltop, 2-upper slope, 3-midslope, 4-flat, 5-lower slope, 6-valley bottom using a 270 × 270 m sliding window. Topographic wetness index is a measure that is used to estimate the influence of topography on hydrological processes. None of our resistance surfaces were strongly correlated with each other as Pearson's correlation coefficients never exceeded r = .70.
We used the "ResistanceGA" package (Peterman, 2014) in R to assess the effects of distance and landscape resistance on pairwise genetic differentiation. ResistanceGA uses a genetic algorithm (GA; Scrucca, 2013) to adaptively optimize resistance surfaces through a series of transformations (continuous resistance surfaces) or by assignment of resistance values (categorical resistance surfaces). At each iteration, the relative support for a landscape resistance surface was assessed using linear mixed-effects models fit with "lme4"  Peterman, 2014). For species where more than one model emerged as describing genetic structure better than distance alone, we created a composite resistance surface for each species using a combination of all resistance surfaces that performed better than distance alone. We optimized the composite resistance surfaces using the same methods we employed for single resistance surfaces and compared the performance of the composite resistance surface and single resistance surfaces using AIC C .  (Table 2 and Table S7), and −0.050 for N. v. louisianensis (Table 2 and

| Landscape analyses
Tests for IBD were significant in all three ambystomatid species but not in N. v. louisianensis (Figure 4). There were significant differences in the relationship between genetic distance and geographic distance among species (F 7, 536 = 14.33, p < .001). The slopes of genetic distance and geographic distance for both A. annulatum (β = .021, 95% CI = 0.016-0.260) and A. opacum (β = .020, 95% CI = 0.008-0.032) were not significantly different from each other but both species have significantly greater slopes than A. maculatum (β = .004, 95% CI = 0.003-0.005) and N. v. louisianensis (β = 0.006, 95% CI = −0.006 to 0.018). In addition, we observed strong support in A. annulatum and marginal support in the other three species for isolation by resistance being a better predictor of genetic structure than distance alone (  (Meirmans, 2006). distance from ravines also being strongly supported (R 2 m = .54; Table 3). Distance from ravine was optimized such that areas closer to the ravines had high resistance, and resistance decreased as distance from the ravine increased (Fig. S2). The composite resistance surface, which combined TPI and distance from ravines, was our top model in which hilltops and upper slopes have lower resistance than lower slopes and valley bottoms (R 2 m = .63). Ambystoma opacum genetic differentiation was best predicted by the eastness resistance surface (R 2 m = .30); moreover, TWI, distance from ravine, and percent slope resistance surfaces all had ΔAIC C < 2.00 (Table 3 and Fig. S3). Using ΔAIC C, the composite resistance surface was not strongly supported for this species although it described approximately 1.5 times the variance of the eastness surface alone (R 2 m = .50). For A. maculatum, our best supported model was northness (R 2 m = .43) with eastness (R 2 m = .42) and percent slope (R 2 m = .46) both having ΔAIC C < 2.00 (Table 3; Fig. S4). Our composite surface, which combined northness, eastness, slope, streams, TWI, and distance from ravines, was well supported (ΔAIC C = 0.40; R 2 m = .50), although it was not identified as the top model (Table 3). For

| DISCUSSION
We observed different patterns of genetic diversity and structure among four salamander species that co-occur on the landscape. A. opacum, and A. maculatum) and these significant relationships were observed over a smaller spatial scale for A. annulatum and A. opacum than for A. maculatum (Figure 4). This suggests that A. annulatum and A. opacum are more dispersal limited than A. maculatum on the same landscape; a phenomenon supported by previous studies that investigated the genetic and demographic dispersal of our study species.
In a 7-year metapopulation study of A. opacum in Massachusetts, Gamble et al. (2007) found that the maximum demographic dispersal distance of juvenile A. opacum was 1,300 m. The average genetic dispersal distance for A. annulatum has been estimated to be 1,693 m which is significantly less than the estimated 2,050 m genetic dispersal distance for A. maculatum . Observed demographic dispersal distances are lacking for both A. annulatum and A. maculatum; however, adult A. maculatum have been observed moving 756 m during postbreeding emigrations (Madison, 1997) and single night total distance movements of 20-50 m for A. annulatum (Osbourn, 2012) and 53.44 m have been documented for A. maculatum (Pittman & Semlitsch, 2013). Furthermore, laboratorybased movement assays of A. annulatum and A. maculatum found that A. annulatum have a greater maximum movement distance but smaller median dispersal distance than A. maculatum suggesting that A. annulatum are capable of moving farther but do so less often than A. maculatum (B. Ousterhout, unpublished data). Similarly, the lack of an IBD effect in N. v. louisianensis is not surprising as this species has been found to be capable of dispersing at least 3 km (Gill, 1978(Gill, , 1979. Thus, the lack of spatial genetic structure and the lower degree of genetic differentiation on our study landscape for A. maculatum and N. v. louisianensis are likely a result of the ability and or propensity of these organisms to disperse over greater distances than either A. annulatum or A. opacum.
Across all species, we found support for landscape resistance describing genetic differentiation better than distance alone. With A. annulatum, we observed the strongest support for our composite resistance surface that combined TPI and distance from ravines. For this, surface ridges, flat areas, and upper slopes had lower resistance values than mid-slopes, lower slopes, and valley bottoms. As ridges and higher slopes in the Missouri Ozarks tend to be warmer and drier than lower slopes and valley bottoms , the results from both of these surfaces suggest that individuals experience lower resistance in these areas, which are typically perceived as suboptimal for amphibians that are highly susceptible to water loss (Spotila & Berman, 1976). Although these results seem counterintuitive, previous work in terrestrial salamanders suggests that individuals will move more quickly and directly through unfavorable areas in which they are physiologically stressed (Peterman, Connette, Semlitsch, & Eggert, 2014;Semlitsch et al., 2012). However, this observation could also be a consequence of pond placement as 19.57% of ponds are constructed on ridgetops and upper slopes and 57.30% of ponds are constructed on flat areas located on the ridgetops of FLW. As all species used in our study are dependent on ponds for breeding and the ponds act as stepping stones for dispersal, the fact that A. annulatum show lower resistance to movement on ridge tops could be an artifact of the pond configuration on our landscape. Additionally, A. annulatum have been observed moving through old field and pasture habitats toward breeding ponds (Briggler, Johnson, & Rambo, 2004) despite this habitat type leading to decreased survival in many species of Ambystoma likely due to increased desiccation risk and predator abundance (Rittenhouse & Semlitsch, 2006;Rothermel, 2004;Rothermel & Semlitsch, 2002) and higher resistance for gene flow than forested habitats (Crawford, Peterman, Kuhns, & Eggert, 2016;Greenwald et al., 2009).
For A. opacum and N. v. louisianensis, eastness emerged as the top model and northness emerged as the top model for A. maculatum.
These surfaces approximate the temperature and soil moisture of a landscape as north-and east-facing aspects are cooler and moister than south-and west-facing aspects. Given that amphibians are prone to desiccation (Peterman, Locke, & Semlitsch, 2013;Spotila & Berman, 1976), a lower resistance of north and east facing aspects would indicate that soil moisture may be driving the increased genetic connectivity for individuals moving through these areas. For A. opacum, resistance was lower on western slopes, which are typically warmer and drier, than on slopes with easterly aspects. Similarly, A. maculatum resistance was higher on southerly aspects than on northerly aspects. Bolded values indicate models with ΔAIC C < 2.0 in either the single surface or multiple surface optimizations. Composite = a combined resistance surface for all surfaces with ΔAIC C < 4.0; Transformation = best performing transformation of continuous resistance values selected by ResistanceGA; shape = optimal value for the shape parameter for the transformation; max = maximum value for the transformation of resistance values; R 2 m = marginal R 2 value; R 2 c = conditional R 2 value; TPI = topographic position index; TWI = topographic wetness index; K = number of parameters in model. that of Peterman et al. (2015) in that our A. maculatum samples were collected from a larger extent and we used TPI as a discrete variable with a 90 × 90 m cell size instead of using TPI as a continuous surface with a 30 × 30 m cell size. Although inference in "ResistanceGA" does not substantially change based on grid cell size (Peterman, 2014), our resistance surfaces likely encompass more heterogeneity in landscape features which can lead to different landscape genetic inferences (Short Bull et al., 2011).
In addition to landscape resistance predicting genetic differentiation, a likely explanation for the observed patterns of genetic diversity in our study is the effect of breeding phenology. We observed higher degrees of genetic differentiation and the presence of genetic clusters for the fall breeding species, A. annulatum and A. opacum, and a lower degree of genetic differentiation and a lack of genetic clustering in the spring breeding species, A. maculatum and N. v. louisianensis. In Missouri, A. annulatum and A. opacum breed and oviposit from September to early November and their larvae overwinter in the ponds before metamorphosing in late April to early June; thus, they require ponds that are continuously inundated and large enough to not freeze solid during that 6-9-month period for successful reproduction Hocking et al., 2008;Urban, 2007). Both A. maculatum and N. v. louisianensis breed and oviposit in February and March and larvae metamorphose and disperse in a large pulse between early June and late July, although metamorphosis can continue into October, meaning that these two species can utilize more ephemeral ponds on the landscape as larvae can metamorphose in as few as 3 months (Gill, 1978;Hocking et al., 2008;Semlitsch & Anderson, 2016). Although metamorphosis is often prolonged, under stressful environmental conditions, such as pond drying, amphibian larvae are able to initiate metamorphosis more quickly if individuals are larger than the threshold size for metamorphosis (Semlitsch, 1987;Semlitsch & Wilbur, 1988). Our observation of differential patterns of genetic structure between spring and fall breeding salamander species concurs with previous genetic studies of ambystomatid species that investigated patterns of genetic differentiations in species with different life histories. In Missouri, Peterman et al. (2015) found that the fall breeding A. annulatum had higher levels of genetic differentiation than the spring breeding A. maculatum on the same landscape. Similarly, Whiteley et al. (2014) observed that the fall breeding A. opacum has stronger genetic differentiation than A. maculatum among the same ponds in Massachusetts.
From an evolutionary standpoint, the observed differences in dispersal distances and genetic structure between fall breeding and spring breeding species may be a result of the trade-off between larval period length and size at metamorphosis (Petranka, 2007). Organisms with longer larval periods are able to reach a larger size at metamorphosis, a relationship that is directly linked to fitness and size at first reproduction (Scott, 1994), but they are able to utilize fewer breeding habitats in a given landscape (i.e., limited to more permanent ponds). Alternatively, the limited availability of suitable breeding habitat for fall breeding species compared to spring breeding species may have led to higher rates of philopatry in fall breeding species, as suggested by Peterman et al. (2015), because returning to successful breeding habitat may convey a fitness advantage, which may lead to decreased metapopulation connectivity (Petranka, 2007). This decreased connectivity leads to decreased gene flow in the metapopulation, as observed in our fall breeding species, which in turn lowers the potential for demographic rescue (Greenwald, 2010). The potential for demographic rescue is especially important for species, such as those used in this study, that have limited dispersal ability and boom or bust population cycles as the long-term persistence and health of the species is reliant on migration and re-colonization from other populations on the landscape in the event of local population crashes or bottlenecks (Greenwald, 2010).
Had our study only used one caudate species as a surrogate for all of the others, our inference would have differed as we would not have had the ability resolve the influence of breeding phenology on genetic differentiation in our system. Similar comparative landscape genetic studies support this idea that genetic inference can vary substantially, even among closely related and widely distributed species, when there are subtle differences in life history (Engler, Balkenhol, Filz, Habel, & Rödder, 2014;. As such, we urge stake holders to make decisions using knowledge of multiple species on the landscape even if decisions are targeted toward a single taxon.

ACKNOWLEDGMENTS
We are grateful K. Lohraff's logistical support and information about ponds. We would like to thank J. Heemeyer and numerous undergraduates for assistance with collecting field data. This research was approved by the University of Missouri Animal Care and Use Committee