Introduction

Land-use change causing habitat loss and fragmentation is one of the greatest threats to biodiversity globally (Mantyka-Pringle et al. 2012; Haddad et al. 2015; Wilson et al. 2016) and is of particular concern in the Prairie Pothole Region (PPR) of the northern prairies of North America. Here, biotic connections which link prairie-pothole wetlands and associated dryland environments (Smith et al. 2019) are now threatened by accelerating land-use change, as agricultural conversion has removed up to 90% of wetlands in some parts of the PPR (Dahl 1990, 2014; Van Meter and Basu 2015).

Fragmentation of the remnant grasslands and wetlands in the PPR has likely impacted biotic connectedness (Wright 2010; Wimberley et al. 2018). Knowledge of biotic connectivity among habitats is critical information for conservation efforts, as connectivity can counteract the negative effects of isolation on the average fitness within populations and metapopulation extinction risk (Simberloff 1988; Saunders et al. 1991; Holyoak and Heath 2016; Mable 2019). Most exploration of landscape connectivity in the PPR has been through generic modeling (Wright 2010; Wimberley et al 2018), however, connectivity of populations can vary greatly depending on the dispersal abilities of different species (Cushman et al. 2013; Mims et al. 2015). Landscape genetic analyses can be powerful tools to empirically assess and improve models of habitat connectivity, revealing patterns of gene flow and identifying environmental factors influencing the migration of organisms (Epps et al. 2007; Segelbacher et al. 2010; Dickson et al. 2018).

Amphibian meta-populations in the PPR are ideal for studying biotic connectedness among prairie-pothole wetlands because of their need for an intact mosaic of grassland and wetland habitats (Mushet et al. 2014; Holm 2020). Because amphibians are often dispersal-limited, habitat fragmentation poses a high risk to amphibian metapopulations and communities (Lehtinen et al. 1999; Funk et al. 2005; Belasen et al. 2019). Further, amphibians can serve as useful bio-indicators of wetland condition (Bowers et al. 1995), as they occupy roles in both aquatic and terrestrial food webs (Rittschof 1975; Benoy 2005; McLean et al. 2021).

The northern leopard frog (Rana pipiens) is common throughout the PPR (Lehtinen et al. 1999; Balas et al. 2012) and relies on numerous habitat types, including semi-permanently ponded, fishless wetlands for breeding (Dole 1971; Herwig et al. 2013), grasslands and meadows for adult foraging (Rittschof 1975; Pope et al. 2000), and deep-water or flowing-water habitats for overwintering (Emery et al. 1972; Cunjack 1986; Mushet 2010). Rana pipiens is also a very vagile amphibian species with individuals often migrating > 800 m between overwintering and breeding habitats each breeding season. Individuals are also capable of long-distance dispersals that can exceed 3 km (Dole 1971; Knutson et al. 2018). Rana pipiens rely on different seasonal habitats, and annual movements through diverse matrix habitats make connectivity especially important for the persistence of populations in the PPR.

Previous genetic studies of R. pipiens in North Dakota used microsatellite markers and indicated that, statewide, populations were structured by prehistoric climate, glacial history, and river basins (Waraniak et al. 2019). We designed this study to investigate factors that may affect gene flow and connectivity within an area undergoing an exceptionally high rate of agricultural conversion.

Methods

Sample site selection and sample collection

We conducted our study in the North Dakota portions of the Lake Oahe and James River (Hydrological Unit Code HUC-6) river basins. Previous genetic analyses of R. pipiens populations in North Dakota indicated that in these two river basins the species formed a distinct genetic cluster (Waraniak et al. 2019). We selected sample sites within the Lake Oahe and James River Basins by screening National Wetland Inventory (NWI) data for permanently ponded wetlands located on public lands. Twenty-eight sample sites at least 10 km apart were selected based on land-use data from the National Land Cover Database (NLCD) in a 5 km buffer around each wetland. Sites were selected to ensure that a wide variety of land use was represented within our sample. Additionally, the study area spans two Level III ecoregions, the relatively flat Northern Glaciated Plains (hereafter referred to as the Drift Prairie) to the east and the hillier Northwestern Glaciated Plains (hereafter referred to as the Missouri Coteau) to the west, each shaped by different geological histories resulting in distinct topographies and configurations of wetlands (Bryce et al. 1996).

