Microhabitat Use by Disjunct Common Five-Lined Skink (Plestiodon fasciatus) Populations at Their Northwestern Range Extent Open Access

MICROHABITAT suitability is particularly critical for organisms at the periphery of their geographic ranges (Quirt et al. 2006) because of its role in biological processes such as thermoregulation, refugia seeking, reproduction, and prey acquisition. Therefore, knowledge of habitat associations is critical for species conservation and management needs. Habitat loss and degradation are driving forces in population extirpation and can intersect with isolation effects and species’ intrinsic factors to place small, disjunct, or peripheral populations at greater risk (Hardie and Hutchings 2010; Heinrichs et al. 2016; Kimmel et al. 2022). These populations may have reduced genetic variability, potentially an additional risk factor (Noss and Cooperrider 1994; Wick 2004; Howes and Lougheed 2008). Yet peripheral populations also may be reservoirs of genetic diversity valuable for species adaptation under climate change (Lesica and Allendorf 1995; Macdonald et al. 2017).

An often broadly used term, “habitat” is not synonymous with vegetative type or land cover, but rather is ultimately organism specific (Krausman and Morrison 2016). Site-management success relies on identified objectives such as responses of a target species as evidenced by their occupancy, selection, or preferences (Török and Helm 2017). Reptiles may be particularly sensitive to changes in habitat because of life-history traits such as ectothermy, microhabitat requirements, and relatively restricted dispersal abilities (Gibbons et al. 2000; Olson and Saenz 2013; Tan et al. 2023). Nevertheless, conservation and research of reptiles are historically underfunded and underrepresented, despite reptilian ecological importance (Gibbons 1988; Todd et al. 2010; Christoffel and Lepczyk 2012). In part, this may be because of traits that complicate their study (Crane et al. 2021). Reptiles also have lesser social capital relative to charismatic bird and mammal species (Olson and Pilliod 2022). Those taxa may not serve as good proxies for reptile conservation, particularly at-risk reptile species (Lindenmayer et al. 2014; Ding et al. 2016; Doherty et al. 2020). Therefore, collecting microhabitat data for reptiles and identifying management thresholds is important for their conservation.

Although stable elsewhere, Common Five-lined Skinks (Plestiodon fasciatus) are considered vulnerable to critically imperiled in their northernmost range and have legal protection in three states and in Ontario, Canada (NatureServe 2022). Among key risk factors are loss, fragmentation, and degradation of its habitat, including microhabitat (Environment Canada 2013; Hecnar et al. 2018; COSEWIC 2021). In the US Upper Midwest, P. fasciatus occurs in disjunct populations in Minnesota, Iowa, and Wisconsin, which exhibit the lowest genetic diversity across this species’ range (Howes and Lougheed 2008). Throughout most of its range, P. fasciatus is associated with partially open hardwood forests. In contrast, it is an endemic representative of rock outcrops in the prairies of western Minnesota. These large granite or gneiss outcrops supply open areas for basking as well as crevices for overwintering habitat (Lang 1982), providing thermal and survival advantages in northern latitudes, as this lizard generally functions best at relatively high temperatures (Fitch 1954; Howes and Lougheed 2007). The Upper Minnesota River Valley population is small, fragmented, and disjunct. It comprises the western edge of the species’ range, separated by distance and geology. Its nearest Minnesota counterpart is approximately 215 km eastward, yet preliminary evidence suggests greater similarity in habitat type with the southern Canadian Shield population, located northeast over 800 km away (Howes and Lougheed 2004; Quirt et al. 2006; Minnesota Department of Natural Resources [MNDNR] 2022). Conservation status, outdated occurrence records, large data gaps, and insufficient knowledge of its response to habitat alteration prompted renewed study of this population.

