The effect of sampling density and study area size on landscape genetics inferences for the Mississippi slimy salamander (Plethodon mississippi)

Abstract In landscape genetics, it is largely unknown how choices regarding sampling density and study area size impact inferences upon which habitat features impede vs. facilitate gene flow. While it is recommended that sampling locations be spaced no further apart than the average individual's dispersal distance, for low‐mobility species, this could lead to a challenging number of sampling locations, or an unrepresentative study area. We assessed the effects of sampling density and study area size on landscape genetic inferences for a dispersal‐limited amphibian, Plethodon mississippi, via analysis of nested datasets. Microsatellite‐based genetic distances among individuals were divided into three datasets representing sparse sampling across a large study area, dense sampling across a small study area, or sparse sampling across the same small study area. These datasets were a proxy for gene flow (i.e., the response variable) in maximum‐likelihood population effects models that assessed the nature and strength of their relationship with each of five land‐use classes (i.e., potential predictor variables). Comparisons of outcomes were based on the rank order of effect, sign of effect (i.e., gene flow resistance vs. facilitation), spatial scale of effect, and functional relationship with gene flow. The best‐fit model for each dataset had the same sign of effect for hardwood forests, manmade structures, and pine forests, indicating the impacts of these land‐use classes on dispersal and gene flow in P. mississippi are robust to sampling scheme. Contrasting sampling densities led to a different inferred functional relationship between agricultural areas and gene flow. Study area size appeared to influence the scale of effect of manmade structures and the sign of effect of pine forests. Our findings provided evidence for an influence of sampling density, study area size, and sampling effort upon inferences. Accordingly, we recommend iterative subsampling of empirical datasets and continued investigation into the sensitivities of landscape genetic analyses using simulations.