Sample sites were visited 1 to 2 times in 2018 and 2019 between late July and early August with the purpose of collecting toe clip tissue samples from 20 individuals from each site. The timing of the sampling in late summer was chosen to primarily target neonates that had recently undergone metamorphosis so we could be reasonably certain that the sampled individuals originated from the target wetland. All animal handling and tissue collection was permitted by the North Dakota Game and Fish Department (License numbers GNF04736156 and GNF04910652 for 2018 and 2019, respectively) and strictly following North Dakota State University Institutional Animal Care and Use Committee (IACUC) approved protocol #18,052.

DNA extraction and BestRAD library construction

DNA was extracted using the Qiagen DNeasy Blood & Tissue kit according to the manufacturer’s protocol. Genomic libraries were constructed using the BestRAD protocol as described in Ali et al. (2016). DNA was digested using SbfI restriction enzyme prior to ligation of P1 RAD adaptors that contained one of 96 unique identifier sequences for each individual within a library. DNA from all individuals within a library was then pooled and incubated with NEBNext dsDNA fragmentase for 45 min at 37 °C. DNA libraries were then concentrated using Ampure XL beads and prepared for Illumina sequencing using the NEBNext Ultra II DNA Library Kit for Illumina, selecting for an average library fragment size of ~ 400 bp. Completed libraries were shipped to the University of Minnesota Genomics Center and sequenced on an Illumina NovaSeq V1 flowcell with 150 bp paired-end reads.

Genomic data processing

Raw genomic data was processed using STACKS v.2.41 (Catchen et al. 2011, 2013). Sequence data was demultiplexed using the ‘process_radtags’ command, assigning each sequence to the correct individual and removing sequences of poor quality. Following demultiplexing, the reads were passed through the ‘clone_filter’ command in STACKS to remove polymerase chain reaction (PCR) duplicates introduced during the library construction process. Next, RAD loci were identified with both a de novo approach and by aligning the reads to the R. catesbeiana reference genome (Hammond et al. 2017) using the bwa mem alignment algorithm (Li 2013). Successfully aligned reads were then run through ‘gstacks’ in STACKS which identifies polymorphic loci. Genomic data was exported from the STACKS pipeline in vcf and genind formats using the ‘populations’ module without any filters. Genotype data filtering was performed and threshold values for quality control filters were set using the ‘radiator’ package in R (Gosselin 2020). Individuals with more than 60% missing data and heterozygosity values that fell outside of 0.02 to 0.65 were removed from the dataset. Loci that failed the following filters were removed from the dataset: (1) minor allele frequencies below 3% (minor allele count less than 12), (2) less than 4X coverage, (3) had greater than eight single nucleotide polymorphism (SNPs), (4) heterozygosity above 0.65, (5) significantly deviated from Hardy–Weinberg Equilibrium within sample sites, or (6) had a global missing data rate of greater than 30%. To avoid effects of linkage, only one SNP with the highest minor allele count was kept for each locus. Putative siblings were identified from individuals collected at the same sample site using COLONY (Jones and Wang 2010). For any groups of siblings identified, the individual with the lowest amount of missing data was retained. The final filtering step was to remove populations that had fewer than eight individuals remaining after the ‘radiator’ filtering steps and putative siblings had been removed (Willing et al. 2012).

Population structure analyses

General population genetic statistics (expected/observed heterozygosity, FIS) were calculated using the ‘hierfstat’ package in R (Goudet and Jombart 2022). Nucleotide diversity (π) for the filtered loci included in the dataset was calculated using the populations module in STACKS and included all SNPs present on the loci. Effective population size (Ne) for each sample site was estimated using the linkage disequilibrium with random mating model, heterozygote excess, and molecular coancestry methods in NeEstimator v2 using markers with > 0.05 minor allele frequency within sample sites (Do et al. 2014).