Historically, open prairie outcrops have become engulfed by encroaching trees and shrubs, including native species like Eastern Red Cedar (Juniperus virginiana) and nonnative species like Common Buckthorn (Rhamnus cathartica), because of changes in land practices (e.g., reduced grazing, fire suppression) and potentially also climate change (Marschner 1974; Lang 1982; MNDNR 2019; McClain et al. 2021). Lang (1982) attributed declines in P. fasciatus in the Upper Minnesota River Valley with increased woody encroachment in openings. Management of public conservation lands and selected private tracts within the project area often include management of woody vegetation encroaching into native prairies and rock outcrops. Management actions typically involve removing trees and brush, creating brush piles, prescribed burning, applying herbicides, and grazing. Reptiles adapted to open landscapes including rock outcrops can benefit from such management (Kingsbury and Gibson 2012); however, management nuances may have species-specific effects on wildlife (Cavitt 2000; Pike et al. 2011; Schrey et al. 2016).

The vulnerability of this small satellite population elevates concerns about risk factors and management effects. Acquiring and integrating improved knowledge of specific habitat associations of underrepresented yet representative species like P. fasciatus is vital to help ensure conservation of native biodiversity and ecosystem function. We identified a need to provide management targets for habitat variables that could be applied to benefit P. fasciatus in this region. We had two objectives: (1) assess site-specific habitat associations of this P. fasciatus population to improve future guidance on preferred timing, location, and application of management and land-use changes; and (2) study P. fasciatus responses to habitat management where we had the opportunity to assess sites before and after treatment.

Our study area was comprised of four sites, including three on public lands (Sites 1, 3, and 4) and one on a private conservation easement (Site 2). We chose sites with P. fasciatus populations within the Upper Minnesota River Valley corridor, which is characterized by prominent crystalline bedrock outcrops (Patterson and Boerboom 1999; Fig. 1). During the mid-1860s, this corridor consisted primarily of native prairie, oak (Quercus sp.) openings (oak savanna), and barrens (open bedrock outcrops) with more heavily forested areas mainly concentrated near the river bottoms (Marschner 1974). Since then, historical and more recent aerial photography revealed major land-cover shifts to forested cover, primarily attributed to fire suppression and reduced grazing by large native and domestic herbivores (Lang 1982; MNDNR 2019). This facilitated expansion of cedars and subsequently invasive species such as buckthorn.

In spring 2015, we established grids consisting of 25–42 sample plots spaced 15–30 m apart across rock outcrops at Sites 1–3. We added Site 4 in spring 2016. Each sample plot consisted of paired 0.5-m² cover objects, one plywood board and one corrugated sheet metal, placed within 2 m of each other. We used two different materials to provide skinks options for thermal and moisture conditions. Although we centered grids on rock outcrops, each represented a range of vegetation, terrain, canopy, and ground-cover variables. We systematically checked cover objects from 2015 through 2017 1–2 times per week throughout the active season, roughly April into October between 0600 and 1800 h CDT. Sites were all checked the same number of times, on the same intervals. We rotated check times among sites to avoid a time bias. Technicians hand-captured skinks and collected morphometric and demographic data. In 2016 and 2017 adults were marked using an 8.5-mm passive integrated transponder (PIT) tag (BioMark®). Juveniles too small for the PIT tag were systematically cauterized with a series of marks partitioned in four quadrants on the ventral surface (LeClere 2017a,b). We used permanent markers to mark skinks in 2015 and found it to be an unreliable method.

We divided the sampling period each year into early and late phases, with the late phase beginning in mid-July. These phases helped us account for changes in vegetation through the season, and changes in behavior, detectability, and population size posthatching.

In 2016 and 2017 we measured habitat characteristics at each sample plot (Table 1), within a 3-m-radius plot centered between the tin and board cover objects. We measured habitat characteristics in each phase each year (2016, 20–26 May and 22–24 August; 2017, 22– 24 May and 15– 18 August). We measured canopy closure using spherical densiometers modified with the Strickler method (Strickler 1959; Jennings et al. 1999).

In 2015 we collected habitat data using photographs of the cover objects and surrounding area, and thus did not include those results in habitat association analyses because of the differences in methodology and lower accuracy. We included 2015 P. fasciatus count data, though, in addition to 2016 and 2017 data, to evaluate responses to habitat management.

