Predicting and mapping Plethodontid salamander abundance using LiDAR-derived terrain and vegetation characteristics

Aim of study: Use LiDAR-derived vegetation and terrain characteristics to develop abundance and occupancy predictions for two terrestrial salamander species, Plethodon glutinosus and P. kentucki, and map abundance to identify vegetation and terrain characteristics affecting their distribution. Area of study: The 1,550-ha Clemons Fork watershed, part of the University of Kentucky’s Robinson Forest in southeastern Kentucky, USA. Material and methods: We quantified the abundance of salamanders using 45 field transects, which were visited three times, placed across varying soil moisture and canopy cover conditions. We created several LiDAR-derived vegetation and terrain layers and used these layers as covariates in zero-inflated Poisson models to predict salamander abundance. Model output was used to map abundance for each species across the study area. Main results: From the184 salamanders observed, 63 and 99 were identified as P. glutinosus and P. kentucki, respectively. LiDAR-derived vegetation height variation and flow accumulation were best predictors of P. glutinosus abundance while canopy cover predicted better the abundance of P. kentucki. Plethodon glutinosus was predicted to be more abundant in sites under dense, closed-canopy cover near streams (2.9 individuals per m2) while P. kentucki was predicted to be found across the study sites except in areas with no vegetation (0.58 individuals per m2). Research highlights: Although models estimates are within the range of values reported by other studies, we envision their application to map abundance across the landscape to help understand vegetation and terrain characteristics influencing salamander distribution and aid future sampling and management efforts.

