Landscape genetics reveal broad and fine‐scale population structure due to landscape features and climate history in the northern leopard frog (Rana pipiens) in North Dakota

Abstract Prehistoric climate and landscape features play large roles structuring wildlife populations. The amphibians of the northern Great Plains of North America present an opportunity to investigate how these factors affect colonization, migration, and current population genetic structure. This study used 11 microsatellite loci to genotype 1,230 northern leopard frogs (Rana pipiens) from 41 wetlands (30 samples/wetland) across North Dakota. Genetic structure of the sampled frogs was evaluated using Bayesian and multivariate clustering methods. All analyses produced concordant results, identifying a major east–west split between two R. pipiens population clusters separated by the Missouri River. Substructuring within the two major identified population clusters was also found. Spatial principal component analysis (sPCA) and variance partitioning analysis identified distance, river basins, and the Missouri River as the most important landscape factors differentiating R. pipiens populations across the state. Bayesian reconstruction of coalescence times suggested the major east–west split occurred ~13–18 kya during a period of glacial retreat in the northern Great Plains and substructuring largely occurred ~5–11 kya during a period of extreme drought cycles. A range‐wide species distribution model (SDM) for R. pipiens was developed and applied to prehistoric climate conditions during the Last Glacial Maximum (21 kya) and the mid‐Holocene (6 kya) from the CCSM4 climate model to identify potential refugia. The SDM indicated potential refugia existed in South Dakota or further south in Nebraska. The ancestral populations of R. pipiens in North Dakota may have inhabited these refugia, but more sampling outside the state is needed to reconstruct the route of colonization. Using microsatellite genotype data, this study determined that colonization from glacial refugia, drought dynamics in the northern Great Plains, and major rivers acting as barriers to gene flow were the defining forces shaping the regional population structure of R. pipiens in North Dakota.