Population genetic structure was identified using several methodologies, including STRUCTURE v2.3.4 (Pritchard et al. 2000), k-means clustering and discriminant analysis of principal components (DAPC), and fineRADstructure (Malinsky et al. 2018). For the initial exploratory STRUCTURE analysis, we ran 10 replicates for each number of clusters, assuming a maximum of 22 clusters (from K = 1 to K = 22) with each replicate consisting of 40,000 burn-in steps and 25,000 iterations. A second set of more robust STRUCTURE runs was performed on the most likely values of K, using 5 replicates for each number of clusters from K = 1 to K = 9, with each replicate consisting of 500,000 burn-in steps and 1,000,000 iterations. The STRUCTURE model was used to infer lambda for three iterations, which is appropriate for this dataset which contains many rare alleles (Falush et al. 2003). The mean of the inferred lambda across these iterations was used as a constant lambda for the full STRUCTURE run. No prior location information was used for the STRUCTURE model. The optimal number of clusters was determined using mean log-likelihood, Evanno’s ∆K, evaluating the change in log likelihood of the dataset given the STRUCTURE model as the number of clusters increased (Evanno et al. 2005) and using the MedMeaK, MaxMeaK, MedMedK, and MaxMedK statistics with a threshold of 0.5 (Puechmaille 2016).

A successive k-means clustering method was used to determine the clusters to be used in the discriminant analysis of principal components. For each value of K (from K = 1 to K = 22), a k-means clustering analysis was performed and evaluated using Bayesian information criterion (BIC) and by comparing the ratio of between cluster sums of squares to the total sums of squares. The number of clusters best supported by these two k-means metrics was used to define the clusters for the DAPC. Sample sites were assigned to clusters based on which cluster a majority of the individuals within the sample site were assigned to by the k-means analysis. The optim.a.score function from the ‘adegenet’ package was used to determine the number of principal component axes to retain to avoid overfitting the discriminant functions (Jombart 2008).

A dataset containing all SNPs from the filtered loci was analyzed using fineRADstructure, a program that models co-ancestry based on linkage of SNPs within RAD loci (Malinsky et al. 2018). We ran the MCMC with 100,000 burn-in steps followed by 100,000 iterations. Assessment of convergence and visualization of results was conducted using the scripts provided with fineRADstructure in R.

Landscape genetics analyses

Three measures of pairwise genetic distances were calculated for each combination of sample sites: an empirical Bayes estimate of FST (Kitada et al. 2007, 2017), Nei’s GST (Nei 1973), and Jost’s DJ (Jost 2008). Pairwise geographic distances were calculated for each pair of sites using Great Circle distances based on the WGS84 ellipsoid. Mantel tests were performed for each set of pairwise genetic distances to test for an effect of isolation by distance. Mantel tests were bootstrapped with 999 permutations to calculate significance. Additionally, because HUC-6 level river basins were found to structure R. pipiens populations on a statewide level (Waraniak et al. 2019), maximum likelihood population effects (MLPE) models were fit to each of the genetic distance measures to see if HUC-8 level subbasins similarly structured R. pipiens populations on this smaller scale. To test for subbasin effects, the MLPE analysis compared the fit of models including a constant factor of random population effects, and two variables (1) pairwise geographic distances and (2) a subbasin variable represented by a pairwise matrix in which sample site pairs had a value of zero if located in the same HUC-8 subbasin and a value of one if located in different subbasins (Fig. S1). The fits of MLPE models were assessed using Akaike’s Information Criterion adjusted for small sample sizes (AICc).