We changed measurement protocols for down woody debris between 2016 and 2017 because we found the 2016 protocol was insufficiently precise. In 2017 we measured lengths and diameters of woody debris >2.54 cm diameter rather than using ocular estimates of cover as in 2016. In addition, in 2017, we divided down woody debris into size categories because of the potential for different sizes of cover to be used differently by P. fasciatus.

We also gathered data using temperature loggers (Thermochron iButton DS1921 Kit, Maxim Integrated; Table 1). We placed temperature loggers under both wood and metal cover objects at 8–9 sample plots in each site, using stratified random sampling so that the number of temperature loggers was approximately evenly distributed among plots in open areas, along edges, and in heavily canopied areas.

Management projects occurred at each of the four sites over the course of this project. Management targeted removal of woody vegetation, particularly shrubs including buckthorn, red cedar, and smooth sumac. Management consisted of prescribed burns, mechanical removal of woody vegetation, herbicide application, or some combination of the three. These projects addressed woody vegetation throughout the entirety of the rock outcrop systems that we sampled. We do not have direct information on the intensity of the management, but we measured changes in shrub density between 2016 and 2017, and thus have a measure of habitat impact for some of the management.

We used two analytical techniques: N-mixture models (Royle 2004) and generalized linear mixed-effects models (GLMMs; Bolker et al. 2009; Zuur et al. 2009; Harrison et al. 2018). This approach allowed us to assess validity of our results based on differences in results between the techniques (Fieberg and Johnson 2015). Our two techniques had contrasting strengths and weaknesses. N-mixture models simultaneously estimate both detection and abundance parameters. Detection in this context is the probability of observing an individual skink under a cover object during a single survey. To provide estimates of detection, the N-mixture model we used required the assumption of a closed population, meaning there should be no deaths, births, emigrations, or immigrations over the span of data collection. GLMMs did not require this assumption; nor did they separate effects of detection from effects of abundance. We ran these two models in each of four year-phases: early 2016, late 2016, early 2017, and late 2017. We ran separate models for each year-phase, rather than one model including all data, to: (1) accommodate changes in some of the habitat measurements between years, (2) reduce our violation of the closed-population assumption for the N-mixture models, and (3) evaluate habitat associations based on their consistency across the year-phases and models.

Two types of covariates can be added to N-mixture models: abundance and detection. For abundance, our covariates were canopy closure, herbaceous cover, shrub stem density, litter depth, and down woody debris (Table 1; for description of dropped variables, see Supplemental Materials I, available online). We included one detection covariate: cover-object temperature. For the GLMMs, we included the same covariates used in the N-mixture models as fixed effects. We included sample plot (n = 125) within site (n = 4) as a nested random effect to account for repeated surveys and spatial autocorrelation at the site level (plot was the only random effect for the 2017 late phase because of convergence problems). We used R v4.2.3 (R Core Team 2023), to run N-mixture models in the unmarked package (Fisk and Chandler 2011; Kellner et al. 2023) using the pcount command, and to run GLMMs in lme4 (Bates et al. 2015) using the glmer command and Poisson link function (for model specifications see Supplemental Materials II, available online).

For each analytical technique and year-phase we developed only one model with one set of covariates, avoiding model selection procedures and reducing potential for overfitting and biased coefficient estimates (Harrell 2001; Giudice et al. 2012; Fieberg and Johnson 2015). We limited the number of covariates in each model to six for 2016 and seven for 2017, following the rule to limit degrees of freedom to 10% of the effective sample size (Harrell 2001; Harrison et al. 2018). We were able to include one more variable in 2017 because of the greater number of skink observations. We included two different size classes of down woody debris in the 2017 models, instead of the single down woody debris variable we used in 2016. Our first step to reduce the pool of variables was to assess multicollinearity. We ran a basic linear model containing all variables and then dropped variables until all had variance inflation factors <3 (Zuur et al. 2009). Second, we prioritized remaining variables that we hypothesized to be most important for P. fasciatus life-history requirements, based on our own knowledge and existing information in the literature (Dormann et al. 2013). Finally, we retained variables that we expected to be important from a management perspective (Fieberg and Johnson 2015). We dropped some variables considered important for P. fasciatus, but that were unlikely to be directly manipulated by land managers; most notably crevices and loose rock (Lang 1982; Howes and Lougheed 2004; Quirt et al. 2006; Feltham 2020). The habitat variables we retained (Table 1) are commonly managed in this region by prescribed fire, tree harvest, herbicide application, and removal of undesirable woody vegetation. Focusing on management-related variables supported our objective of producing information in a management context. We included cover-object temperature because of its potential to influence skink use of cover objects, and thus help account for detectability.

