Environmental gradients in old‐growth Appalachian forest predict fine‐scale distribution, co‐occurrence, and density of woodland salamanders

Abstract Woodland salamanders are among the most abundant vertebrate animals in temperate deciduous forests of eastern North America. Because of their abundance, woodland salamanders are responsible for the transformation of nutrients and translocation of energy between highly disparate levels of trophic organization: detrital food webs and high‐order predators. However, the spatial extent of woodland salamanders’ role in the ecosystem is likely contingent upon the distribution of their biomass throughout the forest. We sought to determine if natural environmental gradients influence the fine‐scale distribution and density of Southern Ravine Salamanders (Plethodon richmondi) and Cumberland Plateau Salamanders (P. kentucki). We addressed this objective by constructing occupancy, co‐occurrence, and abundance models from temporally replicated surveys within an old‐growth forest in the Cumberland Plateau region of Kentucky. We found that Plethodon richmondi had a more restricted fine‐scale distribution than P. kentucki (mean occupancy probability [ψ¯^] = 0.737) and exhibited variable density, from <250 to >1000 individuals per hectare, associated with increased soil moisture and reduced solar exposure due to slope face. While more ubiquitously distributed (ψ¯^ = 0.95), P. kentucki density varied from <400 to >1,000 individuals per hectare and was inversely related to increased solar exposure from canopy disturbance and landscape convexity. Our data suggest co‐occurrence patterns of P. richmondi and P. kentucki are influenced primarily by abiotic conditions within the forest, and that populations likely occur independently and without evidence of biotic interaction. Given the critical role that woodland salamanders play in the maintenance of forest health, regions that support large populations of woodland salamanders, such as those highlighted in this study—mesic forest stands on north‐to‐east facing slopes with dense canopy and abundant natural cover, may provide enhanced ecosystem services and support the stability of the total forest.