Landscape effects on genetic diversity were assessed using linear models with measures of genetic diversity [expected heterozygosity (Hs), nucleotide diversity (π), and inbreeding coefficient (FIS)] as the dependent variables. Four landscape variables previously shown to affect genetic diversity of leopard frogs at larger scales in North Dakota (Stockwell et al. 2016) were included in the full models. These included latitude, longitude, mean annual precipitation (WorldClim), and proportion of wetland land cover in a 15 km buffer around each sample site to replicate the analysis conducted in Stockwell et al. (2016). A fifth variable based on the proportion of cultivated crop land cover in a 15 km buffer around samples was also included. All combinations of variables were fit using the glm function in R and relative variable importance was calculated using AICc weights.

The effects of landscape on genetic connectivity were assessed with circuit theory resistance distances implemented through Circuitscape in Julia (McRae 2008; Anantharaman et al. 2019). Ten landscape variables that could contribute to resistance were included in the analysis (Table 1, Fig. 1). Rasters for each of these landscape variables were clipped to the study area (50-km buffers around the sample sites that passed quality control filters within an extent between 45.5ºN and 48ºN latitude and 97.5ºW and 101.5ºW), aligned with a snap raster to ensure that the spatial dimensions of each landscape variable were identical, and coarsened to a cell size of 120 m × 120 m to reduce computational load. The resistance models for each landscape variable were first optimized as a single resistance surface using custom R scripts and the 'ResistanceGA’ package in R (Peterman 2018). Continuous surfaces were transformed using monomolecular or Ricker functions and the best performing transformation was assessed using AICc calculated from maximum likelihood population effects models with the genetic distance as the dependent variable, sample site as a random effect, and resistance distance as the sole fixed explanatory variable. The optimized resistance surfaces were then combined in combinations of up to four surfaces in an additive process where all remaining surfaces were added to single surfaces or surface combinations with ∆AICc < 2 and optimized. Circuit distances were calculated for each single and combined optimized surface and modeled against Nei’s GST and Jost’s DJ pairwise genetic distances. These resistance models were assessed using AICc calculated from MLPE models as described above.

Table 1 List of landscape variables tested in circuit distance-based landscape resistance analysis
Full size table
Fig. 1
figure 1

Maps of study area showing patterns of the 10 landscape variables used in the Circuitscape-based ResistanceGA analysis, organized by categories according to Table 1

Full size image

Results

SNP dataset

After all STACKS and ‘radiator’ filtering steps, the final dataset contained 2868 loci. Thirty-three individuals were identified as putative siblings by COLONY, so 23 of the individuals were removed to retain only one individual from each sibling group. After individuals had been removed for high levels of missing reads as well as putative sibship, six sample sites had fewer than eight individuals remaining, so these populations were removed from the dataset. The final dataset contained 315 individuals across 22 sample sites with a mean missing data rate of 23.78% (min 4.01%, max 57.81%). Patterns of missingness were consistent across sample sites except for slightly elevated missingness in the sample sites of Chase Lake NWR, Sykeston, and Wing where nearly all individuals had missingness higher than the dataset average (Table S1.)

Population genetic summary statistics and population structure

The sample sites had similar measures of observed heterozygosity (mean ± sd = 0.0748 ± 0.0073), inbreeding coefficient (mean ± sd = 0.0029 ± 0.0135), and π (mean ± sd = 0.178 ± 0.004; Table 2). Effective population size estimates ranged from 72.9 to infinite for the LD method and from 31.0 to infinite estimates based on heterozygosity excess (Table 2). No private alleles were detected within sample sites.

Table 2 Summary of basic population genetic statistics for the northern leopard frog, Rana pipiens
Full size table

Mean log-likelihood of STRUCTURE results increased gradually across numbers of K from K = 1 to K = 9 (Fig. S2A). The Evanno’s ∆K method of interpreting the STRUCTURE analysis indicated K = 2 as the best estimate for the number of genetic clusters, though this method cannot be used to assess the fit of K = 1 (Fig. S2B). MedMeaK, MaxMeaK, MedMedK, and MaxMedK all estimated the number of genetic clusters to be equal to 1. Overall, the STRUCTURE analysis suggests that all sites belong to the same genetic population.