After running the initial models for the two techniques, we ran two post hoc N-mixture models for each year-phase: (1) in place of shrub stem density we added covariates for shrub stem densities of buckthorn, sumac, and other species; and (2) in addition to canopy closure we added the square of canopy closure to allow for nonlinear response. We used this post hoc analysis to explore potential causal mechanisms and to help define management thresholds.

For each plot and year, we calculated skink observation frequency (number of observations per 10 surveys). We then calculated the mean annual change in observation frequency by site. We compared the direction and magnitude of change between years in which there was management and years when there was no management. This analysis included 2015 data in addition to 2016 and 2017.

We assessed the importance of habitat variables based on model coefficients, including their confidence intervals, magnitude, and consistency across analytical techniques and year-phases. We developed management targets first by calculating N-mixture model predictions using target skink abundances. To calculate the target skink abundance for each year-phase, we: (1) averaged the number of observations across surveys for each plot, (2) divided the mean observations for each plot by the detection rate estimated from the N-mixture model, and (3) calculated the 75th percentile of the detection-adjusted observation rates. We then used the predict function in R package unmarked (Fisk and Chandler 2011), where we manipulated the habitat variable of interest, while holding all other variables at their mean, until the prediction equaled the target skink abundance. We averaged the predicted habitat values among the year-phases where the coefficient confidence interval excluded zero (i.e., was statistically significant). We used these mean habitat values as the basis for developing management targets, adjusting them based on the practical realities of habitat management and results from the post hoc models and management response analysis.

Surveyors recorded 976 total P. fasciatus observations across 6,630 checks of cover-object pairs (Table 2). Observation frequency varied from as low as 0.06 per 10 surveys per plot (Site 1, 2015) to as high as 1.78 per 10 surveys per plot (Site 4, 2017).

Across 2016 and 2017 we recaptured 23 skinks a total of 39 times. Mean movement distance was 9.6 m (SD = 17.7 m). Twelve animals were recaptured at the same plot where they were originally captured, nine showed movements ranging from 13.3 to 32.4 m, and two showed comparatively long movements: 76.5 m between captures 53 d apart, and 75.3 m between captures 26 d apart.

Plestiodon fasciatus abundance decreased with canopy closure, increased with down woody debris, and decreased in sites that had over 75% herbaceous cover (Table 3; Figs. 2 and 3; for detailed model outputs, see Supplemental Materials II). In 2017, after we divided down woody debris into size categories, P. fasciatus density was associated with large down woody debris and showed no association with medium woody debris. The direction of these habitat coefficients (positive or negative) was consistent across year-phases, though the strength of the coefficients was usually not consistent. Confidence intervals of coefficients excluded zero (i.e., were statistically significant) in only two or three of the four year-phases for each of these habitat variables. Skink abundance overall changed little in response to shrub stem density and litter depth, with coefficient values at or near zero in three of the four phases (Table 3). Skinks were less likely to be observed under cover objects as temperature increased under those objects (Table 3).

We drew the same conclusions from this model approach as we drew from the N-mixture models. Associations were the same, and coefficient values were often similar (Table 4).