| INTRODUC TI ON
Analyzing the distribution and abundance of species along environmental gradients yields invaluable information about their niche requirements (Costa, Wolfe, Shepard, Caldwell, & Vitt, 2008), population dynamics (Peterman & Semlitsch, 2013), and biotic interactions (Maestre et al., 2010) and can even inform decisions about the management and restoration of landscapes for species conservation (Peterson, 2006). In unaltered landscapes, the distribution of species is a function of natural environmental gradients, which include abiotic factors (e.g., surface temperature, moisture, topographic relief, water and soil chemistry, and solar radiation) and biotic factors (e.g., vegetative structure and the presence of predators, prey, and mates; Van der Putten, Macel, & Visser, 2010). Taxa likely to exhibit strong responses to such natural gradients are those with limited dispersal capabilities (Cushman, 2006), low reproductive success (Elton, 1958), and acute sensitivity to environmental conditions (Buckley & Jetz, 2007).
Despite this sensitivity, amphibians represent a major component of biomass in aquatic (Gibbons et al., 2006), terrestrial (Burton & Likens, 1975b;Petranka & Murray, 2001), and riparian (Peterman, Crawford, & Semlitsch, 2008) ecosystems. Because the life history of many amphibians involves movement between and among aquatic and terrestrial ecosystems (Regester, Whiles, & Taylor, 2006), they are responsible for the transformation (Burton & Likens, 1975a) and translocation (Capps, Berven, & Tiegs, 2014;Luhring, DeLong, & Semlitsch, 2017) of substantial quantities of energy throughout the landscape. However, the role of energy transformation is not unique to biphasic organisms. Terrestrial woodland salamanders (Caudata: Plethodontidae: Plethodon), which lack aquatic larval stages (i.e., have direct development), are among the most abundant vertebrate animals in eastern deciduous forests of North America (Petranka & Murray, 2001;, reaching densities between 0.73 and 18.46 individuals per m 2 Semlitsch et al., 2014). They also act as predators of detrital food webs (Best & Welsh, 2014;Davic & Welsh, 2004;Hutton, Price, & Richter, 2017) and represent a prey resource for a wealth of vertebrate and invertebrate predators (for a taxonomic review of Plethodon predators, see Semlitsch et al., 2014). As such, woodland salamanders are hypothesized to serve as a key energetic intermediary between highly disparate levels of trophic organization in terrestrial ecosystems (detrital communities and high-order vertebrate predators; Burton & Likens, 1975b) and exert a significant, top-down, regulatory force upon detrital food webs, leaf litter decomposition, and organic material retention (Burton & Likens, 1975a;Hairston, 1987). Therefore, woodland salamanders may significantly influence the direction and magnitude of energy flow through ecosystems (Davic & Welsh, 2004). Wyman (1998) suggested that, through predation of detrital food webs, woodland salamanders (Plethodon cinereus, eastern redbacked salamander) can indirectly reduce leaf litter processing rates, aiding in the retention of organic carbon in forests. However, additional studies have found that the significance, strength, and direction of top-down effects on leaf litter decomposition and detrital communities is subject to variation (Best & Welsh, 2014;Hocking & Babbitt, 2014;Homyack, Sucre, Haas, & Fox, 2010;Walton, 2005;Walton & Steckler, 2005;Walton, Tsatiris, & Rivera-Sostre, 2006).
Recent evidence suggests that variation in the effects of woodland salamanders on forest floor dynamics is likely correlated with spatiotemporal variability in environmental conditions (Walton, 2013) and the abundance of salamander predators (Hickerson, Anthony, & Walton, 2017). Therefore, the nature of woodland salamanders' role in terrestrial ecosystem nutrient cycling is likely contingent upon the spatial distribution of their biomass within the ecosystem (Hickerson et al., 2017;Semlitsch et al., 2014), which is influenced by spatial patterns in environmental conditions and resource availability (Milanovich & Peterman, 2016;Peterman & Semlitsch, 2013;Walton, 2013).
Numerous studies have found the distribution of woodland salamanders to be influenced chiefly by terrestrial ecosystem features such as soil moisture (Jaeger, 1971a;Peterman & Semlitsch, 2013;Wyman, 1988) (Gibbs, 1998;Peterman & Semlitsch, 2013). Furthermore, the presence of heterospecifics has been found to influence microhabitat usage (Farallo & Miles, 2016;Keen, 1982), distribution (Hairston, 1950;Jaeger, 1970Jaeger, , 1971a, and abundance (Hairston, 1951) of individual species. Thus, the species-specific contribution of woodland salamanders to terrestrial ecosystem processes may be modified through population-level effects of interspecific competition. Due to the diversity and endemism of woodland salamanders, particularly in Appalachian forests where their diversity is greatest (Dodd, 2004), community structure varies dramatically across physiographic regions. Therefore, community interactions are likely geographically nuanced and not easily generalizable from any single region. within an old-growth forest in the Cumberland Plateau region of Appalachia. In old-growth forests of this region, variation in environmental conditions of the forest floor (e.g., soil moisture, availability of woody debris, and solar exposure) is largely influenced by canopy dynamics (Runkle, 1982). Tree mortality supplies woody debris to the forest floor (Harmon et al., 1986) and provides habitat for woodland salamanders (McKenny et al., 2006;Petranka, Brannon, Hopey, & Smith, 1994); however, resultant canopy gaps increase solar exposure, accelerating evapotranspiration. The size and persistence of canopy gaps represent a natural disturbance regime, which greatly modifies local environmental conditions (Schaetzl, Johnson, Burns, & Small, 1989;Scharenbroch & Bockheim, 2007). The objectives of this study were to determine if environmental gradients associated with the natural disturbance regime of an Appalachian old-growth forest influence the fine-scale distribution and density of P. richmondi and P. kentucki. Furthermore, this study sought to determine if patterns of salamander co-occurrence vary along natural environmental gradients, and if those patterns are modified behaviorally through interspecific competition and territoriality. These objectives are addressed by constructing hierarchical models, which incorporate imperfect detection from temporally replicated surveys within an old-growth forest. With no history of timber harvest, the old-growth forest at LCW has experienced no substantial anthropogenic disturbance with the exception of understory livestock grazing, which ended in the 1950s (Martin, 1975). Canopy disturbances in LCW are primarily stochastic (Davis, Chapman, Wu, & McEwan, 2015), and therefore, the distribution of canopy gaps is predicted to be uniform and resultant from endogenous processes. Of the three tracts of old-growth forest at LCW, one tract, "Shop Hollow," currently experiences little disturbance from human recreation (only guided hiking on a narrow,

| Amphibian sampling
This study relied upon visual encounter surveys (VES) to detect species, and therefore, all observations resulted from hand captures during standardized searching. MacGregor pers. comm.), and therefore, this sample period was chosen as the most representative of true patterns in occurrence and abundance. VES were conducted along a linear 3-m × 32-m transect (96 m 2 ), which intersected the center of each 800-m 2 circular sample plot at a randomly chosen bearing between 0° and 180°. The bearing of each transect was also randomized during each sequential visit, making the likelihood of sampling the same microhabitat at a given sample plot negligible. In LCW, woodland salamanders are found primarily by searching under natural cover (coarse woody debris [CWD] and rocks) on the forest floor. During VES, all CWD and rocky cover within the 96-m 2 transects were flipped, and microhabitats beneath were examined for the presence of salamanders before replacing cover items to their exact position.

| Site covariates
During each visit, soil moisture was measured with a Pro Check moisture probe (Decagon Devices, Inc.) at five equidistant points along the transect within each sample plot, and therefore, estimates were obtained by averaging across transect (N = 5) and visit (N = 4).
Quantification of forest canopy openness was achieved using hemispherical canopy photography (Baldwin, Calhoun, & DeMaynadier, 2006;Frazer, Lertzman, & Trofymow, 1997;Herbert, 1987). A single photograph of canopy structure was captured prior to leaf off with a 24-megapixel digital single-lens reflex camera (Nikon D7100) on automatic settings, fitted with a 180° lens (Nikon AF DX Fisheye-Nikkor 10.5 mm f/2.8G ED; Nikon Instruments, Melville, NY, USA). Percent canopy openness was calculated by converting images into binary color (black pixels = closed canopy, white pixels = open canopy) using a binarization algorithm provided by the Auto Threshold Plugin for ImageJ software (Abramoff, Magalhaes, & Ram, 2004;Rasband, 2014), and then calculating the percent of white pixels in each frame.
A 1.11-m 2 digital elevation model was used to derive the following layers: aspect, slope, Topographic Position Index, and direct solar radiation. Aspect was scaled into a linear variable ranging from 0 (xeric, southwest-facing slopes) to 2 (mesic, northeast-facing slopes) using the Beers transformation (Beers, Dress, & Wensel, 1966;O'Donnell, Thompson, & Semlitsch, 2015). Topographic Position Index (TPI) is a measure of landscape convexity which was calculated by comparing the slope position of individual sample plots relative to a 150-m 2 surrounding landscape area using a neighborhood function (Guisan, Weiss, & Weiss, 1999). During the calculation of TPI, a suite of additional neighborhood scales was considered, beginning with a circular area of 50 m 2 and increasing incrementally by 50-1,500 m 2 .
The most appropriate TPI scale was selected by correlating each TPI calculation with plot averages of raw salamander counts and identifying the scale with the highest correlation coefficient. Direct Solar Radiation-a component of the total solar radiation-represents the quantity of solar radiation remaining after a fraction is absorbed by the atmosphere (diffuse solar radiation) or reflected off of the earth's surface (reflected solar radiation). Normalized Difference Vegetation Index (NDVI) is a measure of vegetative cover (range: −1.0 [barren] to 1.0 [heavily vegetated]) and was derived using imagery from the 2016 National Agriculture Imagery Program. All data were gathered with ArcGIS 10.3 (ESRI, 2011). See Table 1 for a description of all sampling and site covariates.

| Sampling covariates
The quantity of fallen coarse woody debris larger than 20 cm in diameter (Muller & Liu, 1991) Table 1 for a description of all sampling covariates.
Nichols, Hines, Knutson, & Franklin, ) were used to estimate the probability that a species occupied a given site (ψ), while N-mixture models (Royle, 2004) were used to estimate species true population size (λ). Fitting occupancy and N-mixture models followed a stepwise procedure: (a) models were constructed to estimate detection parameters, p, by holding the state parameters, occupancy and abundance, constant; (b) the model-averaged estimated effect size (̂) of each detection covariate was calculated using multi-model inference (Burnham & Anderson, 2002;Mazerolle, 2006)   . Most single-season occupancy models, including the model used in this study, estimate the "conditional capture probability" (p ), defined as the probability of capture, given the individual is present (capture probability|availability). For these terms, availability is defined as 1-(temporary emigration). Single-season N-mixture models estimate a form of detection which combines a term for the ability of the observer to capture an individual that is present (conditional capture probability) with a term for the individual's availability for capture (expressed as: availability × conditional capture probability) and is thus referred to as an "effective detection probability" (p λ ). By exploiting the relationship between effective detection probability and conditional capture probability, estimates of population capture availability and temporary emigration (probability an animal is alive, but unavailable for capture) can be obtained from single-season models mathematically.
Two-species single-season occupancy models (MacKenzie, Bailey, & Nichols, 2004) were used to estimate the probability that P. richmondi occupies a site or sites wherein P. kentucki is known to be present. Under the null hypothesis, the pattern and frequency of co-occurrence does not vary across environmental gradients. This hypothesis was tested by comparing a null model of co-occurrence, wherein the pattern in which species co-occur at sites is unrelated to environmental conditions (essentially random), to models of cooccurrence, which predict co-occurrence patterns relating to environmental gradients. Using the co-occurrence probability (ψ AB ), a "Species Interaction Factor", or φ, can also be obtained (MacKenzie et al., 2004;Richmond, Hines, & Beissinger, 2010). For species A and B, φ is defined as: where ψ A and ψ B are the occupancy probabilities of species A and B, and ψ AB represents the co-occurrence probability of species A and B. Under the null hypothesis, φ = 1, species populations exist independently and the pattern and frequency of species cooccurrence are assumed to be random. If φ > 1, species co-occur more frequently than expected from chance; likewise, φ < 1 indicates species co-occur less frequently than chance. Using the same modeling procedure as with single-species occupancy models, two-species candidate models were fitted within the maximum-likelihood framework provided by program PRESENCE (v. 11.7) under a ψ Ba -parameterization (Richmond et al., 2010) and ranked using AIC. Predictions were obtained from the highestranking models.

| Detection, Availability, and Temporary Emigration
The conditional capture probabilities of P. richmondi and P. kentucki were moderately low (p = 0.36 and 0.24), while effective detection probabilities were much lower (p = 0.06 and 0.05; Table 2).  Table 2). The conditional capture probability of P. kentucki increased sharply with rising quantities of fallen CWD, while the effective detection probability increased only gradually (Figure 1).
The probability that salamanders were alive, on the soil surface, and available for surveys (Availability probability) was low due to frequent vertical temporary emigration into the soil by both P.
richmondi and P. kentucki ( Table 2). The probability of temporary emigration by P. richmondi varied slightly with time of day, moderately increasing, linearly, from morning until evening. The temporary emigration probability of P. kentucki did not vary substantially, but exhibited a hump-shaped relationship with the abundance of coarse woody debris (Figure 1). Table 2), although large confidence intervals provide considerable uncertainty in our assessment. Percent soil moisture ("MST"), NDVI ("VEG"), and canopy openness ("CAN") were all important covariates in estimating occupancy of P. richmondi (Figure 2). Like P. richmondi, P. kentucki occupancy was also correlated with percent soil moisture and NDVI (Figure 3), but the directions of the covariates' effects were heterogeneous (Figure 2). The remaining covariates included in models of occupancy produced heterogeneous effects and were therefore not considered to be reliable predictors of woodland salamander distributions in LCW.