Optim.a.score indicated that 58 principal components should be retained to prevent overfitting of clustering models for use in DAPC clustering. In the successive k-means clustering analysis on the 58 retained principal components, clusters from K = 4 to K = 6 had similar BIC values (Fig. S2C). K-means clustering results were further assessed using the ratio of between cluster sums of squares to the total sum of squares (Fig. S2D). The amount of variation explained by adding new clusters consistently increased up to K = 6. The clusters found by k-means clustering showed clear spatial patterns, though there was also evidence of admixture among the different clusters within study sites (Fig. 2).

Fig. 2
figure 2

Map of sample sites used in population structure analysis. Sample sites are shown as pie graphs with the colors representing the average proportion of assignment probability to one of the K = 6 genetic clusters as estimated by discriminant analysis of principal components (DAPC) for individuals in that sample site. The six sample sites that were dropped due to insufficient numbers of individuals that passed quality control filters are represented by black points. Background shows USGS NLCD land use categories. DAPC results for Discriminant Functions 1 and 2 with K = 6 genetic clusters are presented in a scatterplot graph below the map

Full size image

The analysis of co-ancestry using fineRADstructure indicated little evidence for strong population structure within our samples (Fig. S3). The co-ancestry model converged on 22 unique clusters, however, most of the inferred clusters among the samples were weakly supported. When only clusters that received at least a moderate level of support (appearing in > 30% of MCMC iterations) were considered, the resulting clusters were similar to the clusters determined by k-means analysis.

Landscape genetics analyses

Mantel tests showed significant effects of isolation by distance among the sample sites for two of the three genetic distances tested (Fig. 3; EBFST P = 0.001; GST P = 0.187; DJ P = 0.025). MLPE analysis for all three genetic distances showed that the model including both geographic distance and subbasins best fit the genetic data, but in all cases these models had similar fits to the model with geographic distance alone (∆AICc 2.2–3.7; Table S2).

Fig. 3
figure 3

Plots of three pairwise genetic distances against geographic distance. All genetic distances had a significant relationship with geographic distance, indicating a pattern of isolation by distance among the sample sites. Shaded area is 95% confidence interval

Full size image

Analysis of genetic diversity showed that lower, often negative, values of FIS were associated with areas with higher proportions of agriculture, but did not indicate that any of the landscape variables tested had a substantial effect on He or π in the study area (Table S3). All the models testing the effects of landscape variables on FIS with ∆AICc < 2 included a negative relationship with the proportion of agriculture in a 15 km buffer. Agriculture was the most important variable associated with FIS (variable importance = 0.85), no other landscape variables had a consistent effect (variable importance < 0.5) The best ranked model for expected heterozygosity included longitude as a factor, though this model was only slightly favored over the intercept-only model (∆AICc = 1.09). All model variables had an importance < 0.5.

The analysis of landscape resistance showed there were consistent patterns in landscape variables associated with genetic differentiation among the sample sites (Figs. 4, S4, S5). At least one of the land-use or topographic roughness variables was always included in the top models for both genetic distances (Nei’s GST and Jost’s DJ) used in the Circuitscape analysis (Table 3). Model averaged values for land-use categories showed low-disturbance dryland habitats generally had the lowest resistances, including the following land-use categories: forests, developed—low intensity/open space, scrubland, and grassland (Table 4). Agricultural land uses had moderate resistance values, with pasture having lower resistance than cultivated cropland. The areas with the highest resistance were typical barriers like open water and medium/high intensity urban areas, along with wetland core habitat. The topographic roughness components of all top models show agreement that flat terrestrial areas have the lowest resistance. Areas with 0 roughness index values occurred on bodies of water and had very high resistance values. Variables related to planting frequency of different crops also had some marginal effects as areas where one crop type (either corn, soybeans, or winter wheat) was frequently planted tended to have higher resistance (Table 3). Only winter wheat appeared in the top models based on GST, while all three crops contributed to top models based on Jost’s D. The major road variable did not appear in any of the top models, and variables related to water stress only appeared in top models with land use or topographic roughness variables and only made minor contributions to the overall resistance in these models.