After dividing shrubs into different species, we found that P. fasciatus abundance decreased with buckthorn stem density and increased with smooth sumac stem density, although only one of four phases had 95% confidence intervals that excluded zero (Table 5; Fig. 4). Skink abundance changed little with density of shrub species other than buckthorn and sumac in three of the four phases.

We found mixed support for nonlinear associations with canopy closure. There were no nonlinear associations in either late 2016 or early 2017 (Table 5). In early 2016 and late 2017, there was evidence of a threshold at approximately 25% canopy closure, below which P. fasciatus abundance either declined less sharply or actually decreased (Fig. 4).

Detection estimates from the N-mixture models were 5–11% among the four year-phases. Detection covariates (with 95% confidence intervals) were early 2016, 0.113 (0.079, 0.205); late 2016, 0.088 (0.053, 0.177); early 2017, 0.051 (0.039, 0.075); and late 2017, 0.089 (0.079, 0.122).

Temperature under the cover objects was consistently and negatively associated with skinks’ use of cover objects. Overall, the effect of cover-object temperature was similar among year-phases and between the two model types, with the exception of early 2016 when it had no effect in the N-mixture model.

The only significant changes in P. fasciatus observation frequency between years coincided with management to reduce woody vegetation (Table 6). Management projects occurred at each of the four sites over the course of this project. Management consisted of prescribed burns, mechanical removal of woody vegetation, herbicide application, or a combination of the three. These four projects were followed by increases of P. fasciatus observation frequency of 188%, 106%, 120%, and 97% at Site 1 (2016–17), Site 2 (2015–16), Site 4 (2016–17), and Site 3 (2015–16), respectively (Table 6). We recorded corresponding changes in shrub stem density between 2016 and 2017 at two of these sites, with decreases of 51% and 22% at Site 1 and Site 4 (Table 6). We did not measure shrub stem density in 2015.

We predicted values of five habitat variables given target skink abundances (Table 7). Three variables (canopy closure, large down woody debris, herbaceous cover) had at least two covariates with 95% confidence intervals that excluded zero (Tables 3 and 4). Sumac and buckthorn stem density each had only one covariate with a 95% confidence interval that excluded zero (Table 5), but we included them because of the large effects we observed in response to shrub removal. These predictions served as the basis for the final management targets we identify in the Discussion.

Studies quantifying canopy closure or cover (Jennings et al. 1999) and P. fasciatus abundance in northern latitudes commonly show a positive relationship between open areas and skink abundance (Brazeau 2016; Feltham 2020). Our findings regarding canopy closure were consistent with past research in the same or similar landscapes indicating P. fasciatus were associated with exposed granite outcrops within complexes of mixed coniferous and deciduous forests (Lang 1982; Howes and Lougheed 2004; Quirt et al. 2006; Feltham 2020). Canopy closure typically is a significant factor influencing P. fasciatus occupancy in northern latitudes where the average annual temperature is lower. Northern populations typically occupy sites with a greater proportion of open canopy and exposed bedrock than do southern populations (Fitch 1954; Quirt et al. 2006; Watson 2008; Watson and Gough 2012). Plestiodon fasciatus likely benefits from a habitat mosaic that provides thermal gradients for seasonal adaptation.

Open-grown oaks (savanna) contribute to those mosaics where stands still occur on bedrock outcrop complexes within our study sites, similar to historical conditions (MNDNR 2005). However, cedars began dominating many outcrops following the droughts and fire suppression of the early 1930s, enclosing many formerly open sites (Lang 1982). Extensive invasions of hill prairies by cedar and buckthorn can result in closed communities with elimination of prairie species and oak savanna (Ugarte 1987; Wyckoff et al. 2012). Recent management to reduce canopy cover has restored open outcrop habitat in the Upper Minnesota River Valley, although long-term maintenance is necessary because of regrowth.

Our finding that large down woody debris was important for P. fasciatus is consistent with other studies (Hecnar 1991, 1994; Hecnar and M’Closkey 1998). Large down woody debris offers opportunities for protection and favorable microclimates. It is structurally more stable and persists longer in the environment than medium or small down woody debris. It also buffers environmental fluctuations such as temperature and moisture (Savely 1939; Boddy 1983; Hecnar 1991). Hecnar (1994) found that large, moderately decayed logs were required for P. fasciatus nest sites on that study area. These logs may also provide foraging opportunities for invertebrate prey (Harmon et al. 2004).

