Relationships among white-tailed deer density, harvest, and landscape metrics in TN, USA

Landscape and harvest indices are frequently used to represent white-tailed deer (Odocoileus virginianus) density. However, the relationship between deer density and specific landscape indices is unclear. Harvest is another metric often used to estimate deer density. Our objective was to model the relationship among deer density, landscape metrics, and harvest density of deer in TN, USA. We estimated deer density across 11 regions in 2011 using distance sampling techniques. We developed 18 a priori models to assess relationships among deer density, harvest density, and landscape metrics. Estimates of deer density ranged from 1.85 to 19.99 deer/km2. Deer density was best predicted by harvest density and harvest density + percent woody area. However, harvest density was the only important variable in predicting deer density (Σωi = 0.700). Results of this study emphasize the significance of harvest density in deer management. While the importance of harvest as a management tool for deer is likely to increase as landscapes are fragmented and urbanized, specific management guidelines should be based upon deer densities and landscape metrics when they are important.


Introduction
Human-related causes of increased ungulate populations in Europe and eastern North America include a decrease in the number of hunters, climate change, increased number of people in urban areas, management for increased population sizes in rural areas, and land use changes (Maillard et al. 2010). The number of hunters has declined over time as people have developed other social and cultural interests (Enck et al. 2000). Only 5% of the population in the USA (US Department of the Interior 2016) and 0.5% of the population in Europe hunts (Reimoser and Reimoser 2016). Climate change has produced warmer winters, and more animals survive as a result (Maillard et al. 2010;Davis et al. 2016). A trend to live in urban areas has left much of the rural land undeveloped thereby providing more ungulate habitat (Apollonio et al. 2010). Rural inhabitants manage for increased ungulate populations to provide for economic opportunities given the declining number in people and associated economics (Apollonio et al. 2010). Land use, which is related to the aforementioned causes, can provide a landscape level understanding to match the distributions of ungulates. Landscape varies spatially and temporally, and little work has been conducted (Roseberry and Woolf 1998;Lovely et al. 2013) at the landscape scale to understand these relationships. Previous research indicated the importance of landscape metrics, such as amount of forest edge for deer (Plante et al. Electronic supplementary material The online version of this article (https://doi.org/10.1007/s10344-019-1353-8) contains supplementary material, which is available to authorized users. 2004) and habitat fragmentation and diversity for white-tailed deer (Odocoileus virginianus; Quinn et al. 2013), both of which are influenced by land use, in predicting population densities.
Harvest, or bag, is the most common metric used to track hunted ungulate populations (Roseberry and Woolf 1991;Apollonio et al. 2010). However, most studies have demonstrated this relationship at a local or site-level and not a landscape scale (McCullough 1979;Hagen et al. 2018). Managing for maximum sustained yield and incorporating harvest of females results in a reduction in the total population given a sufficient number of females is harvested (Downing 1981;Demarais and Zaiglin 1988).
Not unlike many big game program objectives, TN's objective for white-tailed deer (Odocoileus virginianus; hereafter "deer") populations across the state is to balance the biological and social carrying capacities (Kelly et al. 2019). Human population growth is increasing in TN, as elsewhere, resulting in urban sprawl and a reduction in habitat for wildlife (Robinson et al. 2005). Concurrent with habitat loss is an expectation by sportsmen to maintain a harvestable population (biological carrying capacity) of deer while minimizing loss of plant diversity, crop damage, and deer-vehicle collisions (social carrying capacity). As landscapes change over time, a model predicting population density of deer to compare with the biological and social carrying capacities is needed to manage this species. Therefore, our objective was to model the relationship among deer density, harvest density, and landscape metrics in TN, USA. We hypothesized landscape characteristics, specifically density of woody edge, percent forested, and harvest served as predictors of deer density at the landscape level.

Study area
TN lies in the southeastern region of the USA (latitude 34°59′ N to 36°41′ N; longitude 81°39′ W to 90°19′ W) covering approximately 109,000 km 2 . TN may be divided into 6 major physiographic regions (west to east), i.e., East Gulf Coastal Plain, Highland Rim, Central Basin, Cumberland Plateau, Ridge and Valley, and Unaka Mountains. Counties were assigned to the physiographic region in which they occurred. Counties that contained > 1 physiographic region were assigned to the region with the greatest area in that county. The resulting sub-regions (Griffith et al. 2012;TN Climatological Service 2013;Tennessee Federal GIS Users Group 2013) were East Gulf Coastal Plain west (EGCP-W) and east (EGCP-E), Highland Rim north (HR-N), central (HR-C), and south (HR-S), Central Basin north (CB-N), central (CB-C), and south (CB-S), Cumberland Plateau (CP), Ridge and Valley (RV), and Unaka Mountains-Ridge and Valley (UMRV; Fig. 1). By making use of the sub-regions based on physiographic regions allowed us to sample areas with similar habitats.

Sampling and field procedures
We used a multivariate regression approach to assess the relationship among residual (i.e., after harvest) deer density (response variable), harvest density, and landscape metrics (predictive variables). Deer density was used as the response variable instead of an index of density because it was the direct result of vital rates, and indices rely on critical and unlikely assumptions regarding detection probabilities (Anderson 2001). We used a random sampling design for density estimation using transects as replicates ) within each sub-region. We randomly selected 3 secondary roads in each county to be used as transects (n = 285). Transects were at least 0.5 km apart to maintain independence. We sampled transects beginning 0.5 h after sunset and ended sampling by 0200 h; transects were not sampled when raining. Each transect was sampled once between 1 February and 4 April 2011. One observer stood in the back of a pickup truck driven 8-16 km/h and surveyed one side of the transect using a ProTech Thermal-Eye 250D thermal imager to locate deer. Two additional people were in the back of the truck; one person operated a spotlight and aided with taking measurements, and one recorded data. When a group of deer was located, we measured the perpendicular distance from the road to the center of the group or the location where the group was first sighted with a rangefinder. We recorded group size, time, and the vehicle location using Universal Transverse Mercator coordinates with a GPS unit for each group observed.
We used Program Distance 6.0 to estimate deer density for each sub-region . We set the sampling fraction to 0.5 because we sampled only one side of each transect. We modeled densities based on sample sizes greater than the minimum recommended by Buckland et al. (2001). Landscapes within sub-regions were different across the state, and different models were likely better suited on estimating density in some sub-regions than others. Therefore, we used the half-normal, uniform, and hazard rate key functions in conjunction with cosine, simple polynomial, hermite polynomial, and no series expansions. We right-truncated 5-10% of the data to improve model fit . We calculated the Akaike Information Criterion value corrected for small sample size (AIC c ; Akaike 1973), ΔAIC c value, a density estimate, upper and lower 95% confidence levels for density estimates, χ 2 goodness of fit statistic, coefficient of variation, and probability of detection for each model. We selected the best model based on the minimum AIC c value.
Statewide harvest seasons were divided into 3 general types, i.e., archery, muzzleloader, and gun. Harvest season, across all 3 types, ran from the third Saturday in September through the first Sunday in January. Two young sportsman hunts were open across 1 weekend before and 1 weekend following the gun season, wherein no other hunting was allowed; harvest from this hunt was traditionally minimal. Three harvest management units having different bag limits were used to provide management relative to the perceived deer density. The western two-thirds of the state were in the liberal management unit and, in addition to a 3-buck limit per season, this management unit allowed 3 does per day to be taken; the other 2 units allowed 4-10 does per season depending upon the unit (Yoest et al. 2012). Each hunter harvested an average of 1.9 deer. The harvest regime, therefore, was effectively the same across the state.
Successful deer hunters in TN are required to check their harvest with the TN Wildlife Resources Agency (TWRA). During the 2010-2011 hunting season, animals were checked using 1 of 2 approved methods: (1) check stations (> 800 locations throughout the state) or (2) internet check. Harvest data collected included date of harvest, location of kill (e.g., county or Wildlife Management Area), type of weapon used, sex of harvested deer, and number of antler points if the animal was male.
We developed a list of landscape variables; we thought would have a strong influence on deer density and were found to be important in other studies (e.g., Bobek et al. 1984;Gaudette and Stauffer 1988;Roseberry and Woolf 1998;Anderson et al. 2001;Plante et al. 2004;Long et al. 2005;Pettorelli et al. 2007;Munro et al. 2012 2014), we calculated the 2010 human population size in rural areas for each county by subtracting the urban population from the total. We calculated percent public land by dividing the area of public land by the total land area of each county and multiplying by 100. Because deer are thought to respond to landscape metrics beyond the extent of the home range (Kie et al. 2002), matching landscape scale density estimates to landscape scale metrics was appropriate.
We determined the total area and total perimeter to calculate perimeter-area ratio representative of an edge metric. Number of patches and median patch size for agricultural areas, pastures and grasslands, urban, deciduous woodlands, coniferous woodlands, and mixed woodlands for each county (ESRI 2014) were used to calculate the ratio of edge between open (agricultural, pasture, and grassland combined) and urban areas, open and woodland (deciduous, coniferous, and mixed combined) areas, and urban and woodland areas; these metrics were also representative of edge at a finer scale than simple perimeter-area ratio. We calculated percent urban land cover and percent rural land cover (open and woodland combined) by dividing each respective land cover class by the total land area of each county and multiplying by 100. Average slope for each county was calculated based on 30 m digital elevation models in ArcGIS.

Statistical analyses
We reduced this predictive variable set by removing one of each correlated variable pair (p > 0.5). We performed linear regression using SAS PROC REG (SAS Institute 2011) to model the relationship among estimated deer density (dependent variable) and landscape metrics, deer harvest (Table 1),  Table 2). We calculated AIC c values for each of 18 models and used the minimum AIC c value to select the best model. Models with ΔAIC c ≤ 2.0 were considered competing models. We calculated the weight (ω i ) of each competing model based on ΔAIC c values. Because there was no clearly superior model (i.e., ω i > 0.9), we model-averaged parameter estimates (Burnham and Anderson 2002).
Harvest density was the primary variable in competing models predicting deer density in 2011 (Table 4). Given the great potential that the system was not operating in a linear fashion but instead in a non-linear and non-parametric fashion, we conducted modeling of the data and found the linear (i.e., parametric) models were ranked, via AIC, much higher than the non-linear (i.e., spline) models (results of this analysis provided in Online resource 1). Moreover, harvest was yet the best predictor of density. The non-linear, non-parametric models ranked far below the linear models. The best model was > 2.5 times more than likely to accurately predict deer density than the second best model (harvest density + percent woody area). The top model also accounted for 50% of the variation predicting density (Fig. 2). Averaged across all models, harvest density was the only important variable explaining deer density (Σω i = 0.700, Table 5) and did not include zero in the 95% confidence intervals.