Fig. 4
figure 4

Maps of the study area showing the weighted model average landscape resistances for models based on GST genetic distances (A) and DJ genetic distances (B)

Full size image
Table 3 List of the best fitting landscape resistance models that contribute at least 1% to average model for Circuitscape resistance distances based on both GST and Jost’s D genetic distances
Full size table
Table 4 Model averaged resistance values for the land use types for both sets of models based on GST and Jost’s D
Full size table

Discussion

Our landscape genetic analysis of Rana pipiens in the PPR suggests that land use and topography are the primary drivers of population connectivity in the region, although there did not appear to be strong effects of fragmentation in the frog populations surveyed. Apart from some localized spatial clustering and a general pattern of isolation by distance, leopard frogs in the James River and Lake Oahe basins did not show signs of strong genetic differentiation. There were also no signs of extensive inbreeding (all FIS values < 0.03). While current genetic patterns suggest that R. pipiens populations in this area are characterized by extensive gene flow, increasing landscape resistance to gene flow due to agricultural expansion may still pose a risk.

Land-use effects on connectivity

Landscape resistance analysis of R. pipiens genetic connectivity produced similar patterns of resistance to different land-use categories as previous studies. Open water, urban areas and wetlands provided the highest resistance to gene flow. Open-water areas and urban areas are clear barriers to movement, as both pose high risk of mortality to leopard frogs. Open-water areas are generally avoided, as the lack of vegetative cover leaves northern leopard frogs exposed to high predation risk from both fish and avian predators (Lundgren et al. 2012; Geluso et al. 2013). Likewise, urban areas are often associated with reduced abundances in R. pipiens populations (Gagné and Fahrig 2007; Johnson et al. 2011). Developed land uses can lead to increased adult mortality due to higher density of roads and traffic (Carr and Fahrig 2001; Glista et al. 2008) and higher tadpole mortality due to the presence of invasive and introduced species (Germaine and Hays 2009; Johnson et al. 2011). Additionally, higher exposure to pesticides, road salts, and other chemicals from roadway and urban runoff have negative effects on larval development and behavioral responses (Croteau et al. 2008; Denoël et al. 2010).

It appears counterintuitive that wetlands, the preferred breeding habitat of R. pipiens, would have resistance similar to that of land uses associated with high mortality risk for amphibians, but this pattern of preferred habitat having high resistance values has been observed in other landscape genetic studies on amphibians (Angelone et al. 2011; Peterman et al. 2014). While the explanations for this phenomenon are speculative, for the landscape context and scale of our study area, wetlands may be acting as settling locations for migrants. The PPR is characterized by abundant wetlands separated by drylands, so frogs would necessarily have to pass through non-wetland matrix habitat when they disperse. Migrating frogs may be likely to settle near wetlands, but move through other habitat types.

Agricultural land uses had intermediate resistance estimates, offering more resistance than low disturbance natural habitats, but less resistance than open water and urban land uses which function more as strict barriers. Experiments in the PPR suggested that cultivated crops did not necessarily impede R. pipiens movement or increase desiccation risk of R. pipiens, though these effects were strongly dependent on weather conditions and seasonality (Inczauskis 2017). Telemetry studies of R. pipiens in similar agricultural environments found that migrating frogs avoided traveling through cropland but could cover long distances through crops under certain conditions (Bartlet and Klaver 2017; Knutson et al. 2018).

In our analysis, pasture and hayfields appeared to have much lower resistance compared to cultivated crops, which most landscape genetic analyses often combine into a single category. Movement studies have recorded R. pipiens using pasture and hayfields as dispersal corridors, though these land uses come with additional risks of mortality from being killed by agricultural machinery (Joseph and Johnson 2012; Knutson et al. 2018). Overall, it appears that R. pipiens in the PPR exhibit a similar relationship to agriculture as other anurans in agricultural landscapes, i.e., cropland is not a preferred matrix habitat but is at least somewhat permeable to gene flow (Goldberg and Waits 2010; Richardson 2012; Youngquist et al. 2017). Additionally, because much of the expansion of cultivated-crop land use has been relatively recent, the resistance values may be somewhat conservative as recent isolation may not be fully reflected by genetic distances, referred to as time lag (Manel and Holderegger 2013; Epps and Keyghobadi 2015).