| Co-occurrence
The overall probability of P. richmondi co-occurring with P. kentucki, ̂r ic|ken was 0.72 (95% CI: 0.53, 0.86). Models of co-occurrence featuring covariates that represent environmental gradients were better at predicting patterns of co-occurrence (cumulative Akaike model weight [Σω ij ] = 0.971) than null models (Σω ij = 0.029). Co-occurrence probabilities were positively influenced by percent soil moisture and NDVI (Figure 4)  Percent soil moisture ("MST") and Beers-transformed aspect ("ASP") were the most important covariates when estimating density of P. richmondi (Figures 2 and 5). Plethodon richmondi density exhibited marked, positive curvilinear responses to percent soil moisture and aspect. Plethodon kentucki density was influenced most by Topographic Position Index ("TPI") and percent canopy openness ("CAN"; Figure 2). The density of P. kentucki exhibited gradually dampened negative responses to both Topographic Position Index and percent canopy openness, with inflated upper limits ( Figure 5).

| D ISCUSS I ON
We found that natural environmental gradients created by dynamic ecosystem processes inherent in old-growth forest influence the fine-scale distribution, co-occurrence, and density of P. richmondi and P. kentucki. Species-specific responses to gradients of soil moisture and temperature, solar exposure from canopy structure, and slope position likely reflect physiological restraints associated with desiccation vulnerability and thermal avoidance. Although patterns in co-occurrence of P. richmondi and P. kentucki do vary along gradients of canopy density and soil moisture, little evidence was found to support the hypothesis that populations of woodland salamanders experience interspecific competition.