| INTRODUC TI ON
In the field of landscape genetics, features of the environment that may affect dispersal and gene flow represent potential predictor variables that, once converted to "ecological distances" among local sites, can be compared to corresponding genetic distances among individuals or local populations sampled at those sites (Manel et al., 2003). These comparisons aim to identify abiotic and biotic characteristics of the landscape that influence genetic connectivity, such as riverine or road barriers (Hartmann et al., 2013), agricultural land-use practices (Costanzi et al., 2018;Goldberg & Waits, 2010;Prunier et al., 2014), or the spatial configuration of preferred or nonpreferred habitat (Vergara et al., 2017). Landscape genetics studies have spanned a broad range of spatial scales (e.g., a study area <40 km 2 for the Natterjack toad, Epidalea calamita [Cox et al., 2017] vs. 250,000 km 2 for the greater sage-grouse, Centrocercus urophasianus [Row et al., 2015]). In addition to study area size, the spacing between locations where samples are collected (i.e., sampling density) can affect inferences about which ecological variables have the greatest impact on gene flow, as can the scale(s) at which landscape data are summarized (Anderson et al., 2010;Richardson et al., 2016).
Thus, study area size and sampling density are important aspects of landscape genetics study design (Anderson et al., 2010;Landguth & Schwartz, 2014;Seaborn et al., 2019). Indeed, questions about sampling design that maximize power for detecting relationships between predictor and response variables apply broadly to field-based ecological research (Legendre et al., 2002).
Decisions about final sampling density used in landscape genetic studies are often made a priori. While pilot studies have great value (particularly those outside of the study area, so as to avoid pseudoreplication), they are rarely conducted. As a result, many sampling decisions are instead based on an understanding of the focal species' dispersal ability, but this may be unavailable to some researchers, and so sampling may simply be conducted opportunistically (Anderson et al., 2010). Typically, species with short dispersal distances are sampled more densely than those for which long-distance dispersal is common. For example, Peterman et al. (2014) sampled salamanders at locations as close as 75 m apart, whereas Parks et al. (2015) sampled mountain goats 10-25 km apart. Given that fieldwork can be labor-intensive and expensive, in the case of species with limited dispersal, researchers may be faced with the decision to either conduct their study within a small area so as to avoid exceeding average dispersal distances, or increase the spacing between sampling locations to span steeper gradients of habitat heterogeneity that may be necessary to answer the research questions at hand. For instance, to test the effects of roads on gene flow, ideally, a study area would encompass a moderate to large number of roads of varying sizes (Atzeni et al., 2020;Keller et al., 2014;Richardson et al., 2016)-a condition that may not be satisfied within a small study area. Such decisions about sampling are not trivial, as it is important to design studies that are large enough to detect relevant spatial patterns, yet not so large that these patterns are lost to noise Legendre et al., 2002). Accordingly, there is a need for investigations that explicitly evaluate the impacts of alternative spatial arrangements of sampling sites, using the same focal species, genetic marker set, and landscape setting. This design would enable an examination of the effects of sampling density and study area size on landscape genetics inferences without introducing confounding variables related to intrinsic differences in life-history traits or mutational processes affecting genetic marker polymorphism, or extrinsic regional differences in abiotic variables and biotic interactions.
Sparse sampling (i.e., where average distances among collection sites far exceed typical movement of individuals) may fail to capture relationships between landscape features and gene flow that are more pronounced for short-distance dispersal events, such as localized effects of small water bodies. For example, this was seen by Angelone et al. (2011) in a landscape genetics study of the European tree frog (Hyla arborea) that assessed the impact of different sampling densities. In their study, the authors separately analyzed pairwise comparisons of breeding ponds that were subdivided into contrasting geographic distance classes. This tiered reanalysis identified different ecological predictor variables as being associated with resistance to gene flow for each distance class examined. The authors considered this outcome to be consistent with the notion that rivers or lakes affected short-distance dispersal, whereas geographic distance, wetlands, hedgerows, and the density of forests more strongly affected long-distance dispersal. Another concern about sparse sampling is that it can result in a weaker relationship between genetic distance and the ecological variables that impact gene flow due to stochastic events (e.g., localized extreme weather, interactions with locally invasive or predatory species, and fine-scale disease spread) that may have accumulated effects between distant sampling sites (Epperson, 2010). For instance, Keller et al. (2013) compared the strength of correlation between genetic and ecologically informed geographic distances among populations of wetland grasshoppers (Stethophyma grossum). The strongest fit was found when pairs of sampling locations within close proximity to one another were included in the analyses (i.e., only those up to 3 km apart, the threshold for minimum population connectivity in that study system). Furthermore, the authors found a decrease in model fit when examining widely separated populations, which they suggested may occur because the rarity-and associated stochasticity-of long-distance dispersal events reduced the ability to detect relationships with ecological variables over such scales. Ultimately, both groups of researchers recommended that sampling locations be spaced no further apart than the average dispersal distance of members of the focal species (also see Anderson et al., 2010).
Just as there may be negative consequences for landscape genetic inferences when sampling is too sparse, this may also be true for very dense, fine-scale sampling. The logistical trade-off in which dense sampling is coupled with a smaller study area may create a situation where a representative range of values of one or more ecological predictor variables is not captured by the study design.
In such cases, it may become difficult to identify an environmental variable's true impact on gene flow (Keller et al., 2014). This idea is supported by outcomes from a study by Haran et al. (2017), who repeatedly subsampled their dataset of individual-based genetic distances among pine sawyer beetles (Monochamus galloprovincialis) from the Iberian Peninsula. Through a series of partial Mantel tests examining the association between gene flow and environmental variables, they found significant relationships were more likely to be detected when using larger study areas. Briefly, the authors assessed the relationship between environmental variables and genetic distances in over 30,000 alternative demarcations of their study area, ranging from 220 to 1,000 km in diameter. The number of significant relationships between gene flow and high elevation, cooler temperatures, and pine forests was greatest when study area sizes were large (1,000 km diameter), whereas the number of relationships between gene flow and cooler temperatures was highest when study areas were smaller (600 km diameter). While that study highlighted the potential importance of large study areas, its findings are also consistent with Angelone et al. (2011) and Keller et al. (2013) in supporting the idea that the impact of environmental variables on gene flow may be scale-dependent. Thus, there are reasons for concern regarding overly small study areas, as these may limit researchers' ability to detect biologically meaningful relationships between ecological predictor variables and gene flow (also see Anderson et al., 2010).
In the present study, we explored the impact of sampling density and study area size upon landscape genetic inferences for a lowmobility amphibian, the Mississippi slimy salamander, Plethodon mississippi (Highton, 1989). Insights into the effects of sampling design and effort allocation were generated by assessing the relationship between each of five ecological predictor variables (i.e., land-use types) and gene flow, as measured using individual-based microsatellite genotypic data, with three sampling schemes: dense sampling across a small study area, sparse sampling across the same small study area (i.e., a reanalyzed subset of the former), and sparse sampling across a large study area that encompassed (and exceeded) the boundaries of the small study area (Figure 1). For each sampling scheme, we created and then ranked by goodness of fit a series of maximum-likelihood population effects models (MLPE), a type of linear mixed-effect model that accounts for the nonindependence of pairwise comparisons (Clarke et al., 2002). By comparing the best-fit models for each of the three sampling schemes, we addressed two questions: (1) Does sparse sampling fail to identify relationships between ecological predictor variables and gene flow that are identified using dense sampling? (2) Are the differences (if any) between inferences obtained from the two contrasting sampling densities a consequence of sampling density alone, or does study area size also play a role? F I G U R E 1 (Preferred choice for graphical table of contents) Sampling locations of Plethodon mississippi within Holly Springs National Forest (HSNF), Mississippi, USA. (a) In the sparse sampling across a large study area (630 km 2 ), 19 sampling locations (squares with black dots and circles with black dots) were spaced approximately 7 km apart across the entirety of HSNF. Within the small study area (256 km 2 ) demarcated by a dashed box, 14 sampling locations (plain black dots) were spaced approximately 3 km apart. The circles with black dots within the small study area indicate locations that were also part of both the large and small study area dataset (inset: map of southeastern United States showing location of HSNF). (b) The dense sampling across a small study area included pairwise genetic distances between individuals from all sampling locations. (c) The sparse sampling across a small study area only included pairwise genetic distances from individuals that were >7 km apart