One of the clearest signals from our landscape resistance analysis was the importance of low-disturbance land uses in maintaining connectivity. Grassland, shrubland, and forests usually had the lowest resistance scores among the best models. Grassland is particularly important in the PPR as it comprised 25.7% of the land use in the study area while shrubland and forest land-uses combined consisted < 1% of the total land area and had more variation in resistance among models. It is well documented that R. pipiens extensively uses grassland habitat for foraging and dispersal (Dole 1965, 1971; Pope et al. 2000), and wetlands networks with large amounts of grassland matrix habitat have been associated with greater occupancy and movement in other amphibian species as well (Balas et al. 2012; Mushet et al. 2014; Mester et al. 2020).

Differences in gene flow patterns between the drift prairie and the Missouri Coteau

The Drift Prairie and Missouri Coteau are characterized by different landscape features. The Missouri Coteau is more topographically complex, with numerous moraines and glacial depressions compared to the relatively flat plains of the Drift Prairie (Bryce et al. 1996). This difference in topographic roughness is reflected in the patterns of gene flow from the Circuitscape models. In the Missouri Coteau, topographic roughness contributed 79% to the total landscape resistance in the GST model averages. As expected, topographic roughness in the Drift Prairie contributed less to resistance (67%) to the GST model averages. Topographic roughness contributed less to resistance overall in DJ models, but showed similar patterns to the GST models when comparing the two ecoregions (Table S4).

These two regions also differed in land-use patterns. The Drift Prairie was dominated by agriculture, with cultivated cropland comprising 55.3% of the total land use. By contrast, the Missouri Coteau had a more varied land use dominated by grassland (45.6%), cultivated crops (9.6%), and pasture (7.9%). Land use contributed 30% and 18% of the total landscape resistance in the Drift Prairie and Missouri Coteau, respectively. Differences in landscape context are often invoked to explain contrasting resistance models from multiple disjunct populations (Trumbo et al. 2013; Micheletti and Storfer 2015; Emel et al. 2019) or in comparisons to the results of other studies on the same or similar organisms (Goldberg and Waits 2010; Coster et al. 2015; Youngquist et al. 2017). Our findings show how landscape context can alter gene flow in different parts of a single metapopulation.

Lack of strong population structure

Despite this study covering an area (maximum distance of 244 km between sample sites) much greater than the dispersal distance of an individual frog (maximum recorded dispersal distance of 5 km; Dole 1971), there was no indication of strong population structure in our genetic data. This result was not entirely unexpected, as a previous statewide survey of genetic diversity of Rana pipiens in North Dakota indicated that the R. pipiens sampled in the James River and Lake Oahe River basins belonged to the same genetic cluster (Waraniak et al. 2019). Likewise, genetic analysis of R. pipiens on a more local scale in this region also found no evidence for strong genetic structure (Mushet et al. 2013). Genetic clustering analyses more sensitive to fine-scale structure (DAPC and fineRADstructure) reported here detected some spatially organized genetic clustering, however, most sample sites contained a mixture of the genetic clusters identified, suggesting high rates of gene flow among these genetic groups.

The dispersal behavior of R. pipiens and the recent climatic history of the eastern PPR may help explain the weak genetic population structure found in this study. Dispersal and colonization rates appear to be exceedingly high and site fidelity relatively low for northern leopard frogs in parts the Upper Midwest with high wetland densities compared to drier parts of the species’ range (Merrell 1970; Lehtinen and Galatowitsch 2001; Mushet 2010). Effective population size estimates of the sampled populations also suggest large breeding populations. Large population sizes along with extensive gene flow would limit genetic differentiation due to genetic drift and are likely explanations for the lack of strong population structure observed in this study. Additionally, this portion of the PPR has been experiencing an unprecedented wet period for the past 30 years (McKenna et al. 2017). During drought years, R. pipiens have limited suitable overwintering habitat, and this occasionally leads to sharp bottlenecks that would increase the likelihood of population differentiation due to genetic drift (Mushet 2010). However, the extended wet period likely prevented any recent bottlenecks and may have allowed increased gene flow between genetic clusters that may have been less connected during typical climate conditions.