The availability and structure of grass tussocks and other vegetation may be more important as microhabitat providing suitable microclimate, food availability, and protection from predators than previously thought, similar to Northern Prairie Skinks, Plestiodon septentrionalis (Danielsen et al. 2014; Hecnar and Brazeau 2017; Brazeau and Hecnar 2018). Hecnar and Brazeau (2017) found that adding ground-cover variables to canopy and microclimate improved predictive ability of their habitat models for P. fasciatus. Our study sites included herbaceous cover, which may have provided these benefits. However, P. fasciatus densities decreased on sites where herbaceous cover exceeded 75%. These declines may be related to reduced availability of bare rock.

The lack of an association with shrubs in our a priori models was due to the contrasting effects of buckthorn versus sumac, which essentially canceled each other out in the overall effect of shrub stem density. We found this disparity of P. fasciatus response to buckthorn versus sumac stem density using an exploratory, post hoc analysis, which supports the need to evaluate effects of shrub stem density in the presence of invasive buckthorn more formally. The finding is worth exploring because of the important role of woody plant removal in management of this ecosystem. Explanations for the difference in P. fasciatus response to buckthorn versus sumac stem density are uncertain. Common buckthorn is an introduced species, whereas smooth sumac is native to Minnesota. However, both species produce suckers and can form large dense stands when uncontrolled. Anecdotally, the structure of buckthorn and sumac thickets appeared to differ on our sites. Knight et al. (2007) and Charles-Dominique et al. (2012) discussed traits that account for invasion and prolific regeneration of common buckthorn. However, we found nothing in the scientific literature directly comparing seedling densities in the understories of common buckthorn with smooth sumac shrub thickets. Other hypotheses potentially explaining the dissimilar effects of the two shrub species include differing effects on soil chemistry and litter layer properties, and potentially major differences in impacts to the arthropod community on which skinks depend for food (Trial and Dimon 1979; Kollmann and Grubb 1999; Heneghan et al. 2002, 2006). Therefore, we speculate that buckthorn thickets may be less attractive to skinks, in part due to reduced prey availability.

Our study of responses to habitat management reinforces the importance of shrub management for P. fasciatus. There were large, positive responses to woody removal, either mechanical, herbicide, or prescribed fire. This is consistent with published justifications for managing both native and nonnative shrubs in prairies, particularly in more rugged terrain more prone to woody encroachment (Bragg and Hulbert 1976; Petranka and McPherson 1979; Werner and Harbeck 1982; Ugarte 1987). Yet in the N-mixture models and GLMMs, we found no effect of shrub stem density until we separated shrubs by species in our post hoc exploratory models. This discrepancy highlighted the difference between studying spatial associations, as we did in our models, versus a comparison of pre- and posttreatment effects. Our findings confirmed those of Franca et al. (2016), who found that using a “space-for-time” approach, where different sites are compared, underestimated the consequences of landscape change and habitat disturbance, compared to the approach of measuring effects pre- and posttreatment. Despite the advantages of measuring skink abundance pre- and posttreatment, there were drawbacks in our study: (1) the habitat management was applied intentionally rather than randomly, so it was applied where it was needed the most—in sites that had the biggest shrub problem; and (2) our measures of observation frequency did not account for variation in detection. This latter issue could be viewed as a major problem because fewer shrubs could mean simply that skinks were more detectable. However, variation in detection likely had little effect on our results given three aspects of our analysis: (1) the near-equivalent results between the N-mixture and GLMM models, (2) the low variation in N-mixture detection estimates (5–11%), and (3) the magnitude of the changes we observed (97–188% increases in detections posttreatment).