| Focal species and landscape setting
Plethodon mississippi is found within bottomland hardwood and wet pine forests of Mississippi, Alabama, and Tennessee, USA, and develops terrestrially without the need to disperse to aquatic environments for reproduction (Petranka, 1998

| Study design
To examine the effect of sampling density and study area size upon landscape genetic inferences, nested datasets were created. These spaced approximately three kilometers apart, across a 16 × 16 km area nested within the aforementioned large study area, with very similar land-use composition ( Figure 1b). Third, the "Sparse-Small" dataset was a subset of the Dense-Small dataset that included individuals from all 14 sampling locations, but only (and all of) those pairwise comparisons between individuals that were greater than seven kilometers apart, thereby matching the distance between sampling locations in the Sparse-Large dataset (leading to 2,375 individualbased pairwise comparisons; Figure 1c).

| Interpretative framework
Comparisons of similarity among landscape genetic inferences drawn from three nested datasets were used to address our research questions about the relative importance of sampling density and study area size (see Introduction). However, even in the case of very strong similarity, there is little reason to require this outcome to be reflected by identical best-fit models. Accordingly, we focused on using four criteria for which comparisons of the best-fit models from the three datasets were possible. These were as follows: (1) the rank order of importance of five land-use classes, (2) the sign of effect upon gene flow, (3) the optimized spatial scale of effect of each landscape variable, and (4) the optimized transformation of each landscape variable. Together, these four criteria enabled objective and meaningful (albeit broad) assessments of similarity among best-fit landscape genetic models.
Given the three alternative sampling strategies that we considered, for any one of the four criteria of comparison, there are five potential outcomes regarding similarity vs. dissimilarity. First, the effects seen for a given ecological predictor variable may be consistent across all three datasets ( White boxes indicate outcomes are similar to other white boxes, and gray boxes indicate outcomes are dissimilar to both white and gray boxes. A) Similar outcomes are obtained for all three datasets. B) Sparse and dense sampling across a small study area yield similar outcomes, but differ from the sparse sampling across a large study area. C) Sparse samplings across a large and small study area yield similar outcomes, but differ from dense sampling across a small study area. D) Sparse sampling across a large study area and dense sampling across a small study area yield similar outcomes, but differ from sparse sampling across a small study area. E) All models differ. A combination of the effects seen in B) and C) may be the cause of D) or E).
TA B L E 1 Different combinations of sampling density and study area size, and hypothetical outcomes relating to similarity of landscape genetics inferences among datasets outcomes are insensitive to sampling density and study area size.
In this case, robust landscape genetic inferences could be obtained from the most logistically tractable study design (i.e., Sparse-Small), thus providing management-relevant insights at relatively low expense. A second potential outcome would be a case where sparse sampling in both large and small study areas gave similar inferences, but these strongly differed from those obtained from dense sampling in a small study area (Table 1B). This would suggest that sampling density primarily drives outcomes. Here, the densitydependent nature of conclusions drawn from a landscape genetic study would mean that the sampling design needs to be tailored to the research question from the outset, based on a priori knowledge.
In the absence of such knowledge, an alternative approach would be to sample densely but then reanalyze successively "thinned" datasets to understand the magnitude of impact. A third potential outcome could be where the small study area generated similar inferences irrespective of sampling density, yet these inferences strongly differed from those based on sparse sampling in the large study area (Table 1C). This would suggest that study area size primarily drives outcomes, and as above, the scale-dependent nature of conclusions would need to be accounted for in the initial sampling design, or explored via reanalysis of successively scaled-down datasets. The fourth possible outcome would be one in which sparse sampling across the large study area and dense sampling in the small study area produced closely matching inferences, but these contrasted with those obtained from sparse sampling in the small area (Table 1D). We would interpret this to indicate a critical influence of overall sample size (i.e., the total number of individual-based pairwise comparisons), for which there may be threshold effects (i.e., above vs. below n pairwise comparisons). If a minimum adequate sampling "effort" exists, this would need to be identified and exceeded, perhaps necessitating a pilot study. The fifth potential outcome would be a scenario where the effects of a given ecological predictor variable strongly differed across all three datasets (Table 1E). This would indicate that a combination of the aforementioned phenomena (i.e., sampling density-, study area size-and/or effort-dependent outcomes) affected landscape genetic inferences. This latter scenario would suggest that our ability to perform meaningful comparisons across studies is limited, as independent studies will inevitably differ in a least one (or more) aspects of sampling design.