Discussion
Our results indicated that residual deer density was best predicted by density of deer harvested. Other studies reporting harvest as a predictive variable for deer density have been disparate. Hansen et al. (1986) determined deer population size in IL was positively correlated with harvest, while Pettorelli et al. (2007) found no correlation between deer density and the number of deer harvested per day in Québec. Roseberry and Woolf (1998) speculated statewide pre-hunt deer density was regulated by harvest, even though density was considered low (e.g., 4-5 deer/km 2 ). Five of our subregions were estimated to have densities < 5 deer/km 2 , and these occurred in eastern and western portions of the state (the EGCP-W, RV, UMRV, CP, and HRS sub-regions had densities of 1.85, 2.54, 2.53, 3.30, and 4.13 deer/km 2 , respectively). Eastern regions of TN are known to contain lower   . 2 Regression between harvest density (deer/km 2 ) and deer density (deer/km 2 ) in TN, USA, 2011. Solid line represents the predicted relationship and the dashed lines represent 95% confidence intervals quality deer habitat, and deer in this part of the state are dependent upon unpredictable mast production (Johnson et al. 1995;Ryan et al. 2004) which would explain the lower observed deer densities. Lack of association between density and any of the selected landscape metrics was surprising. The implication would be that deer respond to harvest more numerically than to habitat. Some areas in TN, those in central TN in particular, experienced much human population growth over the last two decades, and with that growth fragmentation of the landscape occurred. Western and eastern portions of the state experienced little human population growth by comparison. Harvest in some sub-regions of western and eastern TN demonstrated up to 49% harvest of the population during this study. One possible explanation for our lack of finding any landscape metrics as important is that harvest of the population was additive to the natural mortality that occurs through disease (e.g., epizootic hemorrhagic disease) and predation, thereby holding the population low enough such that habitat factors were not quite sufficient to be statistically identified. This suggests the population was below carrying capacity and just above the point of maximum growth. In the central portion of the state where human population was growing and fragmentation was increasing, carrying capacity may be such that the deer population was responding numerically to the improved habitat but had not reached a point where habitat is a factor. These sub-regions demonstrated low proportions of the populations being harvested. In short, we think two possible processes are in play. First, fragmentation has been occurring in the central portion of the state at such a rate that the deer population has not yet reached a level to where the habitat factors are substantially affecting population growth. These sub-regions do not have a high percentage of the population harvested. The eastern and western areas do have a high proportion of the population harvested (approaching 50% in some cases) and these are likely kept well below K, such that the harvest is more important than landscape metrics in explaining the results.
Because deer-not just white-tailed (Alverson et al. 1988;Anderson et al. 2001), but also other species (i.e., black-tailed [Odocoileus hemionus columbianus; Kremsater and Bunnell 1992], sika [Cervus nippon ;Takatsuki 1989], roe [Hansson 1994;Tufto et al. 1996])-are characteristically edge species, it would be expected that the amount of open-woody edge in an area would better predict deer density (Waller and Alverson 1997). We found, however, open-woody edge was secondary to harvest density and did not contribute to the model predicting deer density.
Previous research indicates deer prefer forested edges (Aulak and Babińska-Werka 1990;Marques et al. 2001;Torres et al. 2012). Home range selection is strongly based on landscape metrics, including edge (Beier and McCullough 1990;Kie et al. 2002;Fulbright and Ortega 2006). Moreover, landscape metrics, including edge, have been tied to reproductive fitness (McLoughlin et al. 2007;Miyashita et al. 2008). As edge has increased through landscape, fragmentation deer have likely selected smaller home ranges because of more available forage and responded by increasing population size; this is supported by the highest densities being associated with areas of increasing population growth and associated fragmentation.
Based upon our results, deer populations in TN are still harvest-dependent but not influenced by landscape use. Biologists and managers will be challenged to maintain the harvest tool as the number of sportsmen and sportswomen decline in the USA (US Department of Interior et al. 2016). Models incorporating landscape metrics to track ungulate populations in both Europe and the USA have been many (Foster et al. 1997;Radeloff et al. 1999;Acevedo et al. 2011), and inclusion of landscape metrics (deCalesta and Stout 1997;Putman et al. 2011) in addition to harvest density in models to predict deer population density would be less costly than estimating population density on a regular basis. Harvest is often the only common metric collected for big game by wildlife agencies (Unsworth et al. 2002;Apollonio et al. 2010). However, when landscape metrics are not related to density, as found in this work, harvest remains an important and crucial tool (Milner et al. 2006) to manage population densities. Acknowledgments We also thank the numerous wildlife officers, biologists and technicians from the Tennessee Wildlife Resources Agency who assisted with data collection. Editorial assistance was provided by B. E. Flock. We also wish to thank anonymous reviewers, the Editor and Associate Editor for their recommendations and suggestions to improve this work.
Funding information Financial support for the study was provided by the Tennessee Wildlife Resources Agency using Pittman-Robertson Federal Aid in Wildlife Restoration funds.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.