Amphibians of the northern Great Plains provide an opportunity to evaluate prehistoric climate signatures on genetic variation. This region was partially glaciated and has been characterized by glacial retreat (Mickelson et al., 1983) followed by cycles of drought and wet periods throughout the Holocene (~11 kya -present; Valero-Garcés et al., 1997;Xia, Haskell, Engstrom, & Ito, 1997). While the region southwest of the Missouri River was not glaciated, the region to the north and east was glaciated in the late Pleistocene until a rapid glacial retreat (~13 kya; Mickelson et al., 1983). This led to northward range expansions by a wide variety of species from southern refugia into newly deglaciated habitats (Masta, Laurent, & Routman, 2003;Wisely, Statham, & Fleischer, 2008;Yansa, 2006). The region remained cool and wet in the early Holocene following the glacial retreat (~11 kya-9 kya) before going through a prolonged extreme drought period (~9 kya-6 kya; Valero-Garcés et al., 1997;Xia et al., 1997). Throughout the mid and late Holocene to the present, the climate of the northern Great Plains has been characterized by milder oscillations in precipitation, ending in a current wet period following drought during the Little Ice Age and the Medieval Warm Period (950-750 BP;Fritz, Engstrom, & Haskell, 1994;Xia et al., 1997).
The modern-day northern Great Plains region encompasses a number of distinct ecoregions. The semi-arid northwestern Great Plains is located to the south and west of the Missouri River and includes badlands, steppes, and shortgrass prairie with relatively low wetlands densities (Bryce et al., 1998;Euliss & Mushet, 2004). To the east of the Missouri River, there are the glaciated plains, which include the Prairie Pothole Region (PPR), an area with millions of depressional wetlands embedded primarily in shortgrass prairie that stretches from Saskatchewan to Nebraska and Iowa (Bryce et al., 1998;Tiner, 2003). Farther east, the shortgrass prairie transitions into tallgrass prairie and prairie pothole wetlands become less abundant on the Lake Agassiz Plain, an area formed after the draining of glacial Lake Agassiz ~8 kya (Bryce et al., 1998;Barber et al., 1999).
Much of the native prairie and pothole wetlands in the Lake Agassiz Plain and the glaciated northern plains have been converted into small-grain and row-crop agriculture (~50% in North Dakota; Dahl, 1990), and are further threatened by continued agricultural expansion in the eastern part of the PPR and a potential shift to drier conditions due to climate change throughout much of the western part of the PPR (Carter Johnson et al., 2005, 2010. The northern leopard frog (Rana pipiens) is a ranid frog species that is widely distributed throughout the temperate regions of North America (Hammerson et al. 2004). Rana pipiens use a variety of habitats throughout their life cycle. Breeding ponds are used by adults in the spring and by tadpoles throughout the summer (Dole, 1971). Adults and metamorph juveniles extensively use terrestrial habitats for foraging and migration between other habitats during the summer (Dole, 1971;Pope, Fahrig, & Merriam, 2000). Rana pipiens also require well-oxygenated deep or flowing water habitats for overwintering hibernation (Cunjak, 1986;Emery, Berst, & Kodaira, 1972), which are particularly important for survival during the winter months on the northern Great Plains (Mushet 2010). Such habitats are rare during prolonged droughts, which limit the spatial distribution of northern leopard frog populations (Mushet 2010). Rana pipiens can migrate across relatively large distances among these different habitats, commonly moving 800 m, with reported ranges up to 5 km (Dole, 1971;Knutson, Herner-Thogmartin, Thogmartin, Kapfer, & Nelson, 2018).
The genetic population structure and phylogeography of R. pipiens have been extensively studied and revised since the 1970s.
Numerous described species were synonymized in the 1940s, but since the 1970s many species have been described based on morphology and/or genetics data (Hillis, 1988). Currently, R. pipiens has two distinct evolutionary lineages separated by the Mississippi River that were previously described as separate species (Cope, 1889;Hoffman & Blouin, 2004a;O'Donnell & Mock, 2012). Populations of R. pipiens east of the Mississippi River are typically more stable, have greater genetic diversity, and larger effective population sizes than populations west of the Mississippi River (Hoffman et al. 2004b;Philipsen et al. 2011, but see Mushet et al., 2013). Populations on the western edge of the range have become critically endangered and, in some cases, have already been extirpated (Corn & Fogleman, 1984;Rogers & Peacock, 2012). North Dakota lies in the transition zone between the more secure eastern populations and the imperiled western populations (NatureServe, 2017

| Sample collection, DNA extraction, and microsatellite genotyping
Forty-one populations of R. pipiens were sampled throughout North Dakota (Figure 1). Potential sampling sites were selected a priori as permanent or semipermanent wetlands as classified by the US Fish & Wildlife Service National Wetland Inventory (USFWS, 2012;Fisher 2015). The distance between a sampling site and its nearest neighbor was at least 30 km and no greater than 85 km. At each site, neonate and adult specimens were captured by actively searching the wetland perimeter. Sampled specimens were spaced throughout the wetland perimeter to reduce the likelihood of sampling-related individuals. Rana pipiens toe clips were collected from 30 individuals from each sampling site following NDSU IACUC protocol #A10047.
Toe clips were stored in individually marked vials containing 95% ethanol.

| Analysis of genetic diversity and population structure
The genotype dataset was used to derive metrics of genetic diversity and population structure in North Dakota R. pipiens. Linkage F I G U R E 1 Sampling sites (N = 41) where genetic material was collected around the state of North Dakota. Colors represent the major river basins where Rana pipiens populations clustered together. The Turtle Mountain ecoregion is also outlined as the R. pipiens population that was sampled in that region clustered separately from the other sampled populations in the Souris basin disequilibrium (LD) and deviations in Hardy-Weinberg equilibrium (HWE) were assessed using GENEPOP'007 (Rousset, 2008).
Expected heterozygosity and allelic richness for each sample site were calculated using the adegenet package (Jombart, 2008) in R v.3.4.1 (R Core Team 2017). Nei's genetic distance (G ST ) was calculated between each sampling site and was used to construct an Unweighted Pair Group Method with Arithmetic Mean (UPGMA) tree using the phangorn package (Schliep, 2011).
Population structure was examined using two Bayesian clustering algorithms, STRUCTURE 2.3 (Pritchard, Stephens, & Donnelly, 2000) and BAPS 3.2 (Corander, Sirén, & Arjas, 2008), and a multivariate method using K-means clustering following principal component analysis (PCA) of the microsatellite dataset in R. These three methods represent three statistically distinct approaches to describe population structure, allowing consistent patterns of structure to be distinguished from statistical artifacts of the clustering method. STRUCTURE and BAPS rely on different Bayesian clustering algorithms, the main difference being that BAPS has been optimized for incorporating spatial data into the clustering algorithm as a model parameter (Corander et al., 2008). The multivariate method does not rely on Bayesian inference, but instead uses dimensional reduction through PCA and a successive K-means clustering to determine numbers of populations and classify individuals into populations (Corander et al., 2008). STRUCTURE analysis consisted of an admixture model with correlated allele frequencies for each potential number of clusters (K). Each analysis consisted of 200,000 simulations after an initial burn-in of 20,000 simulations. The analysis was run for K values ranging from 1 to 41 possible clusters with 10 independent runs each. The ΔK method (Evanno, Regnaut, & Goudet, 2005) was used to identify the best-supported K value, which was determined based on the K value with the greatest ratio of change in the posterior probabilities of two sequential K values. If inconsistent results in K values were found compared to BAPS, additional nested STRUCTURE analyses were performed individually within each K group identified by the previous STRUCTURE analysis (Breton, Pinatel, Médail, Bonhomme, & Bervillé, 2008;Pereira-Lorenzo et al., 2010). These additional STRUCTURE analyses allowed for identification of potential substructuring that may have been missed.
BAPS analyses were initially performed by both "clustering of individuals" and "spatially clustering of individuals" models within the population mixture analysis. These analyses were performed using K max values ranging from 1 to 41 with 10,000 iterations per run to estimate the admixture coefficients of each sample. Ten replicates were performed for each value of K max to investigate the consistency of results for all values of K.
The adegenet package in R was used to perform multivariate analysis of population structure (Jombart, 2008). Multivariate analysis of the microsatellite data first required the reduction of dimensions using PCA, which was performed by the dudi.pca function from the ade4 R package (Chessel, Dufour, & Thioulouse, 2004). The first 100 principal components (PCs) explaining >95% of the variation in the microsatellite dataset were retained to use in K-means clustering.
Clustering for each K value was evaluated using BIC (Bayesian information criterion), and the K value with the lowest BIC value was selected for further evaluation. Sample sites were classified into populations based on which cluster the majority (>50%) of individuals from the sample site were assigned. For sample sites where no cluster represented a majority of individuals, the top two clusters within that sample site were combined into one cluster and this was applied to all samples across the dataset. This process was repeated until all sample sites contained a majority of individuals from one cluster. The identified population clusters were used as predefined groups for discriminant analysis of principal components (DAPC; Jombart, Devillard, & Balloux, 2010) using the adegenet package in R (Jombart, 2008). This method reduces the dimensionality of the genetic variation between groups using PCA and then uses the PCs produced from this analysis in a linear discriminant analysis (LDA) to create discriminant functions representing a linear combination of correlated alleles that describe the greatest amount of variation in the genetic dataset. The optim.a.score function was used to determine the optimal number of PCs to retain to best describe the population structure without overfitting the discriminant functions.
Population assignment probabilities were calculated for the optimized DAPC model, the hierarchical STRUCTURE model, and the BAPS model to assess how clearly populations were discriminated using each method.

| Landscape genetic analysis
Pairwise F ST (Nei, 1973) values were calculated for each pair of sampling sites as a measure of genetic distance between the frogs from each sampling site. Linear geographic distance between each pair of sampling sites was also calculated. A Mantel test including the linearized F ST values [F ST /(1-F ST )] and geographic distances was performed to find evidence for isolation by distance in the R. pipiens sampling sites across the state using the mantel function from the vegan package in R (Oksanen et al. 2017). A partial Mantel test was also performed to control for the variation associated with geographic distance and test the effect of the Missouri River as a barrier.
Sample sites were coded into a barrier matrix with a binomial variable representing sites on the same side (0) or opposite sides (1) of the Missouri River.
To test for global and local population spatial structure, a spatial principal component analysis (sPCA) was performed on the R. pipiens genetic dataset (Jombart, Devillard, Dufour, & Pontier, 2008). The sPCA works by maximizing the product of variance in allele frequencies and spatial autocorrelation (Moran's I) to find groups of alleles that are correlated with each other through space. An inverse square distances connection network was established to characterize the spatial relationships between each sampling site in the sPCA. The abilities of the eigenvalues produced by the sPCA to explain spatial population structure were assessed using global and local Monte Carlo tests . Three principal components with the largest positive eigenvalues were retained for further analysis.
Landscape features affecting ecological suitability and the ability of R. pipiens to move between populations were included in a redundancy analysis (RDA; Legendre & Legendre, 2012) using the three retained sPCA principal components as a measure of genetic variation of the populations across the state. Landscape factor variables included in the RDA model were land use type, the Missouri River as a barrier, and the 6-digit hydrologic unit code (HUC-6) basin within which the sampling site was located. The land use type variable was created by drawing a 15 km buffer around each sampling site in ArcMap v10.5 (ESRI). This distance represents an area where high levels of gene flow (≥1 migrant/generation) would be expected between ponds within the buffer based on the distance R. pipiens can disperse within one generation (Dole, 1971;Knutson et al., 2018).
Land use types from the National Land Cover Database (Homer et al., 2015) were reclassified into six classes: open water, urban/developed land, forest, scrubland/grassland, agriculture, and wetlands.
The area of these classes within each 15 km buffer was calculated and converted into proportions to use as the land use variable in the RDA model. The Missouri River variable was created by classifying sampling sites as being east or west of the Missouri River. Finally, the basin variable was created by classifying sampling sites based on their HUC-6 location.
A full RDA model including latitude and longitude positions of each sampling site along with all landscape features was run to determine whether any landscape variables explained a significant amount of variation in the sPCA axes. A partial RDA model conditioned on the latitude and longitude of sampling sites was run with all of the landscape feature variables to partial out variation in the sPCA axes due to isolation by distance to determine whether the selected landscape variables explained a significant amount of the remaining variation. Additional partial RDA models were run with each of the landscape factor as a single constrained variable conditioned on the other landscape factors, partitioning the variance in sPCA scores explained by each landscape factor on its own. A permutated ANOVA (PERMANOVA; Legendre, Oksanen, & Braak, 2011) was performed on each partial RDA model result with 999 permutations to determine whether the single unconditioned landscape factor in the partial RDA model explained a significant amount of the variance in the sPCA axes scores.

| Historical population coalescence and paleoclimate modeling
Population structure is strongly influenced by population demographic history, so the times of coalescence between the identified populations were estimated. DIYABC v2.0.4 (Cornuet et al., 2008),  Figure 1 population j). A stepwise mutation model with a average mutation rate of 10 −3 -10 −4 was used to describe mutation dynamics across the whole set of microsatellite markers, while individual marker mutation rates ranged from 10 −2 to 10 −5 , representing a range of mutation rates commonly seen in vertebrates, including amphibians (Bulut et al., 2009;Guillemaud, Beaumont, Ciosi, Corneut, & Estoup, 2010;Storz & Beaumont, 2002). Uniform priors were set for popula- Model validation for all models was carried out using a k-fold sampling scheme. Points were randomly assigned a K value of one through four, and four separate model optimization trials were conducted. In each trial, a subset of points assigned one of the K values was left out to validate the model after it was fit, and a subset of points with the three other K values was used as a training dataset to build the model. This process was iterated four times for each model-

| Genetic diversity and population structure
No null alleles, significant deviations from Hardy-Weinberg equilibrium, or evidence of linkage disequilibrium was observed for any of the 11 loci. Allelic richness for each locus varied from 3 to 28 with an average of 17.8 alleles per locus (Table 1)

| Landscape effects on population structure
The Mantel test indicated there was significant isolation by distance (IBD) among sample sites (r = 0.515, p < 0.001). Geographic distance was positively correlated with genetic distance ( Figure 5). The partial Mantel test corroborated the east-west divide around the Missouri River seen in the population structuring analyses. There was a significant isolation-by-barrier effect between sites on opposite sides of the Missouri River (r = 0.582, p < 0.001).
Spatial PCA found strong global structuring in the genetic variation across all sampling sites (global Monte Carlo test; r obs =0.146, p < 0.001), and there was no evidence of local structuring within sampling sites (local Monte Carlo test; r obs =0.036; p = 906). The top three sPCA axes that were retained explained 84.6% of the spatial genetic structure. The first sPCA axis explained 47.3% of the spatial genetic structure and indicated a stark split between populations on opposite sides of the Missouri River (Figure 6a).
The second sPCA axis explained 23.0% of spatial variation and showed similarities between the Souris and Sakakawea population cluster with the Little Missouri and Lower Yellowstone cluster.
The second sPCA axis also shows a stark division between the   The full RDA model explained 97.6% of the variation in the first three sPCA axes (pseudo-F = 48.7, p = 0.001). The partial RDA conditioned on the geographic coordinates of sites found a significant effect of landscape variables after removing the effects of isolation by distance (pseudo-F = 20.7, p = 0.001). The IBD effect accounted for 60.8% of the variance in the sPCA axes, 36.8% was explained by the landscape factors, and 2.4% was unexplained. Partial RDA models focusing on individual landscape variables found two variables that explained a significant amount of variation in the sPCA axes (Table 3). HUC-6 basin explained 21.0% of the variance (pseudo-F = 21.0, p = 0.001), and the Missouri River explained 4.0% of the variance (pseudo-F = 35.7, p = 0.001). None of the land use types explained a significant amount of variance.

| Population coalescence
DIYABC indicated that median coalescence times among all 10 population clusters varied from 638 to 18,100 generations for the Beaumont model (Figure 7a; Table 4) and 588 to 13,600 generations for the Cornuet-Miller model ( Figure 7b; Table 4). The median coalescence time for the southwestern and northeastern clusters separated by the Missouri River was 13,600 generations and 18,100 generations for the Cornuet-Miller and Beaumont models, respectively ( Figure 7;

| Species distribution modeling and paleoclimate projections
Sixteen models were included in the final averaged model (Table 5). For the binomial logistic regression models, there was a maximum of two models with ΔBIC<2, which were then averaged together. All of the averaged binomial logistic regression models included the same six bioclimatic variables, three related to temperature and three related to precipitation. Two precipitation F I G U R E 3 Results of hierarchical STRUCTURE population clustering. The initial STRUCTURE run split the populations of Rana pipiens in North Dakota into two main clusters separated by the Missouri River; Populations 1-12 and 14 were to the southwest, whereas all remaining populations were northeast of the Missouri. The populations west of the Missouri River broke into four clusters when analyzed separately. The populations east of the Missouri River clustered into two large groups, one in the northwestern part of the state and one in the southern and eastern parts of the state. Each of these two groups further clustered into three clusters largely separated by basin variables, precipitation during the driest month and precipitation seasonality, were present in all of the best performing models and were most important in defining the suitable range of R. pipiens.
Precipitation during the driest month was positively associated with R. pipiens occurrence, which also had a consistent negative relationship between with precipitation seasonality, suggesting F I G U R E 4 Results from DAPC population clustering analysis. (a) The first two discriminant functions explained 37.6% and 21.5% of the genetic variation in Rana pipiens from the sampled sites. Each node represents the genotype of an individual frog connected to a centroid of the cluster the frog was assigned to based on K-means clustering of the DAPC scores. (b) DAPC determined the sampled individuals were optimally clustered into ten groups, with 99.6% of individuals being assigned to one of these clusters with Q > 0.5.  (Table 5)

| Population structure, landscape effects, and prehistoric climate influences on population differentiation
There is strong evidence of population structuring of R. pipiens in North Dakota. The populations are primarily split into a western clade to the southwest of the Missouri River and an eastern clade to the north and east of the Missouri River. The eastern clade is structured into six F I G U R E 7 Coalescence trees of the ten major clusters of Rana pipiens in North Dakota based on the Beaumont (a) and Cornuet-Miller (b) models. Both models show the major split between the populations east and west of the Missouri River occurred during the late Pleistocene (18-13 kya), while most of the subdivision within these major populations occurring during the dry period of the Holocene (11-7 kya) distinct clusters, and the western clade is structured into four distinct clusters. These patterns of population structure were consistent across three structuring methods with different underlying mechanisms and demonstrate the strengths of using multiple clustering methods to analyze population genetic data. River basins and the Missouri River were the most important landscape features in determining the spatial pattern of population structuring, along with isolation by distance.
Acting as a barrier, the Missouri River has prevented gene flow between R. pipiens populations since the retreat of glaciers from North Dakota at the end of the Wisconsin glaciation. Both a partial Mantel test and the RDA variance partitioning (Table 3) (Mickelson et al., 1983), but the bioclimatic modeling indicated that this region was not suitable during the glacial maximum. However, following glacial retreat, the area west of the Missouri was likely to be the first area colonized by R. pipiens within North Dakota. Subsequently, the Missouri River has allowed the R. pipiens populations that colonized the southwest portion of North Dakota to remain genetically distinct from the R. pipiens populations that colonized the region to the east of the Missouri River.
Hoffman and Blouin (2004a)   Note. E, population located east of the Missouri River; W, population located west of the Missouri River.
following the retreat of glaciers. The colonization route these lineages Most of the population subdivision associated with river basins is most likely to have occurred during the mid-Holocene (11-6 kya), a period when the northern Great Plains region went through extreme drought cycles (Valero-Garcés et al., 1997;Xia et al., 1997).
These drought periods would have extremely restricted the breeding and overwintering aquatic habitats that R. pipiens require, possibly confining them to major riparian areas or areas with comparatively higher annual rainfall (i.e., the Turtle Mountain ecoregion; Bryce et al., 1998). In fact, river basins, which accounted for over 20% of the spatial genetic structure (Table 3), have been found to be an important factor in the genetic structuring of other frog species (Lind, Spinks, Fellers, & Shaffer, 2011;Murphy, Dezzani, Pilliod, & Sorfer, 2010 Additionally, in the prairie pothole region, overwintering habitats impose the greatest constraints on the distribution of R. pipiens (Mushet 2010). During droughts, the permanent deepwater and flowing water habitats required for overwintering are clustered in low-lying areas (Cohen et al., 2016;Van Meter & Basu, 2015). The connectivity of amphibian populations between basins is limited during dry periods when upland temporary wetlands are dry and R. pipiens are confined to overwintering sites in clusters of suitable wetlands in low-lying areas of basins (Mushet et al., 2013). Recent work has shown landscape-level decline of northern leopard frogs during severe droughts when populations are concentrated around the limited remaining winter refugia. Following the droughts, populations rapidly expanded (Mushet 2010

Temperature of Driest Quarter
ABC model indicates it occurred approximately 600 years BP, during the Little Ice Age when the northern Great Plains had a relatively arid climate (Fritz et al., 1994;Laird, Fritz, Grimm, & Mueller, 1996;Xia et al., 1997), though the confidence intervals for this divergence also include other somewhat arid periods during the late Holocene.
The dry period associated with the Little Ice Age lasted until approximately 200 years BP before shifting to the relatively wet current climate in eastern North Dakota (Fritz et al., 1994;Laird et al., 1996), a shift that likely restored connectivity and allowed the high amount of gene flow observed between the two subpopulations of R. pipiens located in the Red River basins.
Isolation by distance was the final landscape factor underlying the genetic variation of R. pipiens populations in North Dakota. Isolation by distance was strongly supported by Mantel tests ( Figure 5) and was the most important single factor explaining the spatial genetic structure of R. pipiens in North Dakota according to the RDA variance partitioning analysis (  (Hoffman, Schueler, & Blouin, 2004). Rana pipiens commonly disperse distances of 800 m from their natal ponds (Bartlet & Klaver, 2017;Knutson et al., 2018) and are capable of dispersing >5 km in some landscapes (Dole, 1971). Rana pipiens are more capable of long-distance dispersal than many other amphibians, so patterns of isolation by distance may be difficult to detect in smaller areas with high wetland density due to high levels of gene flow (Mushet et al., 2013).
Land use did not appear to influence population structure of R. pipiens at the scale analyzed in this study. North Dakota has a fairly homogenous landscape, largely comprised of grassland, pasture, and cultivated crops. Intensive agriculture and cattle grazing can have negative effects on the dispersal ability of some amphibians (Mushet, Euliss, & Stockwell, 2012;Rothermel & Semlitsch, 2002;Vos, Goedhart, Lammertsma, & Spitzen-Van der Sluijs, 2007).
However, there is evidence that R. pipiens can use cultivated crop fields as well as native grasslands if the amount of cover and level of soil moisture in agricultural fields are comparable to the native grasslands allowing frogs to avoid desiccation (Bartlet & Klaver, 2017;Pope et al., 2000). Agricultural development may play a subtler role influencing connectivity within R. pipiens metapopulations than could be detected at the spatial scale used in this study. Finer scale sampling using higher-resolution genetic markers may be able to detect local effects of land use on the biotic connectivity between prairie wetlands.

| CON CLUS ION
Rana pipiens populations in North Dakota are highly structured, typical of amphibian populations at regional geographic scales. The patterns of genetic population structure reflect the colonization Rana pipiens populations southwest of the Missouri River, though currently stable, are characterized by lower genetic diversity (Stockwell et al., 2016) and smaller effective population sizes. These populations are at the highest risk of experiencing population declines similar to those reported on the western edge of the R. pipiens range. Conservation of wetlands and riparian areas across the region will be important, both for securing at risk R. pipiens populations on the Missouri Plateau and Badlands, and for maintaining the genetic diversity found throughout the R. pipiens populations on the glaciated northern Great Plains.
Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

ACK N OWLED G M ENTS
We thank Jacob Mertes for assistance with fieldwork; Patrick

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

AUTH O R CO NTR I B UTI O N S
CAS, JDLF, KP, and DM designed the project. JDLF performed fieldwork and laboratory work. KP and CAS provided guidance on data analyses. JW and JDLF conducted data analyses. CAS supervised the study. All authors participated in writing the manuscript.