| Geographic sampling
To determine where to place the boundaries of the small study area so that it closely represented the land-use composition of the large (i.e., forest-wide) study area, we first created a map partitioned into six land-use classes (i.e., agriculture, hardwood forest, pine forest, manmade structures, water bodies, and wetlands) through a super-  Table S1). The supervised classes were chosen to represent land uses that may impact the dispersal of P. mississippi. Bottomland hardwood and wetland areas are the primary habitat for this species (Petranka, 1998). Agricultural areas, manmade structures, and streams have all been associated with resistance to gene flow among populations of some salamander species (Costanzi et al., 2018;Marsh et al., 2007). Finally, putatively suboptimal habitat, such as that provided by pine forests and ridges, has been found to have a relationship with gene flow in a closely related salamander, P. albagula (Peterman et al., 2014). A 16 × 16 km square polygon shapefile was created in ArcGIS and then moved across the classified raster file until the land-use class percentages were as close as possible (i.e., within 6.3%) to those found in the large study area (Table S2). As a consequence of optimal placement of the 16 × 16 km square, three sampling locations were shared by the large and small study areas. Accordingly, individuals from these three locations were included in all datasets. To identify differences in the configuration of patches in the small and large study areas, FRAGSTATS v. 4.2 (McGarigal et al., 2012) was used to calculate the patch density, patch cohesion, and correlation length (i.e., a measure that represents the average distance that a randomly placed individual could move without leaving a patch) for each land-use class (Table S3).

| Genetic sampling and analyses
At each sampling location, tail tissue was collected from at least five P. mississippi individuals and then stored in 95% ethanol. In total, 183 individuals from 33 locations were sampled ( Figure 1). Landscape genetic analyses were conducted using an individualbased genetic distance measure (Shirk & Cushman, 2014). Following Shirk et al. and's (2017) recommendations for when sample sizes and genetic structure are low, an individual-based principal components analysis (PCA) was performed using the ade4 package (Dray & Dufour, 2007) in R. Sixty-four axes of ordination were chosen to achieve a balance between explaining a large amount of variance in the data while also attempting to avoid issues of high dimensionality, where the cost of additional axes outweighs the benefits (Beaumont et al., 2002). The final pairwise genetic distance between all sampled individuals was then calculated using the Euclidean distance between the 64 PCA axes.
To determine the extent of any spatial autocorrelation, we estimated a Mantel correlogram (Oden & Sokal, 1986) for the pairwise genetic distance matrix described above against a straight-line (Euclidean) distance matrix using the vegan (Oksanen et al., 2019) package in R.