Conservation implications for amphibians in agriculture-dominated landscapes

While land-use change and intensive agriculture have appeared to impact connectivity among R. pipiens in the PPR, leopard frogs have been able to maintain enough gene flow to avoid serious genetic consequences from habitat fragmentation. There was no detectable decline in genetic diversity for R. pipiens populations associated with higher amounts of nearby agricultural land use, and the lack of population structure also suggests that gene flow is frequent among breeding sites throughout the study area. This evidence indicates that R. pipiens have been resilient to agricultural conversion (at least from a population genetic standpoint), likely due to their high vagility and an extended period of favorable climatic conditions. However, it is possible that genetic signals may lag (Epps and Keghobadi 2015) relative to the recent rapid expansion of row crops on the landscapes.

There are also alternative possible explanations for the patterns seen in our connectivity analyses. The patterns of gene flow detected might be reflecting historical land use and habitat types correlated with modern land use, and the effects of modern land-use patterns are too recent to be detected in the genetic data. However, most conversion of grasslands to cultivated cropland in the PPR occurred > 40 years ago (Rashford et al. 2011), which is > 20 generations for R. pipiens (Hoffman et al. 2004). Additionally, patterns in gene flow appear to be relatively stable across time because both genetic distances produced highly similar best-fitting models, with GST expected to be more sensitive to population structure caused by recent changes in gene flow or genetic drift while Jost’s D is expected to reflect population structure due to colonization history (Ma et al. 2015). Finally, there could be other more important variables shaping gene flow that were not considered in our connectivity analysis. This is a shortcoming in all landscape genetic studies using landscape variables as predictors of connectivity. However, the models generated in our analyses outperformed null models of uniform resistance distances and Euclidean distances by wide margins, so at least some of the more important factors influencing gene flow were likely identified in our study or were correlated with important drivers we did not measure.

While the genetic connectivity of R. pipiens populations appears to be resilient to agricultural land-use change, our study also highlights the importance of grassland conservation for amphibians, as grasslands were most often the habitat type with the lowest resistance among the top performing models of connectivity. Suitable matrix habitat may be even more important for amphibian species with lower dispersal capabilities [e.g., wood frogs (Rana sylvatica), and tiger salamanders (Ambystoma tigrinum)] for maintaining genetic connectivity. These species exhibit population structure on smaller geographic scales than leopard frogs and may be more susceptible to fragmentation due to increased landscape resistance from agricultural land use (Spear et al. 2005; Crosby et al. 2009). More connectivity studies of amphibians in the PPR could determine if agriculture is having more severe impacts on other amphibian species and if matrix habitat requirements differ among species.

Grasslands with little elevation change are particularly vulnerable to agricultural conversion (Rashford et al. 2011) and are also the most valuable matrix habitat for R. pipiens as identified in our analysis. By including these habitat parameters when making decisions, our results can help land managers and conservationists prioritize areas for conservation and restoration to aid connectivity among amphibian populations. Conservation organizations already use connectivity models to identify valuable habitat patches, often based on species of conservation concern or species-agnostic expert opinion models (Dickson et al. 2018). Adding empirically-backed modelling, especially for underrepresented taxonomic groups with unique habitat requirements, can help improve these land-management decisions (Krosby et al. 2015; Reed et al. 2017). Fragmentation of amphibian habitat is a threat throughout the Great Plains (Heintzman & McIntyre 2021) and using information from R. pipiens connectivity models can ensure that dispersal needs of amphibians are considered in landscape management decisions.