To effectively manage salamanders, it is essential to understand how habitat characteristics affect their distribution and abundance across the landscape. Plethodontid salamanders are lungless and rely solely on cutaneous respiration. This physiological constraint limits surface activity to cool and moist conditions (Jørgensen, 1997). Furthermore, most species have small home ranges and low vagility (Liebgold et al., 2011). Collectively, most studies have suggested that the distribution and abundances of salamanders is primarily influenced by fine-scale habitat conditions related to temperature and moisture. For example, Peterman & Semlitsch (2013), using a 3-m resolution National Elevation Dataset in a GIS and temperature dataloggers to derive spatial covariates, found that salamander abundance was best predicted by indices related to cooler temperatures and higher moisture, including denser canopy cover, especially in ravine habitats and areas on the landscape with low solar exposure and high topographic wetness. Thus, knowledge of fine-scale habitat attributes is essential to salamander management (Stauffer, 2002).
Light detection and ranging (LiDAR) technology can potentially provide detailed vegetation and terrain information needed for accurately describing fine-scale habitat important for salamanders. LiDAR data consist of three-dimensional point clouds with sub-meter positional accuracy. These data can be processed to segment points into ground and vegetation points. Ground points are used to create high-resolution digital elevation models that represent terrain surfaces and vegetation points can be analyzed to develop vegetation metrics that describe vegetation structure (Reutebuch et al., 2005). LiDAR-derived vegetation metrics such as canopy height, canopy cover, canopy heterogeneity, and understory density have been used to characterize the habitat of a number of birds and bats (i.e., Goetz et al., 2010;Tattoni et al., 2012;Eldergard et al., 2014;Jung et al., 2012;Müller et al., 2013), nonflying mammals (Zhao et al., 2012;Coops et al., 2010;Nelson et al., 2005), invertebrates (Müller et al., 2014;Vierling et al., 2008;Müller & Brandl, 2009), lizards (Sillero & Gonçalves-Seco, 2014) and turtles (Yamamoto et al., 2012;Long et al., 2011). However, there are no studies using LiDAR-derived vegetation metrics to characterize salamander habitats and construct models predicting local salamander population size. One probable reason for this in the Appalachian Mountain region is the complex vegetation conditions consists of dense, close-canopy deciduous forest with numerous tree species and highly dissected terrain conditions (Müller et al., 2014;Hamraz et al., 2017).
In this study, we used LiDAR data to describe fine-scale vegetation (i.e., canopy cover vegetation height, and vegetation height standard deviation) and terrain characteristics (i.e., slope's exposure to light, and water flow accumulation) in the Appalachian mountain region of eastern Kentucky. We used LiDAR-derived vegetation and terrain metrics to develop predictive abundance and presence models of two similar species of terrestrial salamander, the Slimy Salamander (Plethodon glutinosus) and the Cumberland Plateau Salamander (Plethodon kentucki). Lastly, we used these models to map salamander abundance across the study area and gain an understanding of the LiDAR-derived vegetation and terrain characteristics influencing their abundance.

Study area
Research was conducted at The University of Kentucky's Robinson Forest (RF), located in the rugged eastern section of the Cumberland Plateau region of southeastern Kentucky in Breathitt, Perry, and Knott counties. Due to access restrictions, the study area was limited to the 1,550-ha Clemons Fork watershed (Fig. 1). The terrain across the study area, and RF in general, is characterized by a branching drainage pattern, creating narrow ridges with sandstone and siltstone rock formations, curving valleys and benched slopes. The slopes are dissected with many intermittent streams (Carpenter & Rumsey, 1976) and are moderately steep ranging from 10 to over 100 % facings predominately northwest and southeast, and elevations ranging from 252 to 503 m above sea level. Vegetation is composed of a diverse contiguous mixed mesophytic forest made up of approximately 80 tree species with northern red oak (Quercus rubra L.), white oak (Q. alba L.), yellow-poplar (Liriodendron tulipifera L.), American beech (Fagus grandifolia E.), eastern hemlock (Tsuga canadensis (L.) Carr.) and sugar maple (Acer saccharum Marshall) as dominant and co-dominant species. Understory species include eastern redbud (Cercis canadensis L.), flowering dogwood (Cornus florida L.), spicebush (Lindera benzoin L.), pawpaw (Asimina triloba (L.) Dunal), umbrella magnolia (Magnolia tripetala L.), and bigleaf magnolia (M. macrophylla Michx.) (Carpenter & Rumsey, 1976;Overstreet, 1984). Average canopy cover across RF is about 93 % with small opening scattered throughout. Most areas exceed 97 % canopy cover, but recently harvested areas have an average cover as low as 63 %. After being extensively logged in the 1920's, RF is considered second growth forest ranging from 80 to 100 years old and is now protected from commercial logging and mining activities, typical of the area. Seventeen species of salamander, most of which belong to the 3 Mapping Plethodontid salamander abundance using LiDAR-derived information family Plethodontidae (Schneider, 2010;Petranka, 1998), are found at RF. Some of the most abundant terrestrial salamander species are P. glutinosus and P. kentucki, which are the focus of this study. These lungless salamanders prefer cool moist habitats and are most active on the ground surface at night after rain events (Petranka, 1999).

LiDAR derived data
A high-density (~ 40 pt m -2 ) LiDAR dataset was acquired in the summer of 2013 during leaf-on season for collecting detailed vegetation information across RF. The parameters of the LiDAR system and flight are presented in Table 1. A set of five LiDAR-derived variables were created to predict and map salamander abundance across the study area. LiDAR ground points were used to create a 0.6 m resolution digital elevation model (DEM) with average as the cell assignment method and natural neighbor as the void fill method using the "LAS dataset to Raster" tool in ArcMap 10.2. Terrain variables included two raster layers based on the DEM: hillshade (HS) and flow accumulation (FA), which were created using the "Spatial Analyst" tool also in ArcMap 10.2. The HS layer was used as a proxy for direct sun exposure, which considered the average daily position of the sun when field data was collected (175° azimuth and 70° altitude). The FA la-yer represents the number of upslope cells theoretically flowing onto a given cell, which provides an indication of relative humidity. These two layers were resampled into a courser 30.5 m resolution using the average cell value to encompass entire field transects into single cells and consider a more appropriate cell size to meaningfully capture site variations across the study area.
LiDAR vegetation points were normalized using the DEM to calculate elevation above ground and used to create three vegetation variables: canopy cover (CC), vegetation height (VH), and vegetation height standard deviation (VHSD). The CC layer was calculated as the percentage of vegetation points above 5 m from ground level to the total points for all 0.6-m cells covering the study area. The 5 m threshold was selected to avoid considering LiDAR points representing ground vegetation, typically up to 3 m tall across the study area, and to be above the maximum elevation change error of the DEM found to be ~1.5 m . Each cell was considered covered and given a value of one if the percentage was greater than 50 % and not covered and given a value of zero otherwise. For data consistency, the CC layer was resampled into coarser 30.5 m resolution using the average cell values. The VH layer was calculated as height of the tallest LiDAR vegetation point inside the 0.6 m cell size and then resampled to the courser 30.5 m resolution using the maximum cell value. Lastly, the

Forest Systems
August 2020 • Volume 29 • Issue 2 • e005 VHSD layer was calculated as the standard deviation of the vegetation height of 0.6 m cells within coarser 30.5 m cells. This layer was created to represent the variability of vegetation heights, which tends to be higher in recently harvested areas and lower in areas with fully closed canopies.

Sampling design
To quantify the abundance of salamanders, we used a stratified sampling where 45 field transects were surveyed across varying soil moisture and canopy cover conditions throughout the study area (Fig. 2). We created an integrated soil moisture index (SMI) layer and used the CC layer to identify the location of these transects. The GIS-based SMI layer (Iverson et al., 1997) was developed to determine soil moisture, which considers terrain slope, direct sun exposure (hillshade), ground curvature, and soil water holding capacity data from the United States Geological Survey. The SMI layer had a 10-m resolution with each cell representing relative soil moisture across the study area and was resampled into a coarser 30.5-m resolution using the average cell values. The SMI layer was classified into high, medium, and low soil moisture classes selecting threshold values resulting in an equal amount of area in each class (516.7 ha). The CC layer was also classified into three classes: low (0-50 % covered) medium (50-75 % covered), and high (75-100 % covered) canopy cover. Lastly, five transects were randomly located in each soil moisture/canopy cover combination using the center point of the raster cells as the transect location, which were not allowed in areas within 5 m from existing roads and streams to avoid their effects on detected salamanders.

In field data collection
The location of the mid-point of transects was determined using a Trimble Juno SB GPS handheld unit with a 6m-precision. The 30.5 m transects were laid out along the contour line and flags were placed at the ends and mid-point to establish a clear line of sight along their length. We used a visual encounter survey to collect salamander count data. Transects were surveyed at nighttime on days following rain events during May -June of 2014. They were surveyed using a headlamp to search inside a 1 m swath along either side of its length. Encountered salamanders were captured, placed in Ziploc bags, and left at the same place where they were found to minimize disturbance to the site. After transects were searched, caught salamanders were examined and species was recorded.
All transects were sampled three times, as required for presence and abundance modeling to account for imperfect detection (Royle 2004;MacKenzie et al., 2002). Transect locations were grouped so several could be accessed in one night, then groups were randomly surveyed with no transects being revisited within three days of the last survey. We also collected six sample-specific detection variables at each transect during each visit, namely: litter depth (cm), Julian date, wind speed, barometric pressure (mmHg), air temperature (°C) using a Kestrel 2500 weather meter, and soil moisture (%) using an Extech MO750 soil moisture probe.
A 15.2-m buffer was placed around each transect, covering an area of 0.17 ha, to maintain consistency with the resolution of the covariates. This buffer area was used to extract a single value for the covariates associated with each transect for model development purposes.

Data analysis
Before model development, we used a Pearson's correlation matrix to examine the relationship among the five LiDAR-derived variables as well as the SMI. Due to the large amount of transect surveys with no salamander observations, we used zero-inflated abundance models that accounted for imperfect detection of individuals. Specifically, we used the statistical function RunZIA (Wenger, 2007;Wenger & Freeman, 2008) in the R software (R Development Core Team, 2008). The RunZIA function, based on N-mixture models (Royle 2004, Royle et al. 2005) and zero-inflated binomial occupancy model (MacKenzie et al. 2002), simultaneously estimates occupancy (or presence), abundance and incomplete A total of 32 predictive models were then developed to estimate salamander abundance, considering all 31 unique combinations of these five predictive variables (HS, FA, CC, VH, VHSD) plus the SMI. In addition to these six abundance variables, all models included Julian date and the days since last precipitation event squared as detection variables because seasonal and weather variables have shown to greatly influence desiccation rates and affect surface activity, and thus detection probability (Petranka, 1998;Peterman & Semlitsch, 2014). We ranked models based on the small sample size Akaike information criterion (AICc) and the evidence ratio (ER), and the weighted Akiake criterion (W AICc) was used to determine the relative performance of the best model. Lastly, models were run separately for both species to determine if there were differences in site preference.
In order to estimate abundance, separate model parameter estimates (β) for presence and abundance were output for the best model for each species using the Run-ZIA. The parameter estimates were calculated based on a binary function for presence (Eq. 2) and an exponential function for abundance (Eq. 3). Using ArcMap the parameter estimates were applied to the raster file of each Figure 2. Location of field plots within the study area. First two letters in the abbreviated plot categories indicates level of canopy cover (CH = canopy cover high level, CM = canopy cover medium level, CL = canopy cover low level) and the second two letters indicates level of soil moisture (MH = soil moisture high level, MM = soil moisture medium level, ML = soil moisture low level).

Forest Systems
August 2020 • Volume 29 • Issue 2 • e005 covariate (X) using equations 2 and 3 to create a presence raster file and preliminary abundance raster filer for each species. Then equation 1 was applied to the resulting presence and abundance raster files to map the estimated abundance, given their presence.

Results
The Pearson's correlation matrix among LiDAR-derived covariates showed a high correlation between CC and VH (r = 0.89) ( Table 2). Despite this correlation, it is important to consider both variable as VH is required to, in general, distinguish between recently harvested areas and older forests, both often with relatively high CC. Another high correlation was found between SMI and HS (r = -0.81), which is expected as HS is used to compute SMI, for which reason no model included both covariates. Correlation among other pairs covariates were relatively low presenting values lower than 0.27, except for FA and SMI (0.42).
A total of 184 salamanders were observed from the three visits to each transect. P. glutinosus and P. kentucki were the most abundant, with a total of 63 and 99 observations, respectively. There were no salamanders observed in 63% of the visits (85 out of 135), which justified the use of a zero-inflated model to determine salamander abundance (Table 3).
Using AICc to evaluate models, we found that the best-supported model for both species was one with a zero-inflated Poisson (ZIP) assumption on abundance. Results from developing the 32 predictive models of P. glutinosus abundance show that the top-ranked model retained VHSD and FA (Table 4, which shows only the top 16 best-ranked models). This model has evidence of over 26 times of, and about 96 % more likely to perform better than the second-ranked model (as shown by the evidence ratio and weighted Akiake criterion, respectively), which also contains HS. Also, VHSD and FA are present in most of the other high-ranked models, which likely indicate their influence in P. glutinosus abundance.
When examining the complete top-ranked model, none of the 95 % confidence intervals of the abundance parameter estimates overlap zero, which indicates the respective covariates are likely to be important in the model (Table 5). The coefficient estimates show an inverse effect of VHSD on abundance while FA has a direct relationship, but for the presence portion of the model both covariates have a direct effect. A plot showing predicted abundance per transect (sampled area of 61 m 2 ) against ranges of VHSD and FA values found across the study area, helps visualize the effect of these two predictors on P. glutinosus abundance (Fig. 3a). For example, on places with homogeneous vegetation heights (i.e., closed canopy, dense forests), predicted salamander abundance ranges from 1.08 to 12.01 individuals per m 2 based on FA values. Similarly, on places with FA near zero (near ridgetops), salamander abundance ranges from 0.04 to 1.08 individuals per m 2 for varying VHSD values. The abundance model predicts a maximum of 12 salamanders per m 2 at site with high FA and low VHSD. These conditions are likely to occur at sites under dense, closed canopy cover near streams. However, when the abundance and presence models are combined (Fig. 3b), it shows more realistic predictions with a maximum of 2.9 salamanders per m 2 .
Results from running the predictive models of P. kentucki abundance show a less clear best-fit model (Table  4). However, the top-ranked model, which only contains CC as the abundance and presence variable, is about 48 % more likely to perform better than other models. Moreover, CC is also contained in six of the top seven models. In the top-ranked model, the positive CC parameter estimate for abundance indicates a direct relationship and  (Fig. 4). When comparing both salamander species, results show that different LiDAR-derived covariates have a significant effect on their abundance. For example, VHSD is retained in several high-ranked models (AICc < 216) for P. glutinous abundance, while CC is retained in most lowranked models (AICc > 219) ( Table 3). The opposite case can be observed for P. kentucki where CC and VHSD are retained in high-ranked (AICc < 279) and low-ranked models (AICc > 279), respectively. When applying the topranked abundance / presence model to map the abundance across the study area, it can be observed that P. glutinosus is predicted to be present in relatively high numbers near streams while not occupying ridgetops, which offers further evidence of the effect of FA on abundance (Fig. 5). However, P. kentucki is predicted to be more abundant throughout except for roads surfaces and recently harvested areas with low CC value, closely resembling the CC special distribution (Fig. 6).

Discussion
We demonstrate the utility of using LiDAR-derived terrain and vegetation information needed to estimate salamander presence and abundance. LiDAR-derived VHSD and FA, and CC were found to be the best predictors for P. glutinosus and P. kentucki abundances, respectively. This is an important finding due to the need to accurately describe fine-scale habitat important for salamanders over large areas, which it would be difficult with traditional, courser remotely sensed data (Peterman & Semlitsch 2013). LiDAR data acquisition is becoming more affordable and datasets are becoming available at the regional scale. For example, LiDAR datasets are now available for large parts of the states of Kentucky, West Virginia, Virginia, North Carolina, and Tennessee, covering most of the Appalachian region where salamanders are an important part of the ecosystem functioning.
LiDAR point density has been shown to affect DEM accuracy (Balsa- Barreiro & Lerma, 2014;Hodgson & Bresnahan, 2004). However, the LiDAR data used in this study is a high-density dataset with 40 pts m -2 with enough points reaching the ground to create a DEM. A previous study  quantified the DEM accuracy across the study area using our high-density (40 pts m -2 ) dataset collected during leaf-on season and a low-density (1.5 pts m -2 ) dataset collected during leaf-off. They found the mean elevation change error to vary from 23 cm to 146 cm based on terrain slope and ruggedness, and most importantly they found no significant Table 3. Field data collection summary of the 45 transects grouped by canopy cover and soil moisture index category. First two letters in the abbreviated plot categories indicates level of canopy cover (CH = canopy cover high level, CM = canopy cover medium level, CL = canopy cover low level) and the second two letters indicates level of soil moisture (MH = soil moisture high level, MM = soil moisture medium level, ML = soil moisture low level) differences between datasets. Considering that available regional LiDAR datasets have point densities between 1.5 and 40 pts m -2 , similar DEM accuracies can be expected. In addition to point density, cell size can also affect DEM accuracies. However, resampling the ori-ginal 0.6 m resolution of the LiDAR-derived layers, to the courser 30.5 m resolution, to match the size of field transects smoothed and averaged cell values and the associated error. Although out of the scope of this study, future study might focus on evaluating the accuracy of Mapping Plethodontid salamander abundance using LiDAR-derived information vegetation height and canopy cover as a function of point density and cell size across the study area as well as the Appalachian region.
Model development results indicate that both salamander species have different habitat preferences. P. glutinosus was predicted to be more abundant in sites under dense, close-canopy cover near streams. This corresponds with several other studies reporting preference for near humid sites (i.e., Marvin, 1996;Davidson, 1956;Grobman, 1944). P. kentucki was predicted to be found across the study sites except in sites with no vegetation. This also agrees with several studies mentioning rocky outcrops, downed logs, leaf litter and living roots systems as suitable habitat for this species, and which are found across the study area where canopy cover is high (Bowers, 2013;Pauly & Watson, 2005;Marvin, 1996). Abundance estimates also similar to those reported by other studies. For example, Burton & Likens (1975a) reported salamander abundance about one third of our estimates (approximately 0.25 salamanders per m 2 vs our combined average  The parsimonious nature of the developed models can facilitate its use as they include one or two LiDAR-derived covariates. Predictor covariates are in line with known phenomena of desiccation effecting salamander activity and abundance (Peterman & Semlitsch, 2014). Although, forest age also affects salamander abundance and presence (Petranka, 1999), it is difficult to determine in eastern deciduous forests, but basal area or diameter at breast height could be used as surrogate. We did not include such covariates due to the difficulty of retrieving individual tree information from LiDAR data in closed canopy deciduous forests (Hamraz et al., 2016;Koch et al., 2006).
Because we used canopy cover as a surrogate for desiccation, areas with low canopy cover should have been better represented in the random selection of transects. This was difficult to achieve because most of the study area is a considered second growth forest, with almost full canopy closure throughout which is mostly covered. Recently harvested areas were the only areas with medium and low canopy cover. Models for both species contained some form of a vegetation variable predicting lower abundance in those areas, which is likely related to increased desiccation from more direct sunlight via canopy openings. Alternative transect selection methods to select transect locations should be used to reduce the number of field observations with zero counts and thus improve model performance. The field data collection could also be improved by increasing the sample size and limiting the data collection to days closer to rain events to ensure sampling during time periods with more surface activity.
We present the first attempt to quantify salamander abundance using LiDAR-derived fine-scale vegetation and terrain information in the deciduous forest of the Appalachian mountain region of eastern Kentucky.  Mapping Plethodontid salamander abundance using LiDAR-derived information Variation in vegetation height and flow accumulation were important predictors of P. glutinosus abundance, and LiDAR-derived canopy cover was the only important predictor of P. kentucki abundance. Methods could be replicated by land and wildlife managers for different species of terrestrial plethodontid salamanders to identify vegetation and terrain characteristics affecting their distribution across the landscape and to model their relative abundance. The presence and abundance models developed can reasonably predict salamander abundance providing estimates within the range of values reported by other studies. However, we recommend their use to estimate relative abundance, instead of estimating population size or biomass. A straightforward application of these models is to map abundance across the landscape to help understand vegetation and terrain characteristics influencing salamander distribution and assist with future more rigorous sampling and management efforts.