| Landscape analyses and model testing
The classified land-use raster (see Geographic Sampling¸ above) was used to create a series of rasters for each land-use class using a moving window analysis in FRAGSTATS v. 4. 2 (McGarigal et al., 2012).
Moving windows were square, and sizes were designated using the length of one side; thus, a 250 m moving window encompassed 0.0625 km 2 . The value of each pixel in a raster was determined by the percent of the surrounding window that contained a given land-use class using the PLAND function. Five rasters were created for each land-use class, with moving window sizes of 100, 250, 500, 750, and 1,000 m. To test for nonlinear relationships between PLAND and gene flow, each of the five rasters per land-use class was transformed into eight different rasters that encompassed a range of possible positive and negative exponential relationships (see Figure S1) using the ResistanceGA package (Peterman, 2018)  between any two land-use classes was >0.5 (indicating a high level of correlation between the random-walk distances generated for a pair of land-use classes), only one of the two land-use classes was retained. To isolate the effect of each land-use class on gene flow, the random-walk distances for each raster were regressed against the homogenous landscape distance using a simple linear regression.
The residuals were used in model testing, thereby removing the effect of straight-line geographic distance.
We analyzed three nested genetic datasets, each associated with a different sampling scheme (see Study Design¸ above). For each of these, we separately optimized for the best-fit spatial scale and transformation of each land-use class using a series of MLPE models. These are a form of random effects model that are considered robust for individual-based comparisons (Shirk et al., 2018).
Pairwise genetic distances between individuals collected from the same sampling location were removed to prevent skewing models. Using the optimized scale and transformation for each land-use class, five univariate models (i.e., each land-use class alone) and five multivariate models (i.e., combinations of two or more land use classes) were generated to examine how these contributed to variance in genetic distance ( Table 2). Each model also included the homogenous landscape distance, as simulations have found that the inclusion of geographic distance in MLPE models increases their model accuracy (Row et al., 2017). Models were ranked using AICc.
The best-fit models were examined to determine the sign of effect for each land-use class, where a negative sign of effect indicates gene flow facilitation, and a positive sign of effect represents gene flow restriction (Row et al., 2017).