| Plethodon richmondi
Plethodon richmondi density in LCW was positively related to forest soil moisture and reduced solar radiation due to slope face, and their fine-scale distribution (i.e., occupancy) was restricted to mesic for- Difference Vegetation Index, "ASP" = Beers-transformed aspect, "TPI" = Topographic Position Index, "ELV" = Digital Elevation Model, "RAD" = direct solar radiation, "MST" = volumetric soil moisture, "CAN" = percent canopy openness. See Table 1 for a description of covariates. Covariates with ̂ values centered at zero were not estimated due to non-convergent models or model instability and are therefore represented as having "zero" effect sizes pest is presumed to continue. It follows that through alterations to canopy characteristics, A. tsugae, and other invasive pests in LCW (e.g., emerald ash borer, Agrilus planipennis) could negatively impact Plethodon salamanders, which lack the vagility to evacuate habitats that have undergone dramatic transformation (Welsh & Droege, 2001). Further research into the mechanisms responsible for canopy loss in LCW may provide a more meaningful interpretation of P. richmondi occupancy and density dynamics. Future surveys and analyses should incorporate data pertaining to tree age, diameter, canopy density, and prevalence of pest-related damage.

| Plethodon kentucki
Density of P. kentucki was negatively impacted by the presence of canopy disturbance on exposed slope faces. Furthermore, canopy disturbance influenced the density of P. kentucki with a greater magnitude than soil moisture and aspect-gradients which both bol- Factors affecting local density of P. kentucki did not necessarily affect their distribution. Specifically, canopy disturbance and landscape convexity (exposure) were found to be key determinants of density at a given site, but did not necessarily influence the likelihood of that site being occupied. These results suggest that within LCW P. kentucki respond with greater consequence to environmental processes which govern population size (i.e., productivity, recruitment) than those which perhaps govern their occurrence (i.e., colonization, extinction).
It is possible that aspects of environmental gradients assessed in this study were not important in explaining population dynamics of P.
kentucki or P. richmondi due to the coarse scale with which they were evaluated. For instance, environmental variables and salamander counts were aggregated to the extent of the sampling area (800 m 2 ), a scale which could obscure relevant information about the relationship of salamanders with micro-scale variation in environmental conditions.