To summarize a complex situation for shrubs, our results support management of shrubs in savanna and rock outcrop systems, but they also support treating sumacs differently from other nonnative shrub species. Sumacs are potentially problematic species that can interfere with a functioning savanna or prairie rock outcrop habitat. Nevertheless, smooth sumac also is native to these habitats, and at low densities or in scattered patches may benefit P. fasciatus by providing moisture, cover for prey, cover from aerial predators, and temperature regulation. Other native shrubs may also provide these benefits, but data were insufficient to isolate other individual species. However, our results reinforce justifications for reducing buckthorn as much as possible or practical.

Detection varied with cover-object temperature in most models, both N-mixture and GLMMs. Cover-object temperature was negatively associated with detection probably because of how it affected the thermoregulatory responses of skinks. Cover objects may have transferred heat to skinks more than other natural cover, both because of their material (metal or wood) and because they were low to the ground.

We found that using two different modeling techniques reinforced our confidence in these results and provided insights into assumptions of the models as they relate to P. fasciatus ecology. The similar abundance coefficients between the two types of models showed that separation of detection effects from abundance effects, as inherent in N-mixture models, was not a requirement for studying skink–habitat associations (Barker et al. 2018). If there had been a large source of varying detection that we failed to include in these models (e.g., a habitat variable that affected whether skinks used the cover objects), we would expect greater differences in the coefficients between the two model types.

Another key difference between the two techniques was that the N-mixture models we used required the assumption of a closed population. This assumption is often violated because of the demographic processes inherent in any wildlife population. In our case, a clear violation of this assumption was that skinks were hatching while we were surveying. We could have dropped young-of-the-year from the data to reduce violation of the closed-population assumption, but we felt it was important to include them because they were young-of-the-year and their habitat use is an important consideration for management. Because of their low numbers, we were unable to model young-of-the-year separately from adults. Our results were only negligibly affected by violation of the closed-population assumption because results from the N-mixture models were very similar to those of the GLMMs.

For this landscape, our results and management recommendations should be interpreted to co-occur with bedrock outcrops that include deep crevices—a landscape feature with highly restricted distribution in Minnesota and considered important for P. fasciatus. Further, our recommendations reflect current climate conditions for the Upper Minnesota River Valley. Similar to what Feltham and Nocera (2022) predicted for the Southern Canadian Shield population, P. fasciatus may shift from open to more closed-canopy sites adjacent to rock outcrops if average annual temperatures continue to increase because of climate change. Our recommendations for biotic factors complement research supporting the importance of abiotic factors for P. fasciatus (Feltham and Nocera 2022).

Targets or thresholds are important for identifying when, where, and what type of habitat management actions are needed to support conservation. Performance-measure approaches are highly variable, with no single standard for setting management targets. Minimally, targets should be measurable, evidence based, and sufficiently specific, yet practical for managers to apply (Hilton et al. 2022; Hilton and Cook 2022).

The targets we recommend for improving P. fasciatus habitat in the Minnesota River Valley and similar landscapes are given in the following. We derived targets from the N-mixture model predictions (Table 7), except where otherwise noted.

We added the upper bound of 25% based on the post hoc models showing the nonlinear relationship between canopy closure and skink abundance. Skink abundance tended to peak at around 25% (Fig. 4), and we thought it would be easy for managers to visualize in the field.

Derived directly from N-mixture model prediction (Table 7).

As a categorical variable, potential management targets were restricted to values distinguishing factor levels. Herbaceous cover >75% had the strongest effect on skink abundance compared to other levels (Tables 3 and 4). Because it was a consistently negative effect, we recommend less than 75%.

There is a need to be conservative with sumac because of its propensity to overtake habitat, but its potential benefits must be recognized. We lowered the prediction of 22 stems per 10 m² (Table 7) to 15 stems per 10 m², and treated it as an upper bound. We chose two stems per 10 m² as the lower bound because it was the lowest value predicting our target abundance.

N-mixture model predictions supported a buckthorn density of zero. We added an upper bound of five stems per 10 m² as a concession to the challenges of attempting to eradicate the species entirely.