| Genetic sampling and analyses
Two microsatellite loci (i.e., PG_QWZ and PG_241) showed evidence of homozygote excess in the three putative populations that were each composed of 9-10 individuals sampled from the same site.
However, these loci did not consistently depart from HWE using either of the other two grouping schemes (i.e., when treating all 183 individuals as members of a single population, or when recognizing eight putative populations based on pooling of sites <20 km apart).
The same two loci showed some indications of null alleles, but the frequency of nulls was very low (i.e., <0.1). Accordingly, we retained all eight loci in downstream analyses. Based on 183 P. mississippi individuals, the number of alleles per locus ranged from 6 to 32, with a mean of 16.5 alleles per locus (Table S6), and the amount of missing genotypic data was 1.8%. The genetic neighborhood size, as shown by the largest distance class at which genetic distance and geographic distance were positively and significantly correlated, was 11.5 km ( Figure 2; Mantel's r = 0.037, p = .036). Consistent with these results, the semivariogram plateaued once, at approximately 10 km (Figure 3).

| Landscape analyses and model testing
Tests for nonindependence of potential predictor variables revealed a strong positive correlation (r > 0.9) between the wetland and water land-use classes. Because the primary habitat for P. mississippi is wetland areas, the water land-use class was omitted, and we performed all subsequent analyses using the five remaining classes. For all three nested datasets (i.e., Sparse-Large, Dense-Small, and Sparse-Small), the full model had the lowest AICc and was therefore considered the best-fit ( Table 2).
There were three cases consistent with the interpretation that landscape genetic inferences are insensitive to sampling density and study area size (Table 3A), and all of these related to sign of effect.
Specifically, across all three datasets, hardwoods and wetlands were correlated with resistance to gene flow, and manmade structures were correlated with facilitation of gene flow (Tables 3 and 4). The only instance of a sampling density-dependent outcome (Table 3B) was for the optimized transformation of the effect of agricultural areas, which was the same (i.e., inverse-reverse Ricker) for the Sparse-Large and Sparse-Small datasets (Table 4). Our results were consistent with study area size-dependent outcomes (Table 3C)

F I G U R E 3
Semivariogram created using pairwise genetic distances and geographic Euclidean distances using 52 distance classes with a distance interval of 1.5 km. The plateau at approximately 10 km indicates that this is the spatial scale over which spatial autocorrelation is the strongest. N size denotes the number of pairwise comparisons within the given distance class datasets, yet differed from the Sparse-Small dataset (Table 5). Finally, combined effects of sampling density, study area size, and/or sampling effort were the most common outcome observed across our best-fit landscape genetic models (Table 3E). For example, for rank order of effect, all land-use classes except for wetlands showed this pattern.
Likewise, for the optimized spatial scale and transformation criteria, hardwoods and pine land-use classes consistently exhibited unique differences among all nested three datasets (Table 4).

| D ISCUSS I ON
Landscape genetics is a relatively young subdiscipline, and as such, knowledge gaps remain (Richardson et al., 2016). One gap that is common within the field of spatial ecology as well as landscape genetics is an understanding of the impacts that a priori choices regarding sampling density and study area size have upon inferences (Seaborn et al., 2019). Indeed, consideration not only of total sampling effort, but also allocation of that effort in light of the trade-off between overall spatial extent vs. local density, has been a major focus of recent simulation-based research (e.g., Landguth et al., 2012;Landguth & Schwartz, 2014). In this study, we used a comparison of landscape genetic models resulting from three alternative sampling strategies for a low-mobility salamander, P. mississippi, to generate insights from datasets that reflect the inherent "noise" embedded in empirical studies. Given that we used ecological predictor variables that were optimized for spatial scale of effect and transformation separately for each dataset, this potentially leads to final best-fit models with several nonidentical features. As such, there can be many nuanced differences among models, some of which may not directly reflect the impacts of sampling density or study area size.
Accordingly, our comparison of models necessarily focused at a relatively coarse level (i.e., primarily the rank ordering of importance

| Outcomes insensitive to sampling density and study area size
For all three genetic datasets, there was a correlation between the presence of wetlands and resistance to gene flow; the same correlation was also found with hardwood forests. Given that P. mississippi individuals often reside in bottomland hardwood forests and wetland areas (Petranka, 1998), it may seem counterintuitive that preferred habitats are associated with resistance to gene flow among populations. However, based on empirical and simulation studies TA B L E 3 Visual representation of the four criteria: rank, sign (i.e., correlation with gene flow facilitation or resistance), scale, and transformation (labeled "trans."), for each of the five potential outcomes of congruence detailed in Table 1 Each "X" indicates the given landscape variable exhibits congruence at the given axis of comparison in the manner described by the hypothetical outcome listed. For example, hardwoods were found to correlate with resistance to gene flow (i.e., the same sign) in all three datasets; thus, the result is consistent with outcome A, where inferences are found to be insensitive to both density and size. (Keely et al., 2016(Keely et al., , 2017, some researchers have reported this same outcome, presumably attributable to individuals choosing to stop or stay in high-quality habitat areas. In our study, another outcome that was consistent across all three datasets was a correlation between manmade structures and gene flow facilitation. Notably, roads were a significant component of the manmade structure land-use class in our analyses. Although avoidance of road edges by many amphibians has been documented, Marsh and Beckman (2004) found no effect of forest roads on the presence of P. glutinosus (the sister species of P. mississippi), supporting the idea that slimy salamanders may move freely across these types of open substrates. This may extend more broadly, given that Prunier et al. (2014) found a correlation between roads and increased gene flow among populations of the alpine newt (Ichthyosaura alpestris). As a corollary to counterintuitive idea of slow movement through high-quality habitat discussed above, an additional interpretation warrants consideration here: individuals may choose to depart or accelerate through low-quality habitat areas, such that apparent conduits for gene flow may not necessarily be places that also provide protection and/or harbor resident populations.

| Sampling density-dependent outcomes
Based on the design of our study, differences in landscape genetic outcomes that were restricted to comparisons between sparse vs.
dense sampling would indicate density dependence of inferences.
In some species, dispersal can be divided into two categories: a collection of several small, "routine" movements (i.e., those associated with recurrent processes such local mate-seeking, or foraging, etc.), or conversely, a singular "special" movement (i.e., directed dis-

TA B L E 5
Rank of effect of landscape variables in best-fit maximum-likelihood population effects models for each dataset field-based experimental study that focused on the ringed salamander, Ambystoma annulatum, Ousterhout and Semlitsch (2018) identified individuals that traveled only short distances as "residents," or salamanders that traveled larger distances as "dispersers." These groups of individuals were found to move differently across the landscape, with dispersers being impacted by habitat type when moving, while residents were not. These contrasting movement behaviors could manifest as different landscape genetic inferences, depending on whether dense vs. sparse sampling was used. However, in the present study, we found only one relatively minor case where this occurred: For the agriculture land-use class, the Sparse-Large and Sparse-Small datasets were both optimized with an inverse-reverse Ricker transformation, whereas the Dense-Small dataset had an inverse Ricker transformation (Table 4). Interestingly, simulations conducted by Landguth and Schwartz (2014) found that for continuously distributed individuals with limited dispersal, sparse sampling risks overestimating the true magnitude of genetic structure, whereas dense sampling tends to lead to underestimating it. Those authors considered the most likely reasons for this to be either false hierarchical structure in the former case, or inadequate characterization of local variability in the latter. While these issues are certainly noteworthy, given the apparent rarity of sampling densitydependent outcomes for P. mississippi, in our study system, sampling density alone may be less of a concern than initially thought.

| Study area size-dependent outcomes
Within our study's comparative framework, differences between inferences generated over a large study area (i.e., using the Sparse-Large dataset) and those from over a small study area (i.e., using the Dense-Small and Sparse-Small datasets) suggest that contrasts in study area size may be consequential for landscape genetics studies. One reason for this may simply be that the limited extent of any small study area increases the probability that dispersal and gene flow from outside areas (i.e., unsampled populations that make nonnegligible genetic contributions to the observed data) may obscure or even overwrite signatures of facilitation or resistance by landscape features within the focal area (Anderson et al., 2010 it may not only affect dispersal probabilities, but also time to fixation of alleles, and thus, detection of such events van Strien et al., 2015). In a study of American pika (Ochotona princeps), Castillo et al. (2016) found that patch connectivity of a given environmental predictor variable impacted its relationship with gene flow. Other studies have also shown that habitat patch characteristics such as density, cohesion, and correlation length can drive different outcomes in replicated or nested landscape genetic analyses (Burgess & Garrick, 2020;Cushman et al., 2012Cushman et al., , 2013Vergara et al., 2017).
When comparing patch configuration metrics, we found that the patch density of pine forests within the small study area was much lower than within the large study area (Figure 4, Table S6). Although patch density can be affected by the size of the area in which it is calculated (McGarigal et al., 2012), in the present study, the difference in pine patch densities was far greater than the differences among the patch densities of agricultural areas, hardwood forests, manmade structures, and wetlands. These findings indicate the effect of pine patches on gene flow among P. mississippi populations may be dependent upon the density of those patches across the landscape.

| Sampling effort-dependent outcomes
Strong similarities between one or more components of the bestfit landscape genetic models that were seen only in the Sparse-Large and Dense-Small datasets may be indicative of the effect of overall sampling effort. This could be associated with a threshold effect, where being above vs. below a certain total number of pairwise comparisons drives the result. In our study, such a threshold might lie between 2,375 and 3,916 pairwise comparisons (Table 1).
We identified six cases where the pattern of similarity across datasets was consistent with this scenario (Table 5). Of these, the most notable relates to sign of effect: whereas the Sparse-Large and Dense-Small datasets showed a correlation between agriculture and gene flow restriction, this relationship switched to gene flow facilitation in the Sparse-Small dataset (Table 5). However, although the correlation length of agricultural patches in the large and small study areas was similar, composition of agricultural areas differed. Therefore, this is a potentially confounding factor.
Indeed, our agricultural land-use class was composed of an array of crops, including soybeans, cotton, sweet potatoes, and corn. It is possible that the crops grown within the small study area are not representative of those within the large study area, and the relationship between gene flow and agriculture might be dependent upon crop type. That said, given there were five additional examples of sampling effort-dependent outcomes, we consider the potential for threshold effects of total sample size to be an important consideration. Likewise, simulations by Landguth et al. (2012) showed that the total number of sampled individuals affects the power for detecting impacts of landscape features on dispersal and gene flow for continuously distributed species, although that study modeled an organism with moderate dispersal abilities and used a different inferential framework (i.e., causal modeling) than we did here. Accordingly, we refrain from suggesting that total sampling effort and associated threshold effects are a common phenomenon for low-mobility species such as salamanders, but it does warrant further investigation.

| Combined impacts of sampling density, study area size, and/or sampling effort
Dissimilar outcomes across all three datasets were very common when considering rank ordering of land-use classes (four cases), and this was also quite common for the optimized scale and transformation included in best-fit landscape genetic models (two cases each;  . However, we suggest that successive reanalysis of subsampled empirical datasets (i.e., those with "thinned" sampling densities, and/or downscaled study area sizes), coupled with metareplication (e.g., Burgess & Garrick, 2020;Castillo et al., 2016;Short Bull et al., 2011), may alleviate some of these issues.

| Limitations
In our nested design, we did not include dense sampling across a large study area (i.e., a Dense-Large dataset) owing the prohibitive F I G U R E 4 Map illustrating the distribution of the "pine" land-use class across large and small study areas in Holly Springs National Forest. Burgess and Garrick (2020). The effect of sampling density and study area size on landscape genetics inferences for the Mississippi slimy salamander, Plethodon mississippi. Ecology and Evolution size of HSNF. Nonetheless, this would have enabled a more complete exploration of impacts of trade-offs between different kinds of effort allocation. This is something that can be addressed with simulations (e.g., Landguth & Schwartz, 2014), but we also see value in nested analysis of empirical genetic data that are more realistic in their inherent "noise" (i.e., genotype matrices that may include some or all of the following: alleles from de novo mutations, some genotype scoring error, weak nonindependence among loci, individuals that are not the product of strict outcrossing, and/ or departure from 1:1 operational sex ratios). Likewise, we did not include a very fine-scale sampling strategy in which geographic spacing among sampling locations was equal to or less than the average dispersal distance of individuals. This was because a feasible number of sampling locations would have encompassed a very small area that was unlikely to include all five land-use classes.
However, such fine-scale sampling has been recommended (e.g., Anderson et al., 2010). Perhaps the most likely impact of this omission from our study is an inability to make inferences on short, near-contemporary timescales. This is because the inhibiting or facilitating effect of landscape features that arose very recently (or are cyclical) upon direct dispersal of individuals from one location to another requires sampling those spatially proximate locations.
Conversely, when sampling is more sparse than average dispersal distances, most of the signal comes from the accumulated multigenerational effects of dispersal and gene flow (i.e., many contributions from unsampled intermediate locations), such that only more permanent or long-standing landscape features are likely to register (van Strien et al., 2015). Although this may be a missed opportunity in the present study, our exclusive focus on "historical" connectivity (cf. a mixture of contemporary and historical, depending on sampling scheme) may have set up fairer comparisons across different data subsets. Nonetheless, taken together, the aforementioned practical limitations did constrain the number of possible combinations of sampling density, study area size, and sampling effort. We were also limited to only one subset area (cf. multiple replicates), again due to considerations about the logistics of sampling. This is a general restriction of most empirical studies. Indeed, to date, landscape genetic methods have mostly been explored using simulations (e.g., Cushman & Landguth, 2010;Landguth et al., 2012;Landguth & Schwartz, 2014;van Strien et al., 2015), and this includes sampling strategy optimization (Selmoni et al., 2020).
While empirical studies can assess interactive and additive forces (Resasco et al., 2017;Seaborn et al., 2019), their use in testing the effect of different methods and/or sampling strategies does have unavoidable limitations. For example, geographic replicates are never identical, and as was the case in the present study, multivariate models created using optimized predictor variables (e.g., landuse classes) often include predictor variables of different scales and transformations (Castillo et al., 2016;Vergara et al., 2017).
We also limited our analyses to assessing the relationship between five land-use classes and gene flow. There are other environmental variables that may also have relationships with gene flow; however, their addition may potentially increase the number of differences between geographic replicates. While there is certainly scope for improvement, given the importance of optimizing "scale" (see Anderson et al., 2010 for synonyms), the present study serves as a useful case study for identifying real-world study systems with characteristics that are conducive to a nested sampling design, and for which iterative landscape genetic analyses can be informative about effort allocation and inference sensitivity.

| CON CLUS IONS
The present study revealed an encouraging finding: the directionality (sign) of effect of a given land-use class on gene flow-whether it acts as a facilitator vs. inhibitor-was commonly found to be insensitive to sampling density and size of the study area. Indeed, this was the case for three of the five land-use classes considered (i.e., hardwood forests, manmade structures, and wetlands; Table 2).
Promisingly, the use of more densely distributed sampling locations (i.e., the Dense-Small dataset) did not lead to changes in the directionality of the relationship between land-use classes and gene flow.
This suggests that sampling locations for continuously distributed salamander species may be allocated relatively sparsely across large study areas to capture greater environmental variation without failing to identify relationships between environmental variables and gene flow that are seen when using denser sampling. That said, for the three remaining criteria that we used to determine similarity among best-fit landscape genetic models (i.e., rank ordering of importance, and optimized scale and transformation), we saw evidence for a strong influence upon outcomes of either total sampling effort (coupled with potential threshold effects), or combined effects of sampling density, study area size, and/or sampling effort (Tables 1   and 2). Aside from outcomes of maximum incongruence across nested datasets (Table 3E), we found the greatest number of dissimilarities when considering the effect of overall sampling effort (i.e., total number of pairwise comparisons) on model inferences.
Indeed, there may be a need for routine subsampling and reanalysis of empirical datasets (i.e., thinning) to understand the sensitivity of inferences to these aspects of sampling design. #0324164. We also thank the U.S. Forest Service for scientific collection permits.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest. The findings and conclusions in this article are those of the authors and do not necessarily represent the views of the U.S. Fish and Wildlife Service.

DATA AVA I L A B I L I T Y S TAT E M E N T
Genotypic data, R code, and land-use rasters are available from DRYAD Repository entry found at https://doi.org/10.5061/dryad. h9w0v t4gt