| Co-occurrence
The degree of overlap in the fine-scale distributions of P. richmondi and P. kentucki within LCW corresponded strongly with natural environmental gradients. The probability of P. richmondi and P. kentucki co-occurring in a given forest stand at LCW was positively  (Marvin, 1998), it is possible that the spatial scale of this study was too coarse to quantify such fine-scale interactions. Another potential explanation of the observed patterns of cooccurrence may be related to mating behavior of P. kentucki. Marvin (1998) found that populations of P. kentucki in this region exhibit territoriality associated with mate pairing. In southeast Kentucky, the breeding period of P. kentucki begins late June to mid-August and lasts until mid-to-late October (Baecher pers. obs., Marvin & Hutchison, 1996). Although unrelated to interspecific competition, it is possible that territoriality associated with P. kentucki breeding behavior was not observed during the timeframe of this study (15 October 2016 to 13 November 2016).

| Detection, Availability, and Temporary Emigration
Unless all individuals in a population are available for capture during a survey (availability = 1), it is important to distinguish between conditional capture probability (probability of capturing an animal given availability = 1) and effective detection probability (probability of capturing an animal given availability ≤1). Given that Plethodon are known to migrate between surface and subsurface refugia frequently (Bailey, Simons, & Pollock, 2004), their availability-the probability of an individual being alive and present on the soil surface during a survey-should be much <1 (availability = 1 -[temporary emigration]; , and therefore, estimates of effective detection probability should be much less than that of the conditional capture probability. Corroborating the assertion that surface inactivity confounds studies of Plethodon salamanders , this study showed that more than half of the individuals in populations of P. richmondi and P.
kentucki had emigrated into subterranean refugia and were unavailable for surveying.

| CON CLUS IONS
This study found that the pattern of distribution and the abundance of woodland salamanders throughout the landscape can be nonrandom. Given that the nature of woodland salamanders' effects on forest floor dynamics (e.g., detrital food webs, organic material retention) can change due to variation in environmental conditions (Walton, 2005(Walton, , 2013, it is likely that the spatial extent of woodland salamander's influence on the ecosystem is nonrandom and varies dramatically across natural environmental gradients . Thus, the role that woodland salamanders play in the maintenance of forest health, biodiversity, and ecosystem services (Davic & Welsh, 2004) is likely contingent upon the inherent inhabitability of the system. Therefore, regions within a forest that support large populations of woodland salamanders, such as those highlighted in this study-mesic forest stands on north-to-east facing slopes with dense canopy-may provide enhanced ecosystem services and support stability in the total forest ecosystem (Davic & Welsh, 2004). This study took place in a stable F I G U R E 4 Co-occurrence of Southern Ravine Salamander (Plethodon richmondi) and Cumberland Plateau Salamander (P. kentucki) in relation to NDVI (Normalized Difference Vegetation Index), a measure of vegetative cover, and percent soil moisture in an old-growth forest at Lilley Cornett Woods Appalachian Ecological Research Station, Letcher County, Kentucky, USA in Fall 2016 old-growth forest (Martin, 1975), virtually undisturbed by human activity (with the exception of light recreation from guided hiking). Because woodland salamanders in this study exhibited such marked responses to natural disturbances associated with forested ecosystems (e.g., isolated canopy perforation and soil desiccation due to solar exposure), population-level responses to nonnatural disturbances (e.g., timber harvest and residential/commercial development) are hypothesized to be much more substantial. Committee (IACUC protocol # 05-2015).

AUTH O R S' CO NTR I B UTI O N S
JAB and SCR conceived the ideas and designed the methodology; JAB collected and analyzed the data; JAB and SCR wrote the manuscript.

DATA ACCE SS I B I LIT Y
Materials to reproduce analyses are freely available for download by visiting the following repositories.