We developed these targets for this and similar northern landscapes with the understanding that in most cases a manager will not need to promote woody and herbaceous plant growth, but rather will need to control them in this disturbance-dependent habitat. However, we emphasize that woody and herbaceous plants are important and a manager’s objective should not be to eradicate them. A mosaic of trees, herbaceous vegetation, down woody debris, and open, rocky areas is needed to ensure that skinks are able to access habitat that allows them to secure their various needs including cover, prey, and thermoregulation. Management that promotes this habitat patchiness, while increasing the area within a site that meets these management targets, would be a good approach for creating or improving P. fasciatus habitat.

Land cover is dynamic because of factors such as vegetative succession, invasive species, climate change, and anthropogenic habitat alteration. Therefore, we recommend monitoring at intervals to assess skink response and management needs. Assessments should also document site management itself such as treatment areas, types, and dates to help inform adaptive management. Hecnar and Brazeau (2017) noted benefits of management practices that enhance and maintain open habitats through mechanical clearing or prescribed burns and augmenting microhabitat by the placement of suitable woody debris. We anticipate that management of woody encroachment (especially buckthorn) will be ongoing and complicated because of rugged terrain and other factors. Although fire is beneficial for maintaining open habitats, fires over consecutive years may have a negative effect on skink populations by killing individuals or eliminating down woody debris (Klemens 1993). However, prescribed fires often do not carry well on rock outcrops. Terrain may preclude the use of heavy equipment for mechanical removal, and hand cutting is labor intensive. Goat grazing may be a viable option in some circumstances, but still likely must be combined with other management such as cutting and targeted, judicious herbicide application. Better understanding of soil and trophic effects of buckthorn also may inform updated best management practices for P. fasciatus and their associated ecosystems. Comparing relative decay rates and other properties of woody species may help prioritize which to retain as downed debris for skinks (Boddy 1983; Mun and Prewitt 2011). Future research should include improved assessments of skink travel distances and their use of adjacent grasslands and closed-canopy forest to inform habitat management and protection further. Radiotelemetry may reveal additional insights regarding P. fasciatus habitat use (Brazeau and Hecnar 2018; Feltham 2020).

Sustaining native biodiversity should include conserving suitable habitat conditions for underrepresented, at-risk species, as well as implementing effectiveness monitoring as part of an adaptive management strategy. Reptiles are underrepresented in monitoring, conservation planning, and implementation (Tingley et al. 2016). Our study of P. fasciatus microhabitat requirements helps to address this gap by quantifying key thresholds for this disjunct population. This will inform subsequent development of best management practices for P. fasciatus habitat in the Upper Minnesota River Valley, which may be applicable in similar landscapes.

This project was funded by the U.S. Fish and Wildlife Service’s State Wildlife Grants MN T-5-R-5 F14AF01118 and MN T-5-R-6 F16AF00287; Minnesota Environment and Natural Resources Trust Fund (ENRTF) as recommended by the Legislative‐Citizen Commission on Minnesota Resources (LCCMR); the Minnesota Department of Natural Resources, Division of Ecological and Water Resources, Nongame Wildlife Fund, Division of Parks and Trails, Scientific and Natural Areas Program; and the Minnesota River Collaborative. We thank E. Geary, T. Grealish, L. Groebner, K. Leuenberger, N. Schiltz, A. Yen, and other agency staff for their assistance during our study. We thank B. Bolduan and other land managers and landowners with whom we worked to identify potential sites and obtain site access. We also thank T. Arnold, University of Minnesota, for statistical and analytical suggestions, and C. Jennelle, B. Bolduan, and anonymous reviewers for comments on earlier drafts of this manuscript. We conducted all research in compliance with State of Minnesota laws, including under Department of Natural Resources Scientific and Natural Areas Program Special Permits for Snake and Skink Surveying and Monitoring Numbers 2015-22R, 2016-22R, and 2017-23R.

Supplemental material associated with this article can be found online at https://doi.org/10.1655/Herpetologica-D-24-00016.S1.