POPULATION IMPACTS AND DISEASE PATTERNS OF HEMORRHAGIC DISEASE IN
WHITE-TAILED DEER: UNDERSTANDING THE ROLE OF AN EMERGENT DISEASE
By
Sonja A. Christensen

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Fisheries and Wildlife – Doctor of Philosophy
2018

ABSTRACT
POPULATION IMPACTS AND DISEASE PATTERNS OF HEMORRHAGIC DISEASE IN
WHITE-TAILED DEER: UNDERSTANDING THE ROLE OF AN EMERGENT DISEASE
By
Sonja A. Christensen
The hemorrhagic diseases (HD) epizootic hemorrhagic disease (EHD) and bluetongue
(BTV) are the most significant source of viral-related mortality of deer in the United States. The
emergence of HD on the landscape is affecting the population ecology of white-tailed deer
(Odocoileus virginianus) in previously unaffected regions of North America. The forces driving
apparent increases in HD are poorly understood, particularly where the disease has recently been
most severe in northern latitudes. Specifically, Michigan has seen an increase in EHD related
deer mortality since 2006, and 2012 was the largest outbreak of EHD in Michigan history. My
research seeks to address current uncertainties for deer populations with regard to impacts from
emergent HD. In this dissertation I attempt to investigate the local distribution and populationlevel impacts on deer in areas affected by EHD, develop and apply population estimation
techniques appropriate for localized population change assessment, and evaluate landscape-level
drivers of HD across an epidemiological gradient.
I estimated the annual population abundance of deer immediately after a recent EHD
outbreak over a 4-year period and quantified the fine-scale spatial extent of EHD impacts
associated with riparian habitats. I assessed these differences in deer abundance for sites affected
and unaffected by EHD in Michigan. I determined if populations increased over time in response
to local population impacts from EHD. I then validated estimates using an independent aerial
survey and a novel application of N-mixture models. I evaluated the N-mixture model against

distance sampling methods and tested the sensitivity of N-mixture model estimates to variation
in spatial unit size. At a regional scale, I evaluated the role of drought severity in both space and
time on changes in HD reports to determine if drought severity explains changing patterns of HD
presence.
Important findings from my research include 1) abundance estimates in the affected area
were significantly lower along transects near the river, reflecting EHD mortality associated with
wetlands; 2) the opposite was true in the unaffected site; 3) estimates from both distance
sampling and N-mixture modeling approaches were similar and aerial surveys using N-mixture
binomial models are a practical tool for estimating abundance of a free-ranging wildlife
population; 4) the selected spatial unit size in the N-mixture model had a significant effect on
abundance and detection probability estimates, and their associated variance; 5) drought severity
was a significant predictor of HD presence and the significance of this relationship was
dependent on latitude; 6) the effect of drought was reduced or non-existent in most southern
states, where the disease is endemic.
My findings suggest that EHD outbreaks will have population impacts near riparian
corridors and those populations can increase post-outbreak. Further, N-mixture models may
provide an improved approach to estimate populations from aerial surveys in the future, although
managers and researchers should evaluate the spatial unit used for these models prior to survey
implementation. While drought severity does increase the probability that HD will be detected
and reported at a county level, my research points to the underlying role of acquired herd
immunity across the endemic-to-emergent disease gradient

ACKNOWLEDGEMENTS
The culmination of this dissertation represents a series of opportunities and challenges,
and people who were willing to take a chance. To these people, and to all those who have
supported me in this endeavor, thank you. To Dr. Paul Krausman, who steadfastly encouraged
my pursuit of a doctoral degree at a critical time, I am forever grateful.
I would like to extend my utmost appreciation to my co-advisors Dr. William Porter and
Dr. David Williams. Their time, effort, professionalism and positive mentoring approach resulted
in a rigorous yet enjoyable work environment. I thank my advisory committee, Dr. Jean Tsao and
Dr. Elise Zipkin, for their thoughtful guidance, collegiality, and review efforts. I extend sincere
appreciation to Dr. Brent Rudolph, who served as the research lead for this work through the
Michigan Department of Natural Resources and provided valuable feedback, insight, and
friendship throughout my time at Michigan State University.
Many collaborations and opportunities were made possible during my doctoral work,
resulting in a much more fulfilling learning experience than I could have ever anticipated from
an intellectual pursuit alone. I would like to recognize Dr. David Stallknecht, Dr. Mark Ruder,
and Dr. John Fischer for graciously welcoming me as a collaborator and affording me the
opportunity to work on data collected by the Southeastern Cooperative Wildlife Disease Study. I
would like to thank Rebecca Humphries for her mentorship and invitation to work with her on
the National Fish and Wildlife Health Initiative. Much appreciate goes to the Michigan members
of the Boone and Crockett Club and their partners group, the many employees of the Michigan
Department of Natural Resources, and the board members of the Glassen foundation, for
providing financial resources and time in support of my research.

iv

It was during my doctoral program that I met and married my husband, David Ortega.
David has been a true gift in my life and I am eternally grateful for his love, support, patience,
and academic advice during the many ups and downs of graduate school. I greatly appreciate the
unwavering support and encouragement from my family, specially my parents Jean and Kit, my
brother Michael and my sister-in-law Cathy. I owe so much of my work ethic and personal
character to them.
Lastly, I am thankful for the many friendships that I have made during my time at
Michigan State University. It has been a privilege to work with my fellow graduate students in
the Department of Fisheries and Wildlife. Special recognition and thanks goes to the students in
the Boone and Crockett Quantitative Wildlife Center who have provided countless laughs,
insightful discussions, constructive feedback, and constant encouragement.

v

TABLE OF CONTENTS

LIST OF TABLES ....................................................................................................................... viii
LIST OF FIGURES ....................................................................................................................... ix
INTRODUCTION ...........................................................................................................................1
LITERATURE CITED ....................................................................................................................6
CHAPTER 1: EPIZOOTIC HEMORRHAGIC DISEASE AFFECTS ABUNDANCE AND
DISTRIBUTION OF FREE-RANGING DEER POPULATIONS ....................................9
INTRODUCTION ..................................................................................................................... 9
STUDY AREA ........................................................................................................................ 13
METHODS .............................................................................................................................. 15
Data collection and study design ........................................................................................ 15
Distance sampling analysis................................................................................................. 18
RESULTS ................................................................................................................................ 21
Data collection .................................................................................................................... 21
Model selection ................................................................................................................... 21
Deer density and abundance estimates ............................................................................... 24
DISCUSSION .......................................................................................................................... 28
APPENDICES .............................................................................................................................. 33
APPENDIX A: SUPPLEMENTARY EHD SURVEILLANCE MAP OF DATA REPORTED
IN 2012 IN MICHIGAN. MAP PRODUCED BY THE MICHIGAN DEPARTMENT OF
NATURAL RESOURCES. .............................................................................................. 34
APPENDIX B: SUPPLEMENTARY DETECTION PROBABILTY HISTOGRAMS FOR
DISTANCE SAMPLING ESTIMATES IN CHAPTER 1.. ............................................. 35
LITERATURE CITED ................................................................................................................. 36
CHAPTER 2: APPLICATION OF N-MIXTURE MODELS FOR AERIAL SURVEYS OF
WILDLIFE: A COMPARISON STUDY ......................................................................... 41
INTRODUCTION ................................................................................................................... 41
STUDY AREA ........................................................................................................................ 45
METHODS .............................................................................................................................. 46
Aerial survey design ............................................................................................................. 46
Distance sampling analysis.................................................................................................. 49
N-mixture model analysis .................................................................................................... 50
Sensitivity analysis of site size for N-mixture models .......................................................... 52
RESULTS ................................................................................................................................ 53
DISCUSSION .......................................................................................................................... 60
APPENDICES .............................................................................................................................. 65
APPENDIX A: R CODE FOR N-MIXTURE MODEL LISTED IN CHAPTER 2................ 66
APPENDIX B: SENSITIVITY ANALYSIS CODE (2KM) FOR CHAPTER 2. ................... 70
LITERATURE CITED ................................................................................................................. 73

vi

CHAPTER 3: THE ROLE OF DROUGHT AS A PREDICTOR OF EMERGENT
HEMORRHAGIC DISEASE IN THE EASTERN UNITED STATES ........................... 77
INTRODUCTION ................................................................................................................... 77
STUDY AREA ........................................................................................................................ 82
METHODS .............................................................................................................................. 82
Data Collection ................................................................................................................... 82
Data Analysis ...................................................................................................................... 88
RESULTS ................................................................................................................................ 91
DISCUSSION .......................................................................................................................... 97
MANAGEMENT IMPLICATIONS ..................................................................................... 100
APPENDICES ............................................................................................................................ 103
APPENDIX A: MODEL VALIDATION RESULTS FOR CHAPTER 3 ............................ 104
APPENDIX B: MODEL OUTPUT FOR BEST MODEL FOR CHAPTER 3 ..................... 105
APPENDIX C: SUPPLEMENTARY MODEL SELECTION DATA FOR CHAPTER 2 ... 107
LITERATURE CITED ............................................................................................................... 109
POSTFACE................................................................................................................................. 115
LITERATURE CITED ............................................................................................................... 118

vii

LIST OF TABLES

Table 1.1. Notation and description for a priori distance sampling models for estimating
abundance and density estimates of white-tailed deer in the Maple River and Cass River study
areas of Michigan, 2013-2016………………………………………………………….………..20
Table 1.2. Models selected for each stratum in each year and corresponding selection criteria and
data descriptors………………………………………………………………………….……….23
Table 1.3. Parameter estimates from best fitting models in each stratum, study area, and year for
2013-2016 in Michigan. Asterisks represent 95% confidence intervals from bootstrapping..…..25
Table 2.1. Distance sampling models selected in each year and corresponding data
descriptors……………………………………………………………………………………..…56
Table 2.2 Distance sampling parameter estimates and associated 95% confidence intervals from
best fitting models in each year for 2014 and 2016 in Michigan…...……………………….......56
Table 2.3 Summarized posterior distributions for N-mixture model parameter estimates, using
the 2 km spatial unit, in each year for 2014 and 2016 in Michigan………………………..........57
Table 3.1 Notation, description, and number of parameters providing the base structure for a
priori models for estimating the probability of HD presence in a county and year during 20002014. All models are conditional on the random effects structure selected initially using
AIC.……………………………………………………………………………………….…......90
Table 3.2 Average coefficient estimates for the fixed effect parameters for HD mortality
presence by county across 23 states in the eastern US, for 2000-2014.…………………….…...97
Table C.1 Full model selection table using AIC criteria for all model variations predicting the
probability of HD presence in 1830 counties within 23 states during 2000-2014…………..…107

viii

LIST OF FIGURES

Figure 1.1 The location of line transects used to estimate fine-scale deer abundance at varying
distances from a riparian corridor in two study areas in Michigan, from 2013-2016. The Maple
River area (red transects) had previously reported severe EHD outbreaks (case site) and the Cass
River area (black transects) had no record of past EHD outbreaks (control site)……………….15
Figure 1.2 Deer density and 95% confidence interval estimates (deer/sq. km) and for the Maple
River study area (EHD site) for 2013-2016 in Michigan………………………………………..26
Figure 1.3. Deer density and 95% confidence interval estimates (deer/sq. km) and for the Cass
River study area (control site) for 2014-2016 in Michigan……………………………………...27
Figure A.1 A map of Epizootic Hemorrhagic Disease occurrence in Michigan 2012…………..34
Figure B.1 Fitted detection probability histograms of deer observations at binned distance
intervals for the Cass River and Maple River study areas from 2014-2016 and 2013-2016,
respectively, in Michigan……....………………………………………………………………..35
Figure 2.1 Aerial survey transects bisecting habitat in the Maple River Study Area of Michigan,
USA, 2014-2016…………………………………………………………………………………46
Figure 2.2 A comparison of abundance estimates for white-tailed deer from distance sampling
and N-mixture models with associated 95% confidence intervals and 95% credible intervals,
respectively. Abundance represents deer within the Maple River study area of Michigan for 2014
and 2016…………………………………………………………………………………………55
Figure 2.3. Sensitivity of 2014 N-mixture models to spatial unit size for a) population abundance
(N) and b) associated standard deviation (SD) of abundance in Michigan……………………...58
Figure 2.4 Sensitivity of 2014 N-mixture models to spatial unit size for a) detection probability
(p) and b) associated standard deviation (SD) of detection probability…………………………58
Figure 2.5 Sensitivity of 2016 N-mixture models to spatial unit size for a) population abundance
(N) and b) associated standard deviation (SD) of abundance……...…………………………....59
Figure 2.6 Sensitivity of 2014 N-mixture models to spatial unit size for a) detection probability
(p) and b) associated standard deviation (SD) of detection probability in Michigan…………....59
Figure 3.1 Wetland cover from 2006 NLCD data across the 1,830 county study area………….86
Figure 3.2 Physiographic regions of the Eastern US. …………………………………………...87
Figure 3.3 The probability of HD mortality presence with associated standard error (shaded area)
in logit scale on the y-axis and cumulative drought severity on the x-axis. The relationships are
then grouped by latitude. As latitude increases the slope of the line increases………………….94
Figure 3.4 The change in the marginal effect of latitude on the coefficient of drought severity
with associated standard error (shaded area). This demonstrates the statistical significance of
these conditional coefficients and the relationship between the interacting model terms (drought
severity and latitude). With increasing latitude the magnitude increases for the coefficient of
ix

drought severity on HD mortality A superimposed histogram is included at the bottom of the
marginal effect plot to view the distribution of the data…………………………………………95
Figure 3.5 The probability of HD mortality presence with associated standard error (shaded area)
in logit scale given the percent wetland cover by county and the interaction with drought
severity…………………………………………………………………………………………...96

x

INTRODUCTION

The emergence of an epidemic disease on the landscape is affecting the population
ecology of white-tailed deer (Odocoileus virginianus) in previously unaffected regions of North
America. The hemorrhagic diseases epizootic hemorrhagic disease (EHD) and bluetongue (BTV;
a clinically indistinguishable hemorrhagic disease) are the most significant source of viral-related
mortality of deer in the United States (Shope et al. 1960, Thomas et al. 1974). Once the virus has
been transmitted to a naĂŻve deer and onset of disease signs are observed, hemorrhaging typically
leads to shock, unconsciousness, and death within 8 to 36 hours (Shope et al. 1960). Beginning
in the late 1990s, the probability of EHD occurrence began to increase significantly for deer in
the Northeast and Midwest regions of the United States (Stallknecht et al. 2015). In these
regions, population immunity has been historically low and disease outbreaks resulting in
mortality are often severe, potentially causing impacts on population abundance (Fischer et al.
1995, Gaydos et al. 2004, Howerth et al. 2001, Park et al. 2013). In epidemic regions, outbreak
intensity and frequency have continued to increase significantly in the last decade (Stallknecht et
al. 2015). This growing disease risk for deer populations is largely understudied and causing
concern for wildlife managers and the public. My research seeks to address current uncertainties
for deer populations with regard to impacts from emergent EHD.
Understanding the impact of epidemic and emerging infectious disease on wildlife
populations is crucial for conservation and management as disease outbreaks intensify and
spread spatially (Harvell et al. 2002, Ostfeld et al. 2005). Where EHD is epizootic, uncertainty
exists regarding the complex and scale-dependent factors driving the spatial distribution of EHD.
The localized spatial variation in deer mortality may be due to multiple factors including the

1

pathogenicity of the specific virus serotype associated with the outbreak, the ecology and
competence of the disease vector, and acquired or innate host immunity dynamics (Howerth et
al. 2001, Stallknecht et al. 2002, Gaydos et al. 2002a, Gaydos et al. 2002b, Gaydos et al. 2004,
Ruder et al. 2015). For example, the population ecology and distribution of the vector, various
midge species in the genus Culicoides, play a significant role in the epidemiology of
hemorrhagic diseases (Mellor et al. 2000, Wilson and Mellor 2008, Savini et al. 2011). The
development of larval Culicoides is restricted to wetland habitats and moist soils associated with
river or stream banks (Wilkening et al. 1985). Further, most descriptions of detected EHDrelated deer mortalities are associated with water bodies (Gaydos et al. 2004). Despite the logical
connection between an increased risk of EHD for deer and wetland habitats harboring vector
populations, previous research has not quantified the fine-scale spatial extent of EHD impacts
associated with riparian habitats.
Climatic, environmental, and land-cover factors further affect spatial processes of EHD
disease risk at regional scales and likely interact with the factors listed above (Xu et al. 2012,
Berry et al. 2013, Ruder et al. 2015). In a recent review of research gaps in the epidemiology of
hemorrhagic disease in North America, Ruder et al. (2015) discuss the increasing incidence and
poorly understood patterns of hemorrhagic disease in northern latitudes and called for research
addressing the ecology of EHD in changing environments and climates. Further, recent work by
Stallknecht and colleagues (2015) stated that the significance of increases in reported
hemorrhagic disease on deer population health is unknown but potential population impacts are
possible. Trends in increasing disease reports across large epidemic regions of the United States
warrant investigation of landscape-level factors that could be driving these patterns. Because

2

vector populations are sensitive to climate conditions in particular (Purse et al. 2005), changing
climate patterns related to soil moisture availability and temperature should be investigated.
When sudden and severe EHD mortality events in ungulate populations occur, deer
managers face considerable uncertainty in knowing the impact of EHD on the local deer
population. Predicting and collecting pre-epizootic data at a scale relevant to the population
impact of the disease is not possible, so assessments of direct population impacts on free-ranging
deer are rare. Previous research has attempted to assess population impacts from EHD by
conducting seroprevalence surveys (a measure of the level of pathogen in a population) to
measure the post-outbreak immune response of a deer population or by counting deer carcasses
suspected to have succumb to EHD and comparing these values to estimated background deer
population estimates (Fischer et al. 1995, Gaydos et al. 2004). While past studies have provided
some of the only values for free-ranging deer impacts, the background populations used were
often gross calculations from regional harvest-based estimates (Fischer et al. 1995, Gaydos et al.
2004). The ability to assess disease prevalence, mortality, and predicted recovery of deer
populations at appropriate scales for disease, rather than at management-unit scales, presents a
clear need for the application and development of local-scale population estimation for disease
studies.
Michigan is a particularly good example of a Midwestern state facing emerging and
severe EHD impacts on white-tailed deer populations. Michigan has been host to one of the
largest recent outbreaks at the northern range of EHD while also having one of the longest
histories with the disease at northern latitudes. In 1955, researchers first confirmed epizootic
hemorrhagic disease in Michigan during an outbreak event that was mirrored in New Jersey the
same year (Fay et al. 1956, Shope et al. 1960). Until 2006, only one outbreak occurred after this
3

initial discovery, with approximately 100 dead deer reported as EHD mortalities in 1974
(MDNR). However, since 2006, Michigan had an increase in frequency and severity of EHDrelated outbreaks, with EHD reports confirmed in 2006, 2008, 2009, 2010, 2011, 2012, 2013 and
2016 (MDNR). In the summer of 2012, the most severe EHD outbreak occurred in Michigan,
killing at least 15,000 deer based on voluntary disease reporting and agency confirmation alone
(MDNR). In that year, EHD was confirmed in 30 Michigan counties and mortalities were
reported in 21 other counties during this outbreak. Deer were assumed largely naĂŻve and the
effect of EHD mortality on the deer population was a significant concern for the Michigan
Department of Natural Resources and the public, as deer are a highly valued resource in the state.
The overarching goals of my research are to investigate the local distribution and
population-level impacts on deer in areas affected by EHD, develop and apply population
estimation techniques appropriate for localized population change assessment, and evaluate
regional-level drivers of hemorrhagic disease across an epidemiological gradient. These first two
goals were developed to help reduce uncertainty for deer managers regarding population impact
and recovery of deer from EHD at scales relevant to the MDNR’s ability to manage, and where
the public perceive changes in deer density due to disease mortality. In Chapter 1 I used rigorous
distance sampling estimates at localized scales to compare population abundance and distribution
of deer in areas of Michigan that either experienced severe EHD mortality in 2012 or had no
historical EHD evidence. My objectives were to 1) test whether I could detect a change in
abundance of a deer population attributable to an EHD outbreak in 2012, 2) quantify the finescale spatial extent of EHD impacts associated with riparian habitats and corresponding disease
risk for deer populations, and 3) determine if populations increased over time in response to local
population impacts from EHD. In Chapter 2 I provided another tool for localized population
4

assessment using aerial surveys and the novel application of cutting-edge estimation techniques.
Specifically, the goal of this study was to evaluate an N-mixture model applied to an aerial
survey of unmarked, free-ranging white-tailed deer in an area known to have been affected by
EHD-related mortality. Using aerial surveys and a simultaneous data collection approach, the
objectives for this chapter were 1) to compare detection and abundance parameters estimated
from N-mixture models and distance sampling methods and 2) to test the sensitivity of the spatial
unit size inherent in the study design to the N-mixture approach. Developing and applying
practical estimation tools for local populations will help managers better assess impacts
attributable to disease.
At landscape scales, where environmental or landscape factors may drive the occurrence
of EHD in deer in non-endemic regions more broadly, our goal was to better understand recent
landscape-level patterns of EHD and anticipate future disease outbreaks in emergent zones. In
Chapter 3 I used a longitudinal survey of hemorrhagic disease across the eastern United States
for a 15-year period in which disease occurrence has increased (Nettles et al. 1992, Stallknecht et
al. 2015). Specifically, I evaluated spatial models for hemorrhagic disease that included drought
and related landscape factors predicted to affect disease occurrence through mechanisms
affecting the disease vector. My objectives were to 1) develop a spatiotemporal model to
evaluate if drought severity and related environmental variables explain changing patterns of HD
presence; and 2) to determine if this potential risk factor varies in importance over the present
range of HD in the eastern United States. Finally, the epilogue provides a summary of my
findings and general conclusions for research and management of emergent EHD.

5

LITERATURE CITED

6

LITERATURE CITED

Berry, B. S., K. Magori, A. C. Perofsky, D. E. Stallknecht, and A. W. Park. 2013. Wetland cover
dynamics drive hemorrhagic disease patterns in white-tailed deer in the United States.
Journal of Wildlife Diseases 49: 501-509.
Fay, L. D., A. P. Boyce, and W. G. Youatt. 1956. An epizootic in deer in Michigan,
Transactions of the 21st North American Wildlife Conference: 173.
Fischer, J. R., L. P. Hansen, J. R. Turk, M. A. Miller, W. H. Fales, and H. S. Gosser. 1995. An
epizootic of hemorrhagic disease in white-tailed deer (Odocoileus virginianus) in
Missouri: necropsy findings and population impact. Journal of Wildlife Diseases 31: 3036.
Gaydos, J. K., W. R. Davidson, E.W. Howerth, M. Murphy, F. Elvinger, D. E. Stallknecht.
2002a. Cross-protection between epizootic hemorrhagic disease virus serotypes 1 and 2
in white-tailed deer. Journal of Wildlife Diseases 38: 720-728.
Gaydos, J. K., W. R. Davidson, F. Elvinger, D. G. Mead, E. W. Howerth, and D. E. Stallknecht,
2002b. Innate resistance to epizootic hemorrhagic disease in white-tailed deer. Journal of
Wildlife Diseases 38 4:713-719.
Gaydos, J. K., J. M. Crum, W. R. Davidson, S. S. Cross, S. F. Owen, and D. E. Stallknecht.
2004. Epizootiology of an epizootic hemorrhagic disease outbreak in West Virginia.
Journal of Wildlife Diseases 40: 383-393.
Harvell, C. D., C. E. Mitchell, J. R. Ward, S. Altizer, A. P. Dobson, R. S. Ostfeld, and M. D.
Samuel. 2002. Climate warming and disease risks for terrestrial and marine
biota. Science 296: 2158-2162.
Howerth, E. W., D. E. Stallknecht, and P. D. Kirkland. 2001. Bluetongue, epizootic hemorrhagic
disease, and other orbivirus-related diseases. In: Infectious Diseases of Wild Mammals.
3rd edition. Edited by E.S. Williams and I.K. Barker. Ames, Iowa: Iowa State Press;
2001:77–97.
MDNR Report June 2014 Report #3585 http://www.michigan.gov/dnr/0,4570,7-15310363_48664---00.html (Accessed September 17 2014).
Mellor, P.S., J. Boorman, and M. Baylis. 2000. Culicoides biting midges: their role as arbovirus
vectors. Annual Review of Entomology 45: 307-340.
Nettles V. F., W. R. Davidson, and D. E. Stallknecht. 1992. Surveillance for hemorrhagic disease
in white-tailed deer and other wild ruminants. Proceedings of the Southeast Association
of Fish Wildlife Agencies 46:138-146.

7

Ostfeld, R. S., G. E. Glass, and F. Keesing. 2005. Spatial epidemiology: an emerging or reemerging discipline. Trends in Ecology and Evolution 20:328-336.
Park, A. W., K. Magori, B. A. White, and D. E. Stallknecht. 2013. When more transmission
equals less disease: reconciling the disconnect between disease hotspots and parasite
transmission. PloS one 8:4.
Purse, B. V., P. S. Mellor, D. J. Rogers, A. R. Samuel, P. P. Mertens, and M. Baylis. 2005.
Climate change and the recent emergence of bluetongue in Europe. Nature Reviews
Microbiology, 3:171-181.
Ruder, M. G., T. J. Lysyk, D. E. Stallknecht, L. D. Foil, D. J. Johnson, C.C. Chase, D. A.
Dargatz, and E. P. J. Gibbs. 2015. Transmission and epidemiology of bluetongue and
epizootic hemorrhagic disease in North America: Current perspectives, research gaps,
and future directions. Vector-Borne and Zoonotic Diseases 15: 348-363.
Savini, G., A. Afonso, P. Mellor, I. A. O. Aradaib, H. Yadin, M. Sanaa, W. Wilson, F. Monaco,
and M. Domingo. 2011. Epizootic heamorragic disease. Research in Veterinary
Science 91:1-17.
Shope, R. E., L. G. MacNamara, and R. Mangold. 1960. A virus-induced epizootic hemorrhagic
disease of the Virginia white-tailed deer (Odocoileus virginianus). The Journal of
Experimental Medicine 111:155-170.
Stallknecht, D. E., E. W. Howerth, and J. K. Gaydos. 2002. Hemorrhagic disease in white-tailed
deer: Our current understanding of risk. Transactions of the North American Wildlife and
Natural Resources Conference 67: 75-86.
Stallknecht, D. E., A. B. Allison, A.W. Park, J. E. Phillips, V. H. Goekjian, V. F. Nettles, and J.
R. Fischer. 2015. Apparent increase of reported hemorrhagic disease in the midwestern
and northeastern USA. Journal of Wildlife Diseases 51: 348-361.
Thomas, F. C., N. Willis, and G. Ruckerbauer. 1974. Identification of viruses involved in the
1971 outbreak of hemorrhagic disease in southeastern United States white-tailed deer.
Journal of Wildlife Diseases 10:187-189.
Wilkening, A.J., D. L. Kline, and W. W. Wirth. 1985. An annotated checklist of the
Ceratopogonidae (Diptera) of Florida with a new synonymy. Florida Entomologist 511537.
Wilson, A., and P. Mellor. 2008. Bluetongue in Europe: vectors, epidemiology and climate
change. Parasitology Research 103: S69-77.
Xu, B., M. Madden, D. E. Stallknecht, T. W. Hodler, and K. C. Parker. 2012. Spatial-temporal
model of hemorrhagic disease in white-tailed deer in south-east USA, 1983 to 2000.
Veterinary Record 170: 288.

8

CHAPTER 1: EPIZOOTIC HEMORRHAGIC DISEASE AFFECTS ABUNDANCE AND
DISTRIBUTION OF FREE-RANGING DEER POPULATIONS

INTRODUCTION
Hemorrhagic disease (HD) is the most significant source of viral disease-related mortality
in white-tailed deer (Odocoileus virginianus) in the United States (Howerth et al. 2001). An
increase in the frequency and distribution of reported HD in white-tailed deer has occurred
across northern latitudes in the Midwest and eastern United States since 1980 (Stallknecht et al.
2015). In particular, the emergence of the HD serogroup epizootic hemorrhagic disease virus
(EHDV) has resulted in severe clinical disease and localized mortality of white-tailed deer in the
northern range of the disease where deer are immunologically naive (Howerth et al. 2001,
Gaydos et al. 2004). Epizootic hemorrhagic disease (EHD), the disease caused by EHDV, is
considered endemic in the southeastern US, where infection rates are high but mortality and
morbidity are low (Stallknect et al. 1996). However, where outbreaks have been severe in the
north, resulting in hundreds or thousands of deer mortalities, the impacts of EHD on deer
populations are often unknown and difficult to assess. In these incursive and epidemic zones of
disease emergence, the effects of EHD on the population ecology of white-tailed deer are poorly
understood (Ruder et al. 2015). There is a clear need to better understand the population-level
impacts of EHD in these epidemic zones.
Understanding the impact of epidemic and emerging infectious disease on wildlife
populations is crucial for conservation and management as disease outbreaks intensify and
spread spatially (Harvell et al. 2002, Ostfeld et al. 2005). However, knowing the extent and
population-level impact of emergent wildlife diseases is challenging. In captive white-tailed deer
9

herds where animals are monitored extensively and the spatial extent of the herd is defined,
disease mortality can often be more easily assessed than in free-ranging herds. For example, 46%
average herd mortality was documented in captive herds after an EHD outbreak in 2012, which
is much higher than any free-ranging deer mortality reports (Stevens et al. 2015). While limited
studies exist of EHD-related mortality in free-ranging populations, evidence has shown EHD
outbreaks may reduce deer population abundance at least 20%, with more severe reductions
where the disease occurs in immunologically naĂŻve deer populations (Gaydos et al. 2004). The
immunological response to EHD is related to the history and frequency of past disease outbreaks.
This is in part explained by the proportion of EHDV naive animals within a population and the
subsequent ability for animals to develop neutralizing antibodies and acquire nearly life-long
resistance after exposure (Stallknect et al. 1991, Stallknecht et al. 1996, Coleman et al. 2001,
Howerth et al. 2001, Gaydos et al. 2002). Typically, the turnover rate of a deer population is
faster than the average periodicity of EHD outbreaks in epidemic zones, leaving a small
percentage of the population with immunity to the pathogen and a large susceptible population at
the time of the next outbreak. If these patterns persist, deer populations may be prone to
continued die-offs given optimal environmental conditions for disease transmission. Over time,
deer populations will recover as reproduction continues and new susceptible deer are born or
immigrate into an affected area. However, recovery time will be confounded by the maintenance
of seroprevalence within the population if disease outbreaks continue (Park et al. 2013). This
pattern has been shown at broad scales with HD in the US (Park et al. 2013). However, localized
spatio-temporal patterns of EHD impacts on deer populations have not been well studied in
epidemic regions of EHD (Gaydos et al. 2004).

10

EHD is vector-borne and is reliant on the presence and persistence of a vector species
capable of transmitting the virus to the deer. The vector, an adult biting midge in the genus
Culicoides, transmits EHDV in the late summer across the range of the disease and is tied closely
to climate conditions and available water for life cycle completion (Shope et al. 1960, Mellor et
al. 2000). The effect of EHD on deer populations tends to be localized and often is associated
with features such as wetlands or other sources of water where Culicoides are found (Davidson
2006). Consequently, the population ecology and distribution of the vector plays a significant
role in the epidemiology of hemorrhagic diseases (Mellor et al. 2000, Wilson and Mellor 2008,
Savini et al. 2011). HD reports are often associated with warm temperatures and features such as
wetlands or other sources of water suitable for the life cycle of the vector species (Davidson
2006, Berry et al. 2012, Xu et al. 2012). While the specific species responsible for EHDV
infection in the northern and eastern United States is unclear, we do know that Culicoides require
similar wetland habitats where disease risk is thought to be greater for deer through increased
virus transmission (Ruder et al. 2015). To our knowledge, no study has identified the spatial
extent from a wetland at which deer population impacts occur.
Given the temporal and spatial variability of EHD impacts on deer, assessing population
parameters in the context of an EHD outbreak is one way in which disease impacts on deer
populations may be disentangled. The few studies that have quantified EHD mortality rates or
population impacts have relied on opportunistically radio-collared animals or modeling
population indices at deer management unit scales (Fischer et al. 1995, Beringer et al. 2000,
Gaydos et al. 2004). In the absence of marked animals, Fischer and colleagues (1995) and
Gaydos and colleagues (2004) used population metrics from harvest data at the deer management
unit scale, in conjunction with disease reports, to estimate white-tailed deer mortality rates in
11

Missouri and West Virginia, respectively. However, harvest-based population models at the scale
of deer management units (e.g., ~2,000 km2) have limitations for assessing disease impacts. This
limitation is especially true for the highly localized mortality patterns of EHD, as traditional deer
management units would tend to average local population fluctuations across the unit. Linking
the scale of the environmental process with the scale of appropriate management action is
important for effective regulation and maintenance of the resource (Cumming et al. 2006).
Therefore, estimating the impact of EHD on white-tailed deer populations would be best done at
a local scale relevant to the ecology and population impacts of the disease.
Increased frequency and intensity of EHD-related deer mortalities on the Michigan
landscape since 2006 has generated concern from the public and the Michigan Department of
Natural Resources (MDNR) regarding the effects of EHD on local deer populations. In 2012,
EHD was confirmed in 30 Michigan counties and mortalities were reported in 21 other counties
resulting in approximately 15,000 EHD suspected deer mortalities (MDNR). The first cases of
EHD in Michigan occurred late summer of 1955, appearing in four counties across the central
lower peninsula of the state (Shope et al. 1960). Another outbreak occurred in 1974 in 5 counties
in Michigan, and did not reappear until 2006, when it was found in a single county in the
southwest portion of the state (MDNR). Since 2006, EHD was confirmed in 2008-2013, with
2012 being the largest known outbreak (MDNR). Large numbers of dead deer were confirmed to
have died from EHDV through laboratory virus isolation in 2012, although not all deer were
tested and EHD cases were likely underreported (MDNR). The spatial distribution and relative
intensity at which a deer population is affected by EHD is particularly concerning as the sudden
reduction in deer abundance may have long lasting effects at local scales.

12

Due to the emergent nature of this disease in similar northern latitudes over the last 10
years and the expected relationship with midge habitats, there is a clear need for information on
localized population impacts of the disease on deer. Our objectives were to 1) test whether I
would detect a change in abundance of a deer population attributable to an EHD outbreak, 2)
quantify the fine-scale spatial extent of EHD impacts on deer populations associated with
riparian habitats and corresponding disease risk, and 3) determine if populations increased over
time in response to local population impacts from EHD. I hypothesized that EHD has localized
population impacts on white-tailed deer and that those impacts are habitat dependent. Previous
research investigating habitat use of deer in agriculture-forest systems has shown that deer select
for lowland deciduous-forest cover, often associated with wetlands (Hiller et al. 2009). Thus,
deer abundance in a healthy population (as defined by an absence of EHD outbreaks) should
have greater local abundance in forested and wetland areas compared to adjacent agricultural
areas with limited forest cover. However, in an EHD impacted region, deer associated with
wetlands would be more susceptible to disease-related mortality. If that mortality was severe,
and abundance is assessed with sufficient precision and at appropriate scales, local population
abundance of deer should be lower in forested wetland areas compared to adjacent areas
dominated by agriculture.
STUDY AREA
Two study areas in southern Michigan were selected for a natural case-control
comparison based on past disease occurrence, landscape metrics relevant to the disease process,
and deer abundance. Deer hunting occurred every fall in both study areas and deer harvest
regulations were similar between study areas. The Maple River study area (190 km2) within
Montcalm, Ionia, Clinton, and Gratiot counties in Michigan, USA was the site of a major deer
13

mortality event due to EHD in 2012 and served as our case study site (Fig. 1.1 and Appendix A).
This study area had notable EHD mortality reports in 2012 with Ionia county reporting 2,498,
Montcalm county reporting 910, Clinton county reporting 514, and Gratiot county reporting 319
EHD deer mortalities during the outbreak period (July – October, MDNR unpublished data).
Buck harvest, an index of county-level deer abundance prior to data collection, was 2,871, 2,614,
3,400, and 6,442 for Ionia, Clinton, Gratiot, and Montcalm counties, respectively for 2012. This
area is in the Lansing loamy plain ecoregion characterized by a flat landscape supporting a mix
of forest cover and agriculture consisting primarily of soybeans, grain, corn, and fragmented
woodlots. Greater forest cover and water exist near the riparian corridor of the Maple River. The
climate has large seasonal variation in temperatures ranging from an average low of -10 degrees
C to an average high of 28 degrees C, with 90.42 cm of average annual precipitation.
The Cass River study area (203 km2) lies within Tuscola, Sanilac, and Huron counties in
Michigan USA and is similar to the Maple River area in land cover composition and climate
(Fig. 1.2). The area is in the Maumee Lake plain ecoregion and dominated by corn, grain, and
soy-based agriculture beyond the riparian corridor. Although a total of 8 EHD deer mortalities
were reported for Tuscola and Sanilac counties in 2012, none were reported in townships within
our study area and none were reported for Huron county, making this site likely EHD-free. The
buck harvest for Tuscola and Sanilac and Huron counties in 2012 was 4,619, 6,312, and 3,942
respectively. The climate is slightly less variable than in the Maple River study site, with
temperatures ranging from an average low of 1.7 degrees C to 13.3 degrees C, with 85.69 cm of
average annual precipitation.

14

Figure 1.1 The location of line transects used to estimate fine-scale deer abundance at varying
distances from a riparian corridor in two study areas in Michigan, from 2013-2016. The Maple
River area (red transects) had previously reported severe EHD outbreaks (case site) and the Cass
River area (black transects) had no record of past EHD outbreaks (control site).

METHODS
Data collection and study design
To evaluate the role of EHD mortality specifically on local deer abundance and test our ability to
estimate abundance at scales relevant to EHD, I surveyed populations at both the EHD-affected
site and the EHD-free site. To detect heterogeneity and change in abundance of white-tailed deer
at scales relevant to EHD within each site, I designed our study to assess deer abundance within
localized habitat arrangements adjacent to riparian areas and in years following a documented
disease outbreak (MDNR 2014). I used line-transect distance sampling to estimate deer

15

abundance for each stratum within each site (Thomas et al. 2010). I placed 2 line-transects within
approximately 1 km of the river, and 2 line-transects approximately 5km from the river (Fig 1.1).
Survey routes were approximately 24 km long line-transects positioned along low-intensity dirt
roads in an adjusted zig-zag pattern (Marques et al. 2001), following a major river corridor
within the study area. The observers made counts of deer from motor vehicles to allow for
observer detection across maximum distances and habitats. Observers maintained a speed that
would ensure that deer were not double counted if running but slow enough to ensure optimal
detection of deer along the transect (~20 km/h; Zamboni et al. 2015). The high proportion of
private landownership and lack of accessibility for non-motorized vehicles made random transect
placement impractical in this study. I selected roads for line-transects based on characteristics
that minimized the potential of a negative behavioral response by deer across open agricultural
landscapes (Meisingset et al. 2014). All transects were tested for driving conditions and visibility
prior to survey initiation. Transects were then surveyed during a time period in which the
population could be assumed closed, after fawning (mid-June) season and before hunting (midSeptember) season between 2013-2016 for the Maple River study area and 2014-2016 for the
Cass River study area. Further, movement of white-tailed deer resulting in population change
through immigration and emigration was assumed to be minimal during this period (Long et al.
2008). Low deer density may result in fewer deer available for detection; thus I expected low
detection probabilities in the stratum that had recently experienced the greatest disease mortality.
Subsequently, survey replicates were increased in these areas to ensure an appropriate sample
size and encounter rate. The line-transects were grouped as habitat-based strata in each study
area, with each stratum encompassing 1 of the 2 dominant habitat types associated with
proximity to a riparian corridor, hereon designated as “wetland” and “agriculture” strata. The

16

total area for estimation was calculated by creating a 1000 m buffer around each line-transect,
and pooling similar habitat types. The agricultural stratum in the Maple River study area was
surveyed 6 times in 2013, 8 times in 2014, 8 times in 2015, and 7 times in 2016. The wetland
stratum in the Maple River study area was surveyed 8 times in 2013, 10 times in 2014, 10 times
in 2015, and 9 times in 2016. The Cass River study area also contained 2 habitat-based stratum,
agricultural and wetland, which were surveyed 6 times each in 2014 and 2015, and 8 times each
in 2016. Evening surveys were completed within a 3-hour time period and were conducted once
per week throughout the data collection period for replication. Finally, line-transects were placed
approximately 3 km from each other or were separated by a geographical barrier (large river) to
maintain independence between transects at each site.
Free-ranging unmarked deer were counted through ground-based observations on both
sides of the transect line beginning within 3 hours prior to sunset and concluding when the
survey route was complete or when it was too dark to detect deer. Assuming a minimum number
of 2 observers per line-transect were present, all transects in a study area were surveyed
simultaneously on a given night. Deer observations were measured with range finders to obtain
measures of distance from observers to solitary deer or to the center of a group of deer. A group
was defined as deer within approximately 6m of each other. The unit of observation was the deer
group to ensure independence in the analyses (Buckland et al. 2001, Lovely et al. 2013). The use
of range finders was important to meet the assumption of accuracy in distance measurements. A
team of 2-4 observers, including the vehicle driver, searched ahead and perpendicular to the
vehicle, stopping the vehicle to use binoculars to verify detection of deer. Direction of travel was
alternated among replicate survey days for a given transect to account for potential differences in
detection by time of day. Observations of deer or groups of deer were made or confirmed with
17

10x30 binoculars or unimpeded eye only (no spotlights or thermal imagery were used). If the
observer was not perpendicular to the deer observation, a radial angle measurement was taken
from the vehicle with a planar protracting device fitted perpendicularly along the vehicle
window. Data related to the detection process were collected, including group size, vegetation
type (row crop, pasture, forest, residential, other), behavior upon detection (bedded, standing,
feeding, moving), observer (biologist or non-biologist team), and sky condition (clear skies,
cloudy, light precipitation, fog). Time of day, distance along the transect, location coordinates,
and total transect effort were also collected for analysis purposes. Driving conditions were
recorded, with conditions of excessive precipitation, fog, or thunderstorms precluded from
survey nights. Surveys were coordinated with the Michigan Department of Natural Resources
(MDNR), volunteers from local deer-hunting organizations, and students. Observers were
categorized as either biologists or not. Regardless of experience, all observers were trained prior
to conducting surveys, and protocol was reiterated prior to each survey to ensure consistent data
collection. Because an underlying assumption of distance sampling requires that all objects that
occur on the line-transect (at distance zero) are detected, I repeatedly instructed all observers to
ensure that no objects on or directly adjacent to the line transect were missed. Surveys were not
conducted more than once in a 24-hour period as to avoid double counting and maintain
independence.
Distance sampling analysis
I analyzed data with conventional distance sampling and multiple covariate distance
sampling engines in Program Distance v. 6.2 as described by Thomas et al. (2010). I analyzed
data for each stratum in each study area and each year independently using the Program Distance
data filter settings, as detection was not assumed homogeneous across years or strata (Buckland
18

et al. 2009). After conducting exploratory analyses, I truncated all data beyond 450 meters for all
transects in all years for improved model fit (Buckland et al. 2001, Williams and Thomas 2007)
and for consistency in population inference at a fine-scale across years. I modeled the detection
probability function following recommendations of Buckland et al. (2001) and considered
variables that were predicted to influence the detection process for each stratum and year but that
were not highly correlated with distance from the transect line. I evaluated half-normal and
hazard-rate key functions in a step-wise approach and used appropriate expansion terms (cosine
or polynomial) provided by Buckland et al. (2001) in combination with up to 3 adjustment terms.
Key functions and adjustment terms were selected for each stratum and evaluated with Akaike
Information Criterion (AIC) to select the most parsimonious model which best explained
detection probability (Burnham and Anderson 2002). I used program DISTANCE to select the
appropriate adjustment terms based on minimum AIC for either half-normal or hazard rate key
functions with a cosine series expansion (Buckland et al. 2001). Further, I then determined the
best cut points for distance intervals from the transect, given the covariates, using goodness-of-fit
testing of observed and expected values. Once the most robust key function for the data was
identified, I systematically began including covariates in the MCDS engine and bootstrapped the
variance when appropriate (e.g., cluster size; Buckland et al. 2001). The covariates used in our
model set were coded as factors in the MCDS analysis and included observer team (Obs),
weather conditions (Skycode), habitat (Veg), behavior of deer upon detection (Behave), and the
group size of the observed deer group (Cluster). For the purposes of analysis, the initial deer
detection for an individual or group made by an observer was assigned the appropriate covariate
information. Each variable was modeled individually and as a combination of potential
interactions developed with an a priori model set and implemented with the model definition

19

properties (Table 1.1). This suit of multivariate models was then evaluated with AIC to select the
most parsimonious detection probability model for each year and stratum, assuming model
convergence. Of the best models, I assessed the diagnostic histograms, qq-plots, and the %
coefficient of variation to ensure our data met the assumptions of distance sampling and
parameter estimates had adequate precision needed for population inferences. I evaluated
Kolmogorov-Smirnov goodness-of-fit to further verify models selected as the most likely to fit
the data (Buckland et al. 2004). Once the model with greatest likelihood was selected for each
stratum in each year, I estimated abundance, detection probability, and variance associated with
point-estimates. I derived density by dividing the abundance estimate by the respective stratum
area. Because our desired level of inference was at the stratum level within each study site,
indicating the recent spatial distribution of disease mortality impacts, I did not attempt to
combine estimates for the entire study area.
Table 1.1 Notation and description for a priori distance sampling models for estimating
abundance and density estimates of white-tailed deer in the Maple River and Cass River study
areas of Michigan, 2013-2016.
Detection Variables
Cluster

Description
The number of deer counted within a group

Sky

Weather condition on the given survey night (clear, partly cloudy, overcast,
fog, light rain)

Veg

Habitat cover (Row crop, pasture, forest, residential, other)

Obs

Observer team; inclusion of a biologist or not

Behave
Cluster + Behave

Behavior of deer upon initial detection (running, walking, standing, feeding,
bedded)
Interaction between deer group size and deer behavior

Cluster + Obs

Interaction between deer group size and detection by observer

Cluster + Sky

Interaction between deer group size and weather conditions during survey

Null

Hazard rate or Half-normal key function only (depending on minimum AIC)

20

RESULTS
Data collection
Our total survey effort was 3,323 km in our EHD site (Maple River study area) over the
course of 4 summers (2013-2016), and 1,993 km in our control site (Cass River study area) over
the course of 3 summers (2014-2016). Surveys resulted in a total of 1,481 observations and 2,725
observations in the EHD site and the control site, respectively, with observations being a single
deer or a group of deer within approximately 6m of each other. This group classification was
similar to Lovely et al. (2013) distance sampling protocol of white-tailed deer in Virginia and
based on prior group size research of white-tailed deer (Lagory 1986). The range of mean group
(or Cluster) size for the EHD site across years was 1.69 – 2.69 (SE = 0.09 – 0.21) deer and was
2.29 – 2.98 (SE = 0.11 – 0.16) deer for the control site. However, because the expected group
size occasionally suggested size bias, I used the size-bias cluster estimation method for the deer
density and abundance parameter estimation.
Model selection
The half-normal key function was selected using minimum AIC for all except 2 years in a
single stratum when the hazard-rate key function was selected (agriculture stratum in EHD site
2015 and wetland stratum in control site 2015) (Table 1.2). After conducting exploratory
analysis and testing models for convergence, the habitat type covariate was removed from further
model evaluation as a variable that may affect detection function estimation due to high
correlation with distance from the transect line. Habitat type has often been cited as an important
predictor of detection functions for deer (Lovely et al. 2013, Zamboni et al. 2015), but model
convergence rarely occurred when the habitat covariate was included and correlation issues were
21

apparent relative to distance and habitat composition associated with transect lines (low-intensity
roads). The models that garnered the highest support from the data with minimum AIC criteria,
provided reliable model convergence, and best explained the variance in detection probability
included covariates for group size or deer behavior, but varied by year, study area, and stratum
(Table 1.2). The weather condition variable (skycode) was important in 2016 for describing the
variance in detection for the agriculture stratum in both study areas, and the wetland stratum for
Maple River in 2015 (Table 1.2). No other covariates were selected using AIC, indicating
observer effects were not present in our data (Table 1.2).

22

Table 1.2 Models selected for each stratum in each year and corresponding selection criteria and data descriptors.

Year

Study Area

Stratum

Model; Key function and No.
adjustment terms

k

GOF KS pvalue

ESW (est.
strip
width)

5

0.0804

323.49

286

181.8

5

0.2396

311.08

416

186.9

Tot. No.
observations

Tot. Effort

2014

Cass River

Ag

2014

Cass River

Wetland

Cluster + Behave; Half-normal,
0
Cluster + Behave; Half-normal,
0

2015
2015

Cass River
Cass River

Ag
Wetland

Cluster; Half-normal, 0
Cluster; Hazard-rate, 2

3
5

0.1734
0.2416

259.78
373.9

396
572

191.6
180

2016
2016

Cass River
Cass River

Ag
Wetland

Clust + Skycode; Half-normal, 3
Clust + Behave; Half-normal, 3

6
9

0.06
0.1911

306.52
298.95

479
576

241.9
256

2013
2013

Maple River
Maple River

Ag
Wetland

4
3

0.792
0.0152

273.45
192.96

83
59

171.1
227.3

2014
2014

Maple River
Maple River

Ag
Wetland

5
2

0.5838
0.7168

306.65
290.37

251
170

245.8
309.4

2015
2015

Maple River
Maple River

Ag
Wetland

Behavior; Half-normal, 0
Behavior; Half-normal, 0
Cluster + Behave; Half-normal,
0
Cluster; Half-normal, 0
Cluster; Hazard-rate Simple
poly, 0
Clust + Skycode; Half-normal, 0

3
4

0.118
0.0404

260.72
304.51

229
185

246.4
313

2016
2016

Maple River
Maple River

Ag
Wetland

Clust + Skycode; Half-normal, 0
Clust + Behave; Half-normal, 0

5
7

0.1638
0.142

298.16
235.63

259
245

240.8
311.2

23

Deer density and abundance estimates
The best model for each year and stratum based on the minimum AIC and model fit
criteria was used to estimate a detection function that was specific for annual deer population
abundance and density estimates and associated 95% confidence intervals in each stratum (Table
1.2, Fig. 1.2). All estimates had sufficient minimum observations as recommended by Buckland
et al. (2001) and met goodness-of-fit criteria and detection function fit criteria (Table 1.2-1.3,
Appendix B). In the EHD affected area, I found deer abundance was significantly higher in the
agricultural stratum in 2013 (281; CI: 165-401) and 2014 (220; CI: 151.18-301.79) than in the
wetland stratum for 2013 (60; CI: 42.52-78.7) and 2014 (107; CI: 92.59-133.25). Further, deer
population density estimates were lowest in the wetland stratum of the EHD site for all years
compared to either stratum in the unaffected study area, including the minimum density estimate
(Fig. 1.3; 0.63; CI: 0.44-0.82 in 2013 affected area vs. 8.94; CI: 6.68-11.66 in 2014 unaffected
area). I found the wetland stratum in the EHD site was the only stratum to experience consistent
population increases over the course of the 4 year study, with abundance estimates doubling in
2016 compared to 2013 (Table 1.3, Fig. 1.3). Population abundance in the control site (Cass
River study area) did not change over time (Fig. 1.4). Further, the wetland stratum was not
different from the agricultural stratum except in 2014, when deer density was greater in wetland
habitat (580.39; 95% confidence interval (CI): 539.31-628.80) than in the agriculture habitat
(344.68; 95% confidence interval (CI): 257.33-449.30).

24

Table 1.3 Parameter estimates from best fitting models in each stratum, study area, and year for 2013-2016 in Michigan. Asterisks
represent 95% confidence intervals from bootstrapping.
Estimate

Year
2014
2014
2015
2015
2016
2016
2013
2013
2014
2014
2015
2015
2016
2016

Study
Area
Cass
River
Cass
River
Cass
River
Cass
River
Cass
River
Cass
River
Maple
River
Maple
River
Maple
River
Maple
River
Maple
River
Maple
River
Maple
River
Maple
River

Total
area
(km.2) of
estimation

Abundance

Ag

99.82

Wetland

N 95%CI*

N 95%CI*

N 95% CI*

Detection
Prob.

LCI

UCI

4.50

0.72

0.67

0.77

5.21

6.07

0.69

0.65

0.73

5.14

3.59

6.72

0.58

0.54

0.62

854.37

5.58

4.99

8.25

0.83

0.79

0.88

364.98

532.38

4.59

3.66

5.34

0.68

0.65

0.72

569.54

460.29

661.38

5.50

4.44

6.39

0.66

0.63

0.70

94.07

280.68

165.00

401.00

1.49

0.88

2.13

0.61

0.52

0.71

Wetland

95.86

60.24

42.52

78.7

0.63

0.44

0.82

0.46

0.38

0.56

Ag

94.07

219.98

151.18

301.79

2.34

1.61

3.21

0.68

0.63

0.74

Wetland

95.86

106.62

92.59

133.25

1.11

0.97

1.39

0.65

0.58

0.71

Ag

94.07

184.96

109.56

250.92

1.97

1.17

2.67

0.58

0.53

0.63

Wetland

95.86

125.53

105.25

169.16

1.31

1.10

1.77

0.68

0.62

0.74

Ag

94.07

210.53

124.86

306.61

2.24

1.33

3.26

0.66

0.60

0.71

Wetland

95.86

206.07

178.54

252.53

2.15

1.86

2.64

0.52

0.47

0.59

Stratum

LCI

UCI

Density

344.68

257.33

449.30

3.45

2.58

103.6

580.39

539.31

628.80

5.60

Ag

99.82

532.00

371.26

695.99

Wetland

103.6

577.76

516.32

Ag

99.82

458.10

Wetland

103.6

Ag

25

LCI

UCI

3.50
2.24

2.34
3.00
1.97

2.15

2.50

Deer/ km2

1.49
2.00
1.31
1.50

1.11
0.63

1.00

0.50

0.00
2013

2014

2015

2016

2013

Agricultural stratum

2014

2015

2016

Wetland stratum

Figure 1.2 Deer density and 95% confidence interval estimates (deer/sq. km) and for
the Maple River study area (EHD site) for 2013-2016 in Michigan.

26

9.00
5.58

8.00
7.00

5.14
5.60

Deer / sq. km2

6.00

5.50

4.59

5.00
3.45
4.00
3.00
2.00
1.00
0.00
2014

2015

2016

2014

Agricultural stratum

2015

Wetland stratum

Figure 1.3. Deer density and 95% confidence interval estimates (deer/sq. km) and
for the Cass River study area (control site) for 2014-2016 in Michigan.

27

2016

DISCUSSION
Understanding the spatial and temporal dynamics of emerging epizootic hemorrhagic
disease in northern latitudes is becoming increasingly important for wildlife conservation and
management. Disentangling these impacts across space and time as the disease expands
northward is complicated by variation in host infection rates, spatiotemporal host immunity
patterns, vector dynamics, reporting and sampling bias, and host population dynamics.
(Stallknecht et al. 2015, Ruder et al. 2015). Our findings show that I could detect differences in
the population distribution of deer within 1 km from a riparian corridor of known EHD
occurrence in 2012 compared to a population only 5 km from the same riparian corridor, with
reduced abundance nearest the river. In the control site there was no difference in deer
population over time or at varying distance from the riparian corridor. Importantly, the local deer
population most affected by EHD was the only population in our study that had a consistent and
measurable population increases across the 4 subsequent years following an EHD outbreak.
These findings lead us to conclude that EHD was the driving force behind the population
changes I detected in the riparian corridor of Maple River. To our knowledge, this is the first
study to provide empirical evidence that EHD has fine-scale population impacts on white tailed
deer, and those impacts are habitat dependent.
My findings suggest that epizootic hemorrhagic disease had significant impacts on the
deer population in the riparian corridor of the Maple River study area. Populations in this area
were perceived to have been heavily impacted by the 2012 outbreak due to large numbers of
dead deer observed and reported during the 2012 outbreak. Within the same affected study area,
I found an increase in deer abundance from the summer of 2013 to 2016 in the wetland stratum
nearest the riparian corridor were EHD had been most severe in a 2012 outbreak, and slightly
28

higher deer density, but with no population change, for populations measured farther from the
riparian corridor (within 5km) (Fig. 1.3). Further, there was no change in deer density, with
greater variation and higher overall density estimates, in the two stratum with no past EHD
report history (Fig. 1.4). In the stratum nearest the riparian corridor in the control site, mean deer
densities were higher, which corresponded with predictions for increased deer habitat selection
and absence of EHD impacts for this area. Our research suggests that localized EHD mortality
was the driving force behind the observed population variation and subsequent population
change in two otherwise comparable landscapes.
EHD outbreaks have long been assumed to be spatially localized, but only a few previous
studies have attempted to explicitly quantify the distribution and spatial scale of an outbreak for
free-ranging deer at a population level (Fischer et al. 1995, Gaydos et al. 2004). In our study, the
Maple River study area had documented cases of numerous EHD mortalities in 2012 and was
likely a completely naive population to EHDV, as the last reported EHD case in the area was in
1974 (MDNR 2012). I detected differences in deer populations associated with riparian corridors
and wetlands, which have been found to affect EHD dynamics and in endemic regions (Berry et
al. 2013) and are critical for the larval stage of Culicoides spp. (Mellor et al. 2000). This
localized impact of EHD is likely due to a combination of mechanisms affecting observed
mortality in wild deer, including herd immunity dynamics, association of the disease with its
vector (Culicoides), and reporting of detected animals. The abundance and distribution of
Culicoides species capable of transmitting EHDV to deer is a key factor in the presence of EHD
across the range of the disease (Mullens et al. 2004, Savini et al. 2011, Ruder et al. 2015). It is
likely that the distribution and habitat requirements for the Culicoides vector is the limiting
factor driving localized EHD impacts on deer in epizootic disease regions.
29

The differences in deer density between the control site and the EHD site in our study,
particularly when comparing the two wetland stratum representing high-quality habitat, provided
another indication that EHD mortality in 2012 was driving the observed population variation.
Previous research investigating deer habitat use in agricultural-forest landscapes has shown that
deer select for forested riparian corridors (Walter et al. 2011), that greater deer densities exist in
riparian corridors compared to areas with less riparian habitat (Compton et al. 1988), and that
during dispersal movements, deer selectively move within riparian corridors (Clements et al.
2011). Thus, deer abundance in a healthy deer population (as defined by absence of EHD
outbreaks in this case) should have greater local abundance in forested areas associated with
wetlands compared to mostly agricultural areas lacking forest cover. I found the lowest estimated
deer densities across both study areas and all stratum was the wetland stratum nearest the
riparian corridor in otherwise highly suitable habitat (Fig. 1.3). This habitat-density pattern was
not the case in the wetland habitat of the control site, where deer density was greater than or
equal to deer densities in agricultural habitat. Wetland or riparian corridors in agriculture-forest
landscapes thus pose increased EHD risk for deer, given the significant reduction of deer density
in the EHD affected wetland area of Maple River.
Deer population density increased in a consistent manner across the 4 years of the study
in the wetland stratum of the Maple River study area (Fig. 1.3). This was the only stratum that
had a consistent upward trajectory for deer abundance and density parameter estimates. I argue
that the observed population increase was due to the sudden release of intraspecific competition
in a resource abundant habitat due to a severe mortality event from EHD. Assuming EHD
mortality did significantly reduce the deer population in this localized area from pre-outbreak
densities, deer ecology would predict that maximized recruitment rates likely drove the
30

population increase observed in our study. For many species, including white-tailed deer,
population growth rates are optimized when populations are low due to density dependent
mechanisms, which are a function of recruitment relative to population carrying capacity
(Milner-Gulland and Akcakaya 2001, McCullough 1979). In K-selected species, the greatest
opportunity to observe maximum recruitment per female is shortly after a significant population
reduction and subsequent release from intraspecific competition (McCullough 1979). Further,
minimal immigration into the population would be expected within the time interval of this study
and given the low proportion of dispersing individuals (yearling males) in a white-tailed deer
population and the highly philopatric nature of the female portion of the population (Porter et al.
1991).
My findings provide clear spatial distributions associated with EHD impacts at the deer
population level and evidence of population recovery where EHD impacts occurred. Populationlevel effects of EHD observed in my study are likely similar to effects on white-tailed deer
populations in other epidemic regions of EHD occurrence. Recent patterns of EHD have
demonstrated increased frequency and intensity of outbreaks (Stallknecht et al. 2015), resulting
in increasing concern over deer resource management and population declines. I demonstrate
that while disease impacts may result in significant changes in deer populations, severe mortality
will likely remain localized in association with wetland habitats capable of hosting vector
populations. Further, deer populations are expected to increase significantly post-outbreak,
assuming abundant food resources and limited mortality factors exist over time.
Epizootic hemorrhagic disease is a significant source of non-harvest deer mortality across
the geographic range of the disease (Shope et al. 1960, Stallknecht et al. 2015). Although
historically a significant disease for white-tailed deer in the Southeastern and Great Plains
31

regions of the United States, EHD has more recently emerged at increasingly northern latitudes,
affecting deer populations at greater frequency and intensity, generating considerable uncertainty
regarding the population impacts of the disease on deer (Stallknecht et al. 2015). The emergence
of this disease is complex and may have negative consequences for white-tailed deer populations
across much of the northern latitudes of the US where deer are currently naĂŻve to EHDV. This
study provided information on the scale at which EHD impacts should be addressed, and helped
inform expectations for impacts of future outbreaks. If EHD occurrence continues to occur, but
at a frequency not conducive for the development of herd immunity at landscape scales,
significant mortality events will continue and should be addressed when planning deer
management objectives. Further, the collection of consistent seroprevalence data from deer
across the emerging range of the disease would provide managers with better information on the
susceptibility of deer populations to EHDV and expected mortality impacts.

32

APPENDICES

33

APPENDIX A: SUPPLEMENTARY EHD SURVEILLANCE MAP OF DATA REPORTED IN
2012 IN MICHIGAN. MAP PRODUCED BY THE MICHIGAN DEPARTMENT OF
NATURAL RESOURCES.
Figure A.1 A map of Epizootic Hemorrhagic Disease occurrence in Michigan 2012.

34

APPENDIX B: SUPPLEMENTARY DETECTION PROBABILTY HISTOGRAMS FOR
DISTANCE SAMPLING ESTIMATES IN CHAPTER 1.
Figure B.1 Fitted detection probability histograms of deer observations at binned distance
intervals for the Cass River and Maple River study areas from 2014-2016 and 2013-2016,
respectively, in Michigan.

35

LITERATURE CITED

36

LITERATURE CITED

Beringer, J., L. P. Hansen, and D. E. Stallknecht. 2000. An epizootic of hemorrhagic disease in
white-tailed deer in Missouri. Journal of Wildlife Diseases 36: 588-591.
Berry, B. S., K. Magori, A. C. Perofsky, D. E. Stallknecht, and A. W. Park. 2013. Wetland cover
dynamics drive hemorrhagic disease patterns in white-tailed deer in the United States.
Journal of Wildlife Disease 49: 501-509.
Buckland, S. T., D. R. Anderson, K. P. Burnham, J. L. Laake, D. L. Borchers, and L. Thomas.
2001. Introduction to distance sampling: estimating abundance of biological populations.
Oxford University Press, Oxford, United Kingdom.
Buckland, S. T., D. R. Anderson, K. P. Burnham, J. L. Laake, D. L. Borchers, and L. Thomas.
2004. Advanced distance sampling. Oxford University Press, Oxford, United Kingdom.
Buckland, S. T., R. E. Russell, B. G. Dickson, V. A. Saab, D. N. Gorman, and W. M. Block.
2009. Analyzing designed experiments in distance sampling. Journal of Agricultural,
Biological, and Environmental Statistics 14: 432-442.
Burnham, K. P., and D. R. Anderson. 2002. Model selection and inference: a practical
information-theoretic approach. Springer Publishing, Second edition, New York, New
York, USA.
Clements, G. M., S. E. Hygnstrom, J. M. Gilsdorf, D. M. Baasch, M. J. Clements, and K.C.
Vercauteren. 2011. Movements of white‐tailed deer in riparian habitat: Implications for
infectious diseases. The Journal of Wildlife Management 75: 1436-1442.
Coleman, P. G., B. D. Perry, and M. E. J. Woolhouse. 2001. Endemic stability—a veterinary
idea applied to human public health. The Lancet 357: 1284-1286.
Compton, B. B., R. J. Mackie, and G. L. Dusek. 1988. Factors influencing distribution of whitetailed deer in riparian habitats. The Journal of Wildlife Management 52: 544-548.
Cumming, G. S., D. H. M. Cumming, and C. L. Redman. 2006. Scale mismatches in socialecological systems: causes, consequences, and solutions. Ecology and Society 11: 14.
Davidson, W. R. 2006. White-tailed Deer (Odocoileus virginianus). Field Manual of Wildlife
Diseases in the Southeastern United States 3rd ed: Athens, GA, Southeastern Cooperative
Wildlife Disease Study: 26-99.
Fischer, J. R., L. P. Hansen, J. R. Turk, M. A. Miller, W. H. Fales, and H. S. Gosser. 1995. An
epizootic of hemorrhagic disease in white-tailed deer (Odocoileus virginianus) in
Missouri: necropsy findings and population impact. Journal of Wildlife Diseases 31: 3036.
37

Gaydos, J. K., W. R. Davidson, E.W. Howerth, M. Murphy, F. Elvinger, D. E. Stallknecht. 2002.
Cross-protection between epizootic hemorrhagic disease virus serotypes 1 and 2 in whitetailed deer. Journal of Wildlife Diseases 38: 720-728.
Gaydos, J. K., J. M. Crum, W. R. Davidson, S. S. Cross, S. F. Owen, and D. E. Stallknecht.
2004. Epizootiology of an epizootic hemorrhagic disease outbreak in West Virginia.
Journal of Wildlife Diseases 40: 383-393.
Harvell, C. D., C. E. Mitchell, J. R. Ward, S. Altizer, A. P. Dobson, R. S. Ostfeld, and M. D.
Samuel. 2002. Climate warming and disease risks for terrestrial and marine
biota. Science 296: 2158-2162.
Hiller, T.L., H. Campa III, and S. R. Winterstein. 2009. Estimation and implications of space use
for white-tailed deer management in southern Michigan. Journal of Wildlife
Management 73: 201-209.
Howerth, E. W., D. E. Stallknecht, and P. D. Kirkland. 2001. Bluetongue, epizootic hemorrhagic
disease, and other orbivirus-related diseases. In: Infectious Diseases of Wild Mammals.
3rd edition. Edited by E.S. Williams and I.K. Barker. Ames, Iowa: Iowa State Press;
2001: 77–97.
Lagory, K. E. 1986. Habitat, group size, and the behaviour of white-tailed deer. Behaviour 98:
168-179.
Long, E. S., D. R. Diefenbach, C. S. Rosenberry, and B. D. Wallingford. 2008. Multiple
proximate and ultimate causes of natal dispersal in white-tailed deer. Behavioral Ecology
19: 1235-1242.
Lovely, K. R., W. J. Mcshea, N. W. Lafon, and D. E. Carr. 2013. Land parcelization and deer
population densities in a rural county of Virginia. Wildlife Society Bulletin 37: 360-367.
Marques, F. F., S. T. Buckland, D. Goffin, C. E. Dixon, D. L. Borchers, B. A. Mayle, and A. J.
Peace. 2001. Estimating deer abundance from line transect surveys of dung: sika deer in
southern Scotland. Journal of Applied Ecology 38: 349-363.
McCullough, D. R. 1979. The George Reserve deer herd: population ecology of a K-selected
species. University of Michigan Press.
MDNR Report June 2014 Report #3585 http://www.michigan.gov/dnr/0,4570,7-15310363_48664---00.html (Accessed September 17 2014).
Meisingset, E. L., L. E. Loe, Ø. Brekkum, B. Van Moorter, and A. Mysterud. 2013. Red deer
habitat selection and movements in relation to roads. The Journal of Wildlife
Management 77: 181-191.

38

Mellor, P. S., J. Boorman, and M. Baylis. 2000. Culicoides biting midges: their role as arbovirus
vectors. Annual review of entomology 45: 307-340.
Milner-Gulland, E. J., and H. R. Akçakaya. 2001. Sustainability indices for exploited
populations. Trends in Ecology & Evolution 16: 686-692.
Mullens, B. A., A. C. Gerry, T. J. Lysyk, and E. T. Schmidtmann. 2004. Environmental effects
on vector competence and virogenesis of bluetongue virus in Culicoides: interpreting
laboratory data in the field. Veterinaria Italiana 40: 160-163.
Nettles V. F., W. R. Davidson, and D. E. Stallknecht. 1992. Surveillance for hemorrhagic disease
in white-tailed deer and other wild ruminants. Proceedings of the Southeast Association
of Fish Wildlife Agencies 46: 138-146.
Ostfeld, R. S., G. E. Glass, and F. Keesing. 2005. Spatial epidemiology: an emerging or reemerging discipline. Trends in Ecology and Evolution 20: 328-336.
Park, A. W., K. Magori, B. A. White, and D. E. Stallknecht. 2013. When more transmission
equals less disease: reconciling the disconnect between disease hotspots and parasite
transmission. PloS one 8: 4.
Porter, W. F., N. E. Mathews, H. B. Underwood, R. W. Sage, and D. F. Behrend. 1991. Social
organization in deer: implications for localized management. Environmental
Management 15: 809-814.
Ruder, M. G., T. J. Lysyk, D. E. Stallknecht, L. D. Foil, D. J. Johnson, C.C. Chase, D. A.
Dargatz, and E. P. J. Gibbs. 2015. Transmission and epidemiology of bluetongue and
epizootic hemorrhagic disease in North America: Current perspectives, research gaps,
and future directions. Vector-Borne and Zoonotic Diseases 15: 348-363.
Savini, G., A. Afonso, P. Mellor, I. A. O. Aradaib, H. Yadin, M. Sanaa, W. Wilson, F. Monaco,
and M. Domingo. 2011. Epizootic heamorragic disease. Research in Veterinary
Science 91: 1-17.
Shope, R. E., L. G. MacNamara, and R. Mangold. 1960. A virus-induced epizootic hemorrhagic
disease of the Virginia white-tailed deer (Odocoileus virginianus). The Journal of
Experimental Medicine 111: 155-170.
Stallknecht, D. E., M. P. Luttrell, K. E. Smith, and V. F. Nettles. 1996. Hemorrhagic disease in
white-tailed deer in Texas: a case for enzootic stability. Journal of wildlife diseases 32:
695-700.
Stallknecht, D. E., M. L. Kellogg, J. L. Blue, and J. E. Pearson, 1991. Antibodies to bluetongue
and epizootic hemorrhagic disease viruses in a barrier island white-tailed deer
population. Journal of Wildlife Diseases 27: 668-674.
Stallknecht, D. E., A. B. Allison, A.W. Park, J. E. Phillips, V. H. Goekjian, V. F. Nettles, and J.
R. Fischer. 2015. Apparent increase of reported hemorrhagic disease in the midwestern
and northeastern USA. Journal of Wildlife Diseases 51: 348-361.
39

Stevens, G., B. McCluskey, A. King, E. O’Hearn, and G. Mayr. 2015. Review of the 2012
epizootic hemorrhagic disease outbreak in domestic ruminants in the United States. PloS
one 10: e0133359
Thomas, L., S. T. Buckland, E. A. Rexstad, J. L. Laake, S. Strindberg, S. L. Hedley, J. R. B.
Bishop, T. A. Marques, and K. P. Burnham. 2010. Distance software: design and analysis
of distance sampling surveys for estimating population size. Journal of Applied Ecology
47: 5–14.
Walter, W. D., D. M. Baasch, S. E. Hygnstrom, B. D. Trindle, A. J. Tyre, J. J. Millspaugh, C. J.
Frost, J. R. Boner, and K. C. VerCauteren. 2011. Space use of sympatric deer in a
riparian ecosystem in an area where chronic wasting disease is endemic. Wildlife
biology 17: 191-209.
Williams, R., and L. Thomas. 2007. Distribution and abundance of marine mammals in the
coastal waters of British Columbia, Canada. Journal of Cetacean Research and
Management 9: 15.
Wilson, A., and P. Mellor. 2008. Bluetongue in Europe: vectors, epidemiology and climate
change. Parasitology Research 103: S69-77.
Xu, B., M. Madden, D. E. Stallknecht, T. W. Hodler, and K. C. Parker. 2012. Spatial-temporal
model of hemorrhagic disease in white-tailed deer in south-east USA, 1983 to 2000.
Veterinary Record 170: 288.

40

CHAPTER 2: APPLICATION OF N-MIXTURE MODELS FOR AERIAL SURVEYS OF
WILDLIFE: A COMPARISON STUDY

INTRODUCTION
Aerial surveys are often the most practical methods for estimating wildlife population
abundance and have been developed to address issues of survey feasibility and reliability of
parameter estimates (Caughley 1974, Burnham et al. 1980). Accounting for imperfect detection
of abundance is important for obtaining precise and unbiased estimates (Denes et al. 2015).
However, methods to account for imperfect detection using aerial surveys are limited and not
always practical for wildlife managers. Distance sampling (Buckland et al. 2001) is a wellestablished method for obtaining abundance that allows accounting for imperfect detection, and
has been applied to aerial survey approaches numerous times in wildlife research. Unfortunately,
this method relies on assumptions that can be difficult to accommodate in aerial surveys.
Alternatively, N-mixture models (Royle 2004) are a new analytic tool that have recently been
applied to wildlife population estimation (Chandler and King 2011, Zipkin et al. 2014, Lyet et al.
2016, Keever et al. 2017). Similar to distance sampling, N-mixture models explicitly estimate
detection probability and abundance, yet have not been applied to aerial surveys of terrestrial
free-ranging wildlife to our knowledge. Evaluation of this emerging method is needed in the
context of traditional aerial survey methodologies.
Distance sampling has been one of the few methodological options that account for
imperfect detection while estimating abundance of unmarked, free-ranging wildlife populations.
This method has a long history of use for abundance estimation using count data and distance
measurements, and has proven applicable to a wide range of taxa and environments (Buckland et
41

al. 2004). Distance sampling originated from quadrat and line transect sampling in the early
1900’s and was developed further to utilize perpendicular distances (Buckland et al. 2000,
Buckland et al. 2001), ultimately resulting in a methodological extension of plot sampling where
not all animals are detected within the covered area (Buckland et al. 2010). The framework of
distance sampling (Buckland et al. 2001, 2004) typically does not present a fully model-based
analysis (i.e. not using a stochastic, non-randomized sample). Rather, the detection function is
derived from the observed data, relying on the assumption that detection of animals degrades as
distance from the transect line increases and these observations are a randomized sample of the
true population (Buckland et. al. 2001). Detection functions can then be calculated and are
informative for abundance or density estimates given the spatial extent of the study region. Using
this methodology assumes that the population of interest is closed, that distance measurements
are accurately collected, and that animals are distributed independently of a survey line or point
(Buckland et al. 2001).
Challenges of applying distance sampling in aerial surveys vary in degree and impact on
density estimates (White et al. 1989). Among the critical assumptions is that all animals at
distance zero are detected with certainty and that animals are detected and measured at their
initial location upon detection (Buckland et al. 2001).Violation of this particular assumption
results in biasing the density estimate low, which can easily occur in aerial surveys where
animals on the line-transect are missed when directly under the aircraft and not otherwise
counted or corrected for (Caughley 1974). Further, decisions about appropriate distance intervals
or cut-points and assignment of distance interval for observed animals may fail to meet the
assumption of accurate distance measurement for detected animals (Buckland et al 2001). Minor
shifts in altitude and angle affect distance measurements at ground level taken from a moving
42

aircraft. If an aerial survey is able to accommodate these assumptions, distance sampling can
effectively incorporate model covariates for estimating detection probability and is a useful
method for abundance estimation. However, procedures needed to meet these assumptions may
be challenging.
N-mixture models are a relatively new approach for estimating abundance and detection
processes in wildlife populations (Royle 2004). Like distance sampling, this analytical method
explicitly accounts for imperfect detection of sampled animals, but does so using a hierarchical
model for detection probability and abundance parameters (Royle 2004). This method appears to
be flexible for studies across taxa and for a wide range of ecological questions, allowing the
inclusion of habitat covariates, multispecies models, single species models, and estimation of
unmarked animals (Kery et al. 2005, Chandler and King 2011, Zipkin et al. 2014, Lyet et al.
2016, Keever et al. 2017). N-mixture models use count data to estimate population abundance
and associated detection probabilities for a closed population of interest. These models are a
hierarchical extension of a Poisson generalized linear model (GLM) and use temporal replicates
of counts that are sampled from multiple independent spatial locations for site-specific
abundance (Royle 2004). A principal contribution of this class of models is that the population
sizes, N, are specific to each spatial unit (site) and are independent random variables with a
mixing distribution. These models combine a binomial generalized linear model (GLM) and a
standard model (e.g., Poisson) for Ni (population abundance at site i). The observed counts at a
unit i are assumed to be independent realizations of a binomial process dependent on the
detection probability at a given site and time, 𝑝𝑖𝑗 (the observation process; Kery et al. 2005,
Royle 2004). Further, abundance at each site, Ni, has a Poisson or negative binomial distribution,
and is a random latent variable specified by parameters in the model as the ecological process
43

(Denes et al. 2015). The resulting estimated abundance parameters across sites can then generate
a posterior estimate of total abundance given some unit of space, which is determined by the
study design (Royle 2004). The spatial replication employed in these models and flexibility to
incorporate covariates has been shown to produce statistically and ecologically reasonable
estimates of abundance (Kery et al. 2005, Joseph et al. 2009).
While N-mixture models are increasingly applied to wildlife populations, they have
rarely been applied to aerial surveys for terrestrial wildlife. One study by Lyet and others (2016)
estimated abundance of Nile Crocodiles using aerial surveys along the Nile River and N-mixture
models. This method was useful for an aquatic obligate available for detection in a spatially
discrete water body (Lyet et al. 2016). This research may have been successful in part due to the
defined spatial area inherent in the study design for a water body. Spatial unit size selection for
N-mixture models is unclear for free-ranging animals with unknown spatial-use boundaries. No
previous study has used N-mixture models for aerial surveys of free-ranging terrestrial wildlife
and further evaluation is needed to understand the utility of this method for aerial surveys of
wildlife. The goal of my study is to evaluate whether N-mixture models could be applied to an
aerial survey of unmarked, free-ranging white-tailed deer (Odocoileus virginianus) for localscale population estimation. I seek to simultaneously apply distance sampling and N-mixture
models to a multi-year aerial survey of free-ranging deer to make direct comparisons of
parameter estimates between analytical methods. In this context, I am interested in designing and
specifying distance sampling methods for line-transect sampling approaches (Buckland et al.
2001, Buckland et al. 2004, Thomas et al. 2010), and adapting these line-transect surveys to
include data needed for N-mixture models. My specific objectives were to: 1) compare detection
and abundance parameters estimated from N-mixture models and distance sampling methods and
44

2) to test the sensitivity of the spatial unit size inherent in the study design to the N-mixture
approach.
STUDY AREA
The Maple River corridor in Montcalm, Ionia and Clinton, and Gratiot counties in south
central Michigan was the site of a major deer mortality event due to epizootic hemorrhagic
disease (MDNR) in 2012. My study area (190 km2) bisected this corridor and was generally
representative of adjacent landscapes (Fig. 1.1). Buck harvest, an index of county-level deer
population abundance prior to data collection, was 2,871, 2,614, 3,400, and 6,442 for Ionia,
Clinton, Gratiot, and Montcalm counties, respectively for 2012. This area is in the Lansing
loamy plain ecoregion characterized by a flat landscape supporting a mix of forest cover and
agriculture consisting primarily of soybeans, grain, and corn and fragmented woodlots. Greater
forest cover and water exist near the riparian corridor of the Maple River. The landcover
composition aggregated across our survey areas included 2.3 % developed land, 70 % forage
cropland, 9.7% forest, 1.1% shrubland, 16.3% wetland, and .6% bare ground or other landcover
(Fig. 1.1). The climate has large seasonal variation in temperatures ranging from an average low
of -10 degrees C to an average high of 28 degrees C, with 90.42 cm of average annual
precipitation, including snowfall. The Maple River area had notable EHD mortality reports in
2012 with Ionia county reporting 2,498 EHD deer mortalities and Clinton County reporting 514
EHD deer mortalities during the outbreak period (MDNR report 2014). A related research
project established independent deer abundance estimates each summer for this study area in the
same time period (Christensen 2017).

45

Figure 2.1 Aerial survey transects bisecting habitat in the Maple River Study Area of Michigan,
USA, 2014-2016.

METHODS
Aerial survey design
I designed an aerial survey protocol to collect data relevant to both distance sampling and
N-mixture model methods (Royle 2004). Survey flights were coordinated with Michigan
Department of Natural Resources (MDNR) pilots using a fixed-wing aircraft. Deer were
surveyed using a double-count approach where observers viewed deer from opposite sides of the
aircraft and then repeated the survey in the opposite direction, effectively counting both sides of
46

the transect line twice. We did this by using a line-transect method for distance sampling and
including a replicate survey in the opposite direction immediately following the first survey.
Observers recorded the location of deer or groups of deer along the transect using a GPS unit,
allowing the assignment of deer observations to spatial units (i) needed for analysis with Nmixture models. To standardize the effective survey area of estimation for both methods, we
selected a strip width of 900 m that corresponded to a significant decrease in detection
probability provided by the collected distance sampling data and graphed in program
DISTANCE (Thomas et al. 2010). This value was determined using distance data collected
during the flight and was specific to the species and conditions for this survey. Subsequently, the
900 m distance provided a boundary for all spatial units at a distance perpendicular to the
transect line. We were interested in population abundance estimates for the covered survey
region only and did not extrapolate to a larger study area using traditional distance sampling
approaches. We conducted a total of 6 winter surveys, 3 in 2014 and 3 in 2016, that provided
data for both analytical frameworks. Each of the 3 surveys in a given year were conducted >24
hours apart to ensure independence between surveys.
We used line-transect distance sampling to estimate deer abundance (Thomas et al.
2010). Transects were selected using a random starting point and were then assigned at
standardized distances from that point for a total of 5 north-south transects equally spaced at
4.8km apart and 21 km in length, covering our area of interest. The transect placement was
designed to exceed the average size of winter home range for white-tailed deer in a forestagriculture landscape to ensure we did not count the same deer on multiple transects or spatial
units (Hiller et al. 2009). Transects were placed over roads, which were effectively not surveyed
because they were not visible under the plane. I placed transects to bisect a major riparian
47

corridor running southwest to northeast within the study area to sample the deer population
across potential habitat-related population density gradients (Marques et al. 2013). Distance
intervals were demarcated on the window and strut of the aircraft and were calculated based on
guidelines for aerial surveys in Buckland et al. (2001) using relative interval cut-points for line
transects. These were 0-100m = A, 101-200m = B, 201-350m = C, 351-500m = D, and 501m1000m = E. To maintain accurate distance cut-points for each observer, we recalibrated the
appropriate angles for the desired intervals for each observer using a clinometer (Buckland et al.
2001). Surveys were consistently flown at near 153m altitude and a flight speed of 167 km/hr for
each survey (Caughley 1976), and surveys only occurred when full snow coverage existed for
maximum detection of deer.
The same 2 primary observers were maintained in the rear of the plane across years to
reduce variation across surveys by individual observer. Both primary observers were focused on
a single side of the aircraft when heading in a given direction during the survey. At the
completion of a survey pass, the aircraft immediately turned around and the transect was flown
in again in the opposite direction so each observer sampled both sides of the transect line once in
a survey day, allowing a replicated count in each day. No communication occurred between
observers on either side of the plane, but observers were able to communicate with 2 secondary
observers (pilot and forward passenger) who were sitting in front of each primary observer and
not responsible for collecting data aside from GPS coordinates when requested. Observers
collected the perpendicular distance category (distance categories A-E) for each deer
observation. Further, observers recorded covariate data of interest for each observation of deer
including: the initial observer (front or rear of aircraft), the group size, the activity of the deer at
time of observation (bedded, standing, or moving), and the habitat cover type associated with
48

each recorded observation (wooded, open, or residential). These covariates were included in our
analysis.
Distance sampling analysis
We used the conventional distance sampling and multiple covariate distance sampling
engines with a line-transect framework in Program Distance v. 6.2 as described by Thomas et al.
(2010) to estimate abundance and detection probability parameters. The standard formulation of
the population density is given with a function of abundance for a line transect:
̂

̂ = 𝑁 = 𝑛/2𝑤𝐿𝑃̂𝑎 ,
𝐷
𝐴
where constants include A = the size of the region of interest, n = the number of animals or
clusters detected, 𝐿 = ∑ 𝑙𝑗 where 𝑙𝑗 = the length of jth line and j = 1…k (a measure of survey
effort) where k is the number of line transects. The truncation distance for measurements from
̂ is the
the line is w, and 2w accounts for measurements on both sides of a transect line. 𝐷
̂ A).
estimated density of individuals per unit area (𝑁/

The distance sampling method relies on fitting a detection function from observed
distances to an object. Distance measurements were collected for each deer observation during
the survey, providing a frequency of detections across distances. The detection function,
parameter g(x), is estimated by fitting a model to our data. Selecting a key function to select a
good model of the true detection function (g(y)) will determine the basic model shape, reflecting
the data. This detection function contains the observation process, and assumes a degrading,
conditional relationship between detection probability and distance. The most commonly used
key functions have the best shape criterion and model robustness and were statistically derived
49

by Buckland and colleagues (2001). The half-normal and hazard rate are most commonly used,
and were evaluated for fit in all of our models. Key functions are important for estimating a
detection function given the data and should be evaluated with model selection and visual
inspection of distance frequency histograms (Buckland et al. 2001). Appropriate key functions
are robust, fit a variety of detection function shapes, and have a wide shoulder where detection
probability is nearly 1 close to the line transect with a precipitously declining detection
probability at some distance from the transect (Buckland 1992). To improve the estimation of the
detection probability using these functions, we truncated all data beyond 900 meters for both
2014 and 2016. This facilitated a comparable sampling area between years and across estimation
methods. We used Akaike Information Criterion (AIC) to select the most parsimonious model,
given available key functions and covariate terms, that best explained detection probability
(Burnham and Anderson 2002). We assessed the diagnostic histograms, qq-plots, and the percent
coefficient of variation for the most competitive model to ensure our data reasonably fit the
model with appropriate precision (Buckland et al. 2004). Once the best model was selected in
each year, I estimated abundance, detection probability, and 95% confidence intervals associated
with point-estimates.

N-mixture model analysis
As a class of state-space models, N-mixture models explicitly estimate the observation
process (the observed state) and the abundance process (the latent ecological state) assuming that
animals are detected and counted at spatial units (sites) with imperfect detection. Detection
probabilities are site-specific and derived from repeated counts at survey sites. Basic model
structure assumes that detections at a site are independent and that all individuals recorded at a
given site and time have the same detection probability, although detection can vary by
50

observation and with covariates (Royle 2004). We used a simple form of the N-mixture model to
establish parameter estimation and convergence for our data set, to then be compared to distance
sampling estimates using the same data. Because a replicate survey was conducted immediately,
we were able to assume population closure and establish parameter estimates for each survey
day. However, we were interested in population abundance and detection across survey days in a
given year, when deer may have moved into or out of our effective sampling area. To
accommodate variation at this time scale, we used an open population modification of the Nmixture model as described by Kery and Schaub (2012), and Dail and Madsen (2011). This
approach allowed us to consider each day as a primary occasion when populations could
fluctuate and the two replicate observations for each of the spatial units within a day as
secondary occasions when populations were assumed closed (Kery and Schaub 2012). Our
model used constant parameters for abundance and detection in each sampling period (day). The
observation model is shown as:
𝑦𝑖𝑗𝑘 | 𝑁𝑖 ~ 𝐵𝑖𝑛𝑜𝑚𝑖𝑎𝑙 (𝑁𝑖 , 𝑝𝑖,𝑘 ),
where 𝑦𝑖𝑗𝑘 is the total count for each sampling unit i (i = 1…x) during replicate j (j = 1, 2) and
day k (k = 1…3). Detections of deer are thus considered binomial random events and are not
dependent on individual identification (Royle and Dorazio 2008).

Mean abundance is estimated as the latent state and the count data are assumed to be
independent and identically distributed. The basic form of the abundance model is:
𝑁𝑖 ~ 𝑃𝑜𝑖𝑠𝑠𝑜𝑛 (𝜆𝑖 ),

51

where Ni is the true abundance at sampling location i. Although the same observation-level
covariates were collected for data used in the N-mixture model, we did not include covariates in
this model. We were interested in evaluating a simple form of the observation and process model
with count data alone, given the structure of our survey protocol. A Bayesian framework was
used to estimate parameters for the N-mixture model using software program JAGS (Plummer et
al. 2003) in program R (R team). We assumed vague prior distributions for estimated parameters
and evaluated posterior distributions, MCMC mixing plots, and Rhat for model convergence and
potential scale reduction factors (Gelman et al. 2006). We used a single, simple form of the
abundance and detection model for evaluation so did not rely on model selection as methods for
model selection in a Bayesian framework are lacking. Model code is available in Appendix A.

Sensitivity analysis of site size for N mixture models
One challenge posed when using N-mixture models is knowing the effective survey size of a
sampled area in continuous habitat for estimation of a species of interest (Kery and Royle 2015).
Some studies using the N-mixture model approach have focused on aquatic obligate species that
use spatially discrete water bodies, thus avoiding the problem (Zipkin et al. 2014, Lyet et al.
2016). Animal home range size can vary substantially given the behavior and movement of a
species, so assumptions must be made for the analysis or a previously estimated home range size
used as a boundary reference (Kery and Royle 2015). Distance sampling is exempt from these
problems and explicitly addresses the sampling area. This makes simple extrapolations or
estimates of animal density difficult to determine from abundance estimates with N-mixture
models, specifically. We resolved this issue in two ways. First, we used distance measurements
collected during the survey to demarcate an effective strip width for our transect survey. This
provided a boundary for distance sampling estimates (a standard truncation width) and for N52

mixture model estimates, where the boundary marked a single side of a rectangular-shaped
sampling unit for the latter. The unit side closest to the survey transect was equivalent to the “0”
measurement for the distance sampling line transect survey. Then, we selected the total length of
the spatial unit along the transect in the north-south direction, completing the sampling unit
rectangle. However, as this would have been an arbitrary value based on literature values for
deer movement, we developed a sensitivity analysis to evaluate how differing unit lengths (and
therefore total spatial unit size) would affect parameter estimates of abundance and detection.
We allowed the total number of sites i (spatial units) to vary as the length of each site varied, to
maintain full coverage of the survey transect area. Given the total length of the survey transect
was 21 km, spatial unit length scenarios included 1 km, 2 km, 4 km, 5 km, 6 km, 8 km, and 10
km, allowing a minimum of 2 units along a line transect side and a total of 20 spatial units across
the 5 transects. We used spatial coordinates for the start and end points of each survey transect to
further constrain the spatial units. Because deer observations were spatially referenced when
collected during the aerial survey, observations could be assigned to the appropriate spatial unit,
regardless of changing unit length. For each unit size scenario, I ran 1000 N-mixture models with
our collected replicated count data and mean parameter estimates.
RESULTS
I collected 172 observations of deer or deer groups in 2014 across 3 surveys, and
341observations across 3 surveys in 2016. Observers flew 5 transects in each survey day for a
total transect effort of 125.5 km for each transect, totaling 627.5 km of effort in each year. The
encounter rate for deer groups was 0.27 deer/km and 0.54 deer/km in 2014 and 2016,
respectively. Conditions were similar across surveys and years, as selected timing and snow
cover conditions allowed for clear visibility during all survey days. Observers surveyed deer
53

within the study area using a fixed-wing aircraft between 24 February and 4 March in both years.
Surveys were completed between 1100 and 1400 and generally took 2 hours to complete. All
surveys covered a total area of 189 km2, which was the transect length (21 km * 5 transects)
multiplied by 900 m *2, which was the truncated distance for analysis at 900m on both sides of
the transect.
Using distance sampling methods, the most parsimonious model selected by AIC was the
half-normal cosine function in 2014 and the hazard rate cosine function in 2016. Given our suite
of covariates, the model including observer effects was selected using AIC as the best model for
2014 (Table 2.1). This was not consistent for 2016, when the best model included habitat cover
effects rather than observer effects (Table 2.1). Despite differences in covariate inclusion across
years, estimates of the probability of detection were essentially the same for both years (p = 0.66;
CI: 0.6-0.73 for 2014; p = 0.67; CI: 0.62-0.73 for 2016; Table 2.2). The average abundance
estimate across all transects indicated an increase in deer population between years (215; CI:
129-359 for 2014; 425; CI: 277-652 for 2016; Fig. 2.2).
N-mixture models performed well with Rhat values below 1.05 and adequate mixing and
convergence of parameters with visual inspection (Gelman and Rubin 1992). We used a spatial
unit length of 2 km for all summarized posterior distribution parameter estimates, resulting in
100 spatial units total across the 5 transect area. The mean detection probability as determined
from the posterior distribution was higher in 2014 (0.48; CI: 0.443-53) than in 2016 (0.39; CI:
0.35-0.435), although this difference was not substantial (Table 3). Mean deer abundance
increased in 2016 from 2014, with approximately twice the deer population estimated in 2016
(699; CI: 630.67 - 787) than in 2014 (311; CI: 277.25 - 356.75, Figure 2.2).

54

900
800
700

Deer Abundance

600
2014 - Nmix

500

2014 - Distance
400

2016 - Nmix
2016 - Distance

300
200
100
0

2014

2016

Figure 2.2 A comparison of abundance estimates for white-tailed deer from distance sampling
and N-mixture models with associated 95% confidence intervals and 95% credible intervals,
respectively. Abundance represents deer within the Maple River study area of Michigan for 2014
and 2016.

55

The sensitivity analysis using real count data indicated that selected spatial unit size had a
significant effect on abundance and detection probability estimates, and their associated
variance. Generally, smaller unit sizes that matched home range sizes of white-tailed deer
performed well (2 km – 6km), although too small of a unit (1 km) greatly underestimated
detection probability and subsequently overestimated abundance in 2014 (Fig. 2.3 and 2.4).
Large spatial units consistently provided lower detection probabilities and increased standard
deviation in both years (Fig. 2.3-2.6). The 2 km unit size consistently performed the best across
years and among unit size scenarios, providing relatively high detection and low deviation, so
was chosen for our final parameter estimates (Fig. 2.3-2.6).
Table 2.1 Distance sampling models selected in each year and corresponding data descriptors.
Tot. Effort
Year

Model; Key function and # adj. terms

k

Tot. No. observations
(km)

2014

Observer; Half-normal, 0

2

172

627.5

2016

Cover; Hazard rate, 1

3

341

627.5

Table 2.2 Distance sampling parameter estimates and associated 95% confidence intervals from
best fitting models in each year for 2014 and 2016 in Michigan.
95%CI

Total area

Detection
Abundance

lower

upper

(km.2) of

Year

95%CI

lower

upper

limit

limit

Prob.
(N)

limit

limit

estimation

(P)

2014

189

215

129

359

0.66

0.6

0.73

2016

189

425

277

652

0.67

0.62

0.73

56

Table 2.3 Summarized posterior distributions for N-mixture model parameter estimates, using
the 2 km spatial unit, in each year for 2014 and 2016 in Michigan.
Credible interval quartiles for each parameter
Year

Parameter

Mean

SD

2.50%

25%

50%

75%

97.50%

Rhat

2014

Abundance (N)

311.01

20.14

277.25

297

309

323

356.75

1

2014

Detection prob. (p)

0.49

0.02

0.44

0.47

0.48

0.5

0.53

1

2016

Abundance (N)

699.40

39.69

630.67

672

695.67

723.33

787

1

2016

Detection prob. (p)

0.39

0.02

0.349

0.38

0.39

0.41

0.44

1

57

a)

b)

Figure 2.3 Sensitivity of 2014 N-mixture models to spatial unit size for a) population abundance
(N) and b) associated standard deviation (SD) of abundance in Michigan.

a)

b)

Figure 2.4 Sensitivity of 2014 N-mixture models to spatial unit size for a) detection probability
(p) and b) associated standard deviation (SD) of detection probability.

58

a)

b)

Figure 2.5 Sensitivity of 2016 N-mixture models to spatial unit size for a) population abundance
(N) and b) associated standard deviation (SD) of abundance.
a)

b)

Figure 2.6 Sensitivity of 2014 N-mixture models to spatial unit size for a) detection probability
(p) and b) associated standard deviation (SD) of detection probability in Michigan

59

DISCUSSION
We found that application of distance sampling and N-mixture models both produced
realistic abundance estimates for white-tailed deer in our study. Distance sampling methods were
used according to specifications for aerial application that have been well developed in the
literature (Buckland et al. 2001, Fewster and Pople 2008). Aerial survey protocols were easily
adjusted for the application of N-mixture models. N-mixture models were effective in estimating
abundance and appear to be a practical tool for estimating abundance of a free-ranging wildlife
population. Specifically, N-mixture models offered good precision in abundance parameter
estimates relative to distance sampling in this case, and met model assumptions. However, we
found that estimates from N-mixture models were sensitive to the size of the survey unit
selected. While distance sampling has been applied to aerial surveys in the past, to our
knowledge, this is the first application of an N-mixture model using an aerial survey for a
terrestrial wildlife species. N-mixture models appear to be robust and easily applied to aerial
surveys for many species.
Using the distance sampling approach, we estimated detection probability and abundance
for deer, providing an independent analysis using a recognized population estimation method.
Detection probability estimates were consistently greater than N-mixture model detection
probability estimates, although precision was similar to N-mixture models. This parameter may
have benefited from the inclusion of covariate data to explain variation in deer observations,
although competing models in program DISTANCE with no included covariates had similar
detection probability estimates. Abundance estimates were reasonable for our study although
they did not demonstrate an increase in population in 2016 as expected.

60

Distance sampling has been an important tool for wildlife professionals to assess
abundance while assessing imperfect detection and has been used with aerial surveys many times
over. This is one of the few methods available that explicitly defines the limits of a study area,
which is a significant challenge in other estimation methods (Kery and Royle 2015). However,
the assumptions of distance sampling may be difficult to meet given environmental or design
constraints for some species. For example, if an animal moves in response to a low-flying
aircraft and an accurate distance is not recorded at the animal’s original location, or if the cover
is such that very few detections can be made, reliable estimates cannot be generated. Further, if
the aircraft obstructs the transect line, the assumption of perfect detection at the line (distance 0)
is violated. This is difficult to address with survey methods, and post-hoc analysis are often used
to correct for imperfect detection or surveillance at the transect line. Another important
assumption is that animals are statistically independent, are arranged under some stochastic
process in the environment, and are distributed independently of a line or point (Buckland et al.
2001). In many landscapes, the distribution of roads, private land, or geophysical barriers make a
fully randomized study design difficult, violating this assumption. Finally, aerial surveys often
rely on binned distance data from demarcated distance intervals on the aircraft. This presents two
challenges, the first being the accuracy of the distance bin assignment for the animal observation
at high speeds and varying aircraft yaw. The second occurs in the analytical process when fitting
a detection function to the observed distances, as too few bins may result in a poor model fit.
Modern field equipment such as range finders may allow improve accuracy, although data
should be examined for rounding or heaping issues (where observers record measurements to
convenient cut-offs, like 50, 100, 200). If estimating these distance intervals for aerial surveys at
a known altitude, each observer must be individually calibrated for sight angle with a clinometer

61

to ensure standard distances despite observer height. Distance sampling can be accomplished in
single or multiple sampling visits with counts and distance measurements, which makes it easy
to implement and relatively cost effective, if the above mentioned assumptions are met or easily
dealt with (Denes et al. 2015).

The N-mixture model provided reasonable abundance estimates and detection
probabilities (Kellner and Swihart 2014). The most significant difference with distance sampling
results was in the precision of estimated abundance, with the N-mixture model providing a more
precise estimate. We show that N-mixture models were effectively applied to survey a
population of deer that increased in abundance over time. This detected population increase was
supported by concurrent research from another independent survey in the same study area
(Christensen 2017). Population abundance is more often the parameter of interest for wildlife
managers and was similar by year for each analytical method. In this case, parameter precision
was important as it suggested a significant population increase using the N-mixture model, and
no population increase with the distance sampling model. Both modeling approaches tracked an
increased mean abundance in 2016, but the interpretation and potential management implications
thus varied. The precision afforded by this method allowed a clear interpretation of population
change.
We found the abundance and detection parameters estimated with the N-mixture model
approach were highly sensitive to the spatial unit size selected in the study design. This
addressed an important assumption for N-mixture models. The effective area and shape of a
surveyed population in an N-mixture modeling framework is usually not known and subjectively
selected by the researcher, revealing a significant source of uncertainty when implementing and
interpreting population parameters for any species (Kery and Royle 2015). While some
62

assumptions can be made about the area in question or the availability of animals within that area
in a continuous habitat, we found that directly assessing spatial unit size associated with
abundance and detection parameters for deer was critical for our analysis. Ultimately, researchers
or managers need to define a population of interest and subsequent survey area the area in some
way. It is important to address the spatial extent of a spatial unit for this method, and ideally
match these units to ecologically relevant space-use by the animal in study. Evaluation of 7
different spatial unit scenarios indicated that the 2 km spatial selection was ideal. This
corresponds closely with a reported value of 140 ha by Hiller et al. (2009) who estimated spaceuse during the winter for white-tailed deer in south-central Michigan. In this regard, distance
sampling is explicit and may have an advantage over N-mixture models. Regarding our distance
sampling estimates, we were not interested in extrapolating our estimates to a larger study area
(which is uncommon for distance sampling), reducing the likelihood of biasing estimates with
our selected study design
Aerial surveys are a pillar of wildlife population assessment and have been developed
and applied in countless wildlife studies. The advancement of N-mixture models provides an
exciting opportunity to improve population estimates using aerial surveys. N-mixture models are
less ambiguous when modeling the imperfect detection known to be present in sampling count
data and are considered a model-based approach. Detection is not predicated on a single
relationship with distance and comes from some common distribution allowing good precision
and a clear understanding of relationships with other covariates, if included (Kery et al. 2005).
This technique may be a better path forward when faced with uncertainty of fundamental
assumptions with traditional methods.

63

N-mixture models have advantages over traditional distance sampling methods if there is
interest in reducing the need for auxiliary data collection (i.e. distance measurements) or other
challenges associated with the assumptions required for reliable estimation with distance
sampling methods exist. In our case, these models provided a practical field approach and
analytical methodology that resulted in precise estimation of unmarked deer. We suggest more
applications of these models are needed with field data across taxa and attention must be given to
the spatial unit associated with the species’ space use. While our current study design required an
ad hoc approach to sub-sample collected count data in a manner consistent with N-mixture
models, it may provide a good example of application for a species of interest with well-studied
population parameters. Recent developments have in fact merged distance sampling and Nmixture binomial modeling approaches as hierarchical distance sampling, which may be
including the best of both model classes (Sollmann et al. 2015). These newly developed models
use the Dail and Madsen (2011) process model with the standard distance sampling (Buckland et
al. 2001) observation model. N-mixture models have been demonstrated to be applicable to zero
inflation data and continued to be developed for open and dynamic, stage-structured populations
(Dail and Madsen 2011, Zipkin et al. 2014, Denes et al. 2015). This is clearly an area of research
interest and current development and application to readily available survey methods, such as
aerial surveys, will be important.

64

APPENDICES

65

APPENDIX A: R CODE FOR N-MIXTURE MODEL LISTED IN CHAPTER 2
library(R2jags)
library(rjags)
load("Aerial2016sims8000m.RData") #load in Rdata sim.dat object
dimen <- dim(sim.dat)
output <- data.frame(Iteration=integer(),
N=numeric(),
sdN=numeric(),
p=numeric(),
sdp=numeric(),
stringsAsFactors=FALSE)
for (m in 1:dimen[3]){ #dimen[3]){
R <-max(sim.dat[,1,m]) #number of sites: read from sim.dat because it will change for different
size units
t <-2 #rep. counts (always stays the same)
d <-max(sim.dat[,2,m]) #days: read from sim.dat because it will change for different years
# Create structure to contain counts
y <- array(dim=c(R,t,d))
for(l in 1:d){
sel.rows <- sim.dat[,2,m] == l
y[,,l] <- sim.dat[sel.rows,3:4,m]
}

# Specify model in BUGSlanguage
modelstring = "
model {
# Priors
for (k in 1:days) {
alpha.lam[k] ~dnorm(0,0.01)
66

p[k] ~ dunif(0,1)
}

# Likelihood

# Biological model for true abundance
for (k in 1:days) {

#loop over d days

lambda[k] <- exp(alpha.lam[k])
for (i in 1:R) {

#loop over R sites

N[i,k] ~ dpois(lambda[k]) #abundance

# Observation model for replicated counts
for (j in 1:t){
y[i,j,k] ~ dbin (p[k], N[i,k]) #detection, looped over 2 secondary reps
}
}
}
# Derived and other quantities
for (k in 1:days) {
totalN[k] <- sum(N[,k]) # Total pop size across all sites
mean.abundance[k] <- exp(alpha.lam[k])
}
mean.N <- mean(totalN[])
mean.p <- mean(p[])
}
sink()
tmpf=tempfile()
tmps=file(tmpf,"w")
cat(modelstring,file=tmps)
close(tmps)
67

# Bundle data
# code for single day
win.data <-list(y=y, R=nrow(y), t=ncol(y), days=d) #y="C" in AHM book

# Initial values
Nst <- apply(y, c(1,3), max) + 1
#Nst[is.na(Nst)] <- 1
inits <-function(){list(N=Nst, alpha.lam = runif(d,-1,1))}

# Parameters monitored
params <-c ("mean.N", "mean.p")

# MCMCsettings
ni <-10000
nt <-2
nb <-1000
nc <-3

# CallWinBUGSfromR(BRT0.1min)
#out <-jags(win.data,inits,params,"modelcount2.26.2015.txt",n.chains=nc,
# n.thin =nt,n.iter=ni,n.burnin=nb,debug=TRUE, working.directory=getwd())
#possible using rjags
#model=jags.model(tmpf,data=win.data, inits=inits, n.chains = nc)
#update(model,n.iter=nb) #may be incorrect
#output <- coda.samples(model = model,variable.names = params, n.iter = ni, thin=nt)

out16.8test <- jags(win.data, inits, params, model = tmpf, n.chains = nc, n.thin = nt,
n.iter = ni, n.burnin = nb)
print(out14, dig = 2)
68

#jagsout <- as.mcmc.list(out14$BUGSoutput)
#plot(jagsout)
#summary(jagsout)

#extract mean and sd for mean.N and mean.p
output[m,1] <- m
output[m,2] <- out16.8test$BUGSoutput$mean$mean.N[1]
output[m,3] <- out16.8test$BUGSoutput$sd$mean.N[1]
output[m,4] <- out16.8test$BUGSoutput$mean$mean.p[1]
output[m,5] <- out16.8test$BUGSoutput$sd$mean.p[1]

rm(out16.8test)
print(m)
}
write.csv(output,"output16.8ktest.csv")

69

APPENDIX B: SENSITIVITY ANALYSIS CODE (2KM) FOR CHAPTER 2

transect.endpoints <- read.csv("StartPoint.csv")
aerial.dat <- read.csv("Aerial2016xy.csv")
nrep = 1000 #some number of simulations for each unit assignment
unit.size = 2000 #define unit size (km)
n.unit <- floor(21000/unit.size) #calculate number of units of size unit.size that will fit in a
transect
len.units <- n.unit * unit.size #total transect dist contained by units
extra.transect <- 21000 - len.units #remaining transect length (not in units)
spacing <- floor(extra.transect/n.unit) #calculate length of breaks between units
#start.point <- floor(runif(1, 0, spacing)) #moved into loop

unitrows <- n.unit * 2 * 5 #determine total number of units for the study site N,S and 5 transects
vdays <- unique(aerial.dat$Date) #use the data to determine the number of survey days
ndays <- length(vdays)

sim.dat <- array(NA, c(unitrows*ndays,4,nrep)) #create an 3D array full of NAs to hold the data
for all the iterations
sim.dat[,1,] <- rep(1:unitrows,ndays) #fill the unit numbers in the data array

vector.length <- n.unit * 2 #determine length of vector that will hold break points for units
for (h in 1:nrep){
start.point <- floor(runif(1, 0, spacing))
for (m in 1:ndays){
#loop through i study areas (we have 2 so I've hard coded that number)
i = 1 #for (i in 1:2){ #maybe unneeded because each year is only one site
#loop through j transects (we have 5 so I've hard coded that number)
for (j in 1:5){
dat <- subset(transect.endpoints, StudyArea==i & Transect==j & Start_Sout==1) #code for
study site
70

ycoord <- dat$y #read in Y coord for beginning of transect i
breaks <- vector(mode = "integer",vector.length) #create vector to hold start, breaks, and end
(full of zeros)
breaks[1] <- ycoord + start.point
for (k in seq(2,vector.length-2,2)){
breaks[k] <- breaks[k-1] + unit.size
breaks[k+1] <- breaks[k] + spacing
}
breaks[vector.length] <- breaks[vector.length-1] + unit.size

#fill in data table (unit, day, count1, count2) based on breakpoints
aerial.sub <- subset(aerial.dat, Transect==j)
north.sub1 <- subset(aerial.sub, Heading_Di == "N" & OBS == 1 & Date == vdays[m])
north.sub2 <- subset(aerial.sub, Heading_Di == "N" & OBS == 2 & Date == vdays[m])
south.sub1 <- subset(aerial.sub, Heading_Di == "S" & OBS == 1 & Date == vdays[m])
south.sub2 <- subset(aerial.sub, Heading_Di == "S" & OBS == 2 & Date == vdays[m])

for (l in 1:(n.unit * 2)){
if (l < (n.unit+1)){
C1sub <- subset(north.sub2, Y >= breaks[(l*2)-1] & Y <= breaks[l*2])
C2sub <- subset(south.sub1, Y >= breaks[(l*2)-1] & Y <= breaks[l*2])
}else{
C1sub <- subset(south.sub2, Y >= breaks[((l-n.unit)*2)-1] & Y <= breaks[(l-n.unit)*2])
C2sub <- subset(north.sub1, Y >= breaks[((l-n.unit)*2)-1] & Y <= breaks[(l-n.unit)*2])
}
sim.dat[l + (n.unit * 2 * (j-1)) + ((n.unit*2*5) * (m-1)), 3, h] <- sum(C1sub$Group_Size)
sim.dat[l + (n.unit * 2 * (j-1)) + ((n.unit*2*5) * (m-1)), 4, h] <- sum(C2sub$Group_Size)
sim.dat[l + (n.unit * 2 * (j-1)) + ((n.unit*2*5) * (m-1)), 2, h] <- vdays[m]
}
}
71

}
print(h) #prints iteration to see progress
}
save(sim.dat,file = "Aerial2016sims2000m.RData")
sim.dat[,,2]

72

LITERATURE CITED

73

LITERATURE CITED

Buckland, S. T. 1992. Fitting density functions with polynomials. Applied Statistics 41: 63-76.
Buckland, S. T., I. B. J. Goudie, and D. L. Borchers. 2000. Wildlife population assessment: past
developments and future directions. Biometrics 56: 1-12.
Buckland, S. T., D. R. Anderson, K. P. Burnham, J. L. Laake, D. L. Borchers, and L. Thomas.
2001. Introduction to distance sampling: estimating abundance of biological populations.
Oxford University Press, Oxford, United Kingdom.
Buckland, S. T., D. R. Anderson, K. P. Burnham, J. L. Laake, D. L. Borchers, and L. Thomas.
2004. Advanced distance sampling. Oxford University Press, Oxford, United Kingdom.
Buckland, S. T., A. J. Plumptre, L. Thomas, and E. A. Rexstad. 2010. Design and analysis of line
transect surveys for primates. International Journal of Primatology 31: 833-847.
Burnham, K. P., D. R. Anderson, and J. L. Laake. 1980. Estimation of density from line transect
sampling of biological populations. Wildlife monographs 72: 3-202.
Burnham, K. P., and D. R. Anderson. 2002. Model selection and inference: a practical
information-theoretic approach. Springer Publishing, Second edition, New York, New
York, USA.
Caughley, G. 1974. Bias in aerial survey. The Journal of Wildlife Management 38: 921-933.
Caughley, G., R. Sinclair, and D. Scott-Kemmis. 1976. Experiments in aerial survey. The Journal
of Wildlife Management 40: 290-300.
Chandler, R. B. and D. I. King. 2011. Habitat quality and habitat selection of golden‐winged
warblers in Costa Rica: an application of hierarchical models for open
populations. Journal of Applied Ecology 48: 1038-1047.
Christensen, S.A. 2017. Epizootic hemorrhagic disease affects abundance and distribution of
free-ranging deer populations. Doctoral dissertation, chapter 1.
Dail, D., and L. Madsen. 2011. Models for estimating abundance from repeated counts of an
open metapopulation. Biometrics 67: 577-587.
Denes, F. V., L. F. Silveira, and S. R., Beissinger. 2015. Estimating abundance of unmarked
animal populations: accounting for imperfect detection and other sources of zero
inflation. Methods in Ecology and Evolution 6: 543-556.

74

Fewster, R. M. and A. R., Pople. 2008. A comparison of mark–recapture distance-sampling
methods applied to aerial surveys of eastern grey kangaroos. Wildlife Research 35: 320330.
Gelman, A. 2006. Prior distributions for variance parameters in hierarchical models. Bayesian
Analysis. 1: 515–534.
Gelman, A., and D. B. Rubin. 1992. Inference from iterative simulation using multiple
sequences. Statistical science 7: 457-472.
Hiller, T. L., H. Campa III, and S. R. Winterstein, 2009. Estimation and implications of space
use for white-tailed deer management in southern Michigan. Journal of Wildlife
Management 73: 201-209.
Joseph, L. N., C. Elkin, T. G. Martin, and H. P. Possingham. 2009. Modeling abundance using
N‐mixture models: the importance of considering ecological mechanisms. Ecological
Applications 19: 631-642.
Keever, A. C., C. P. McGowan, S. S. Ditchkoff, P. K. Acker, J. B. Grand, and C. H. Newbolt.
2017. Efficacy of N-mixture models for surveying and monitoring white-tailed deer
populations. Mammal Research 62: 413-422.
Kellner, K. F., and R. K. Swihart. 2014. Accounting for imperfect detection in ecology: a
quantitative review. PLoS One, 9: e111436.
KĂŠry, M., J. A. Royle, and H. Schmid. 2005. Modeling avian abundance from replicated counts
using binomial mixture models. Ecological Applications 15: 450-1461.
KĂŠry, M., and M. Schaub. 2011. Bayesian population analysis using WinBUGS: a hierarchical
perspective. Academic Press.
KĂŠry, M., and J. A. Royle. 2015. Applied Hierarchical Modeling in Ecology: Analysis of
distribution, abundance and species richness in R and BUGS: Volume 1: Prelude and
Static Models. Academic Press.
Lyet, A., R. Slabbert, W. F. Versfeld, A. J. Leslie, P. C. Beytell, and P. D. Preez. 2016. Using a
binomial mixture model and aerial counts for an accurate estimate of Nile crocodile
abundance and population size in the Kunene River, Namibia. South African Journal of
Wildlife Research 46: 71-86.
Marques, T.A., S. T. Buckland, R. Bispo, and B. Howland. 2013. Accounting for animal density
gradients using independent information in distance sampling surveys. Statistical
Methods & Applications 22: 1-14.
MDNR Report June 2014 Report #3585 http://www.michigan.gov/dnr/0,4570,7-15310363_48664---,00.html (Accessed September 17 2014).
75

Plummer, M. 2003. March. JAGS: A program for analysis of Bayesian graphical models using
Gibbs sampling. In Proceedings of the 3rd international workshop on distributed
statistical computing 124: 125.
Royle, J. A. 2004. N‐mixture models for estimating population size from spatially replicated
counts. Biometrics 60: 108-115.
Royle, J. A., and R. M. Dorazio. 2008. Hierarchical modeling and inference in ecology: the
analysis of data from populations, metapopulations and communities. Academic Press.
R Development Core Team. 2011. R: A language and environment for statistical computing. R
Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL
http://www.R-project.org.
Sollmann, R., B. Gardner, R. B. Chandler, J. A. Royle, and T. S. Sillett. 2015. An open‐
population hierarchical distance sampling model. Ecology 96: 325-331.
Thomas, L., S. T. Buckland, E. A. Rexstad, J. L. Laake, S. Strindberg, S. L. Hedley, J. R. B.
Bishop, T. A. Marques, and K. P. Burnham. 2010. Distance software: design and analysis
of distance sampling surveys for estimating population size. Journal of Applied Ecology
47: 5-14.
White, G.C., R. M. Bartmann, L. H. Carpenter, and R. A Garrott. 1989. Evaluation of aerial line
transects for estimating mule deer densities. The Journal of Wildlife Management 53:
625-635.
Zipkin, E. F., J. T. Thorson, K. See, H. J. Lynch, E. H. C. Grant, Y. Kanno, R. B. Chandler, B.
H. Letcher, and J. A. Royle. 2014. Modeling structured population dynamics using data
from unmarked individuals. Ecology 95: 22-29.

76

CHAPTER 3: THE ROLE OF DROUGHT AS A PREDICTOR OF EMERGENT
HEMORRHAGIC DISEASE IN THE EASTERN UNITED STATES.

INTRODUCTION
Vector-borne diseases are responsible for nearly 30% of all emerging infectious disease
events in the last decade (Jones et al. 2008). The association between climate change and the
ecology of vector populations is important for understanding disease risk as shifting climatic
conditions are often associated with increased disease emergence (Guis et al. 2011, Altizer et al.
2013, Paull et al. 2017). Changing climate patterns will affect vector-borne diseases
disproportionately more than other disease types, as vector populations and pathogen
transmission are sensitive to climatic and environmental factors (Jones et al. 2008, Harvell et al.
2002, Guis et al. 2011). This impact is particularly true in regions where cooler temperatures
have traditionally been the limiting threshold for a vector community (Rogers and Randolph
2000). However, the epidemiological dynamics of vector-borne disease are often complicated by
other factors, including host dynamics and physical characteristics of the environment, making
climatic relationships difficult to evaluate directly (Altizer et al. 2013, Koelle et al. 2005).
Hemorrhagic disease (HD), which includes serogroups in the genus Orbivirus for both
bluetongue virus (BTV) and epizootic hemorrhagic disease virus (EHDV), is a globally
distributed vector-borne disease affecting domestic and wild ruminant species (Nettles and
Stallknecht 1992, Mellor and Boorman 1995). The vectors are biting midges in the genus
Culicoides, which vary in their species distribution, vector competence, and subsequent roles in
HD transmission (Gibbs and Greiner 1989, Howerth et al. 2001, Ruder et al. 2015a). In Europe,

77

significant outbreaks of blue tongue (BT) have caused high mortality in sheep and cattle (> 1
million) as the disease shifted further north and became emergent in the region (Purse et al.
2005). These outbreaks of BT have been linked to climate change and shifting wind patterns
affecting Culicoides distribution and disease transmission (Wittmann and Baylis 2000, Purse et
al. 2005, Wilson and Mellor 2008, Hartemink et al. 2009, Guis et al. 2012). In North America,
BT does affect some cattle and sheep annually, but has not resulted in significant domestic
animal mortality events, as seen in Europe. Instead, BTV and EHDV serotypes are often the
source of acute mortality for wild ruminants, particularly white-tailed deer, although EHDV
serotypes are responsible for the majority of clinical sign in these populations (Shope et al. 1960,
Ruder et al. 2015a). Because these virus serogroups result in clinically indistinguishable disease
in deer, epizootic hemorrhagic disease (EHD) and BT are referred to jointly as HD, which is
considered the most significant viral disease for white-tailed deer in North America (Thomas et
al. 1974, Nettles and Stallknecht 1992).
Hemorrhagic disease had previously been restricted to the southeastern United States
(US) and portions of the Great Plains but recently has expanded north in range and increased in
frequency across the Midwest and northeastern US (Stallknecht et al. 2015). The increased
frequency and spatial distribution of EHD-related mortality in deer populations at northern
latitudes is likely associated with processes affecting vector transmission and host characteristics.
Vector transmission and the ecology of BTV- and EHDV-competent Culicoides species are
sensitive to climatic and environmental conditions, leading researchers to suspect a role of
climate and midge habitat on the increase in HD reports for the US (Berry et al. 2013, Ruder et
al. 2015a, Ruder et al. 2015b, Stallknecht et al. 2015). The majority of disease in white-tailed
deer occurs seasonally when the adult midge transmits the virus after taking a blood meal from

78

the deer during the summer and early autumn (Mellor et al. 2000). The midge then completes its
reproductive cycle and lays eggs in a wet soil substrate, typically near rivers and streams
(Wilkening et al. 1985, Purse et al. 2005). Accordingly, wetland features tend to support
abundant vector populations in late summer when transmission of the virus to deer occurs
(Savini et al. 2011). Further, replication of EHDV serotypes 1, 2, and 7 in a confirmed HD vector
(C. sonorensis) has been shown to increase with increasing temperature (Ruder et al. 2015b).
Warm temperatures in association with high soil moisture appear to be ideal conditions for both
vector and pathogen persistence, although limited data exists on confirmed vector species and
transmission rates where HD is emergent (Ruder et al. 2015a). The link between the life-history
requirements of the vector and disease occurrence at large scales was supported by correlative
research that found long-term (~25 year) increases in wetland cover at a county level helped
explain increases in HD morbidity reports across the continental US (Berry et al. 2013).
However, wetland features are sensitive to extreme temperatures and changes in precipitation,
such as severe drought conditions, suggesting some interaction between wetlands and climatic
factors in the vector-host-pathogen dynamic.
Few studies have investigated the role of climate on the spatial and temporal patterns of
hemorrhagic disease in the US. Sleeman et al. (2009) found that higher than average winter and
summer temperatures, and lower than average spring precipitation, explained annual HD
incidence over 5 years in Virginia, where HD is enzootic. In a five-state spatiotemporal analysis,
Xu and others (2012) found wind speed, dew point, land surface temperature, rainfall, and NDVI
(a measure of vegetative productivity) were important factors explaining HD occurrence in the
enzootic Southeast between1980 and 2000. A recent study by Stallknecht and others (2015)
found that the number of counties reporting HD had increased between 1980 and 2012 and had

79

expanded northerly. Further, this study assessed drought severity temporally, using US Drought
Monitor (DM) data, which is a widely available indicator of drought intensity incorporating soil
moisture and temperature (Svoboda et al. 2002, Stallknecht et al. 2015a). The cumulative DM
intensity data for each county and year were positively correlated with years of greater HD
reports across the eastern US (Stallknecht et al. 2015). Drought conditions are inherently climate
related and commonly considered to occur when environmental conditions are drier than average
resulting from warmer than average temperature and lower than average precipitation for some
time period. Previous research has alluded to the role of drought on HD but to our knowledge, no
study has directly evaluated the spatiotemporal role of drought severity on emerging HD
patterns.
The observed heterogeneity of HD mortality reports across the range of the disease and
the ability to predict patterns of disease is further complicated by herd immunity dynamics in
white-tailed deer (Howerth et al. 2001, Gaydos et al. 2002b). Previous research has indicated that
HD meets the criteria for enzootic stability, which occurs when disease-incidence is maximized
at moderate levels of disease transmission and seroprevalence within a host population, resulting
in high levels of disease infection and low clinical sign (Coleman et al. 2001, Howerth et al.
2001, Stallknecht et al. 1996, Park et al. 2013). The degree of stability refers to the level of
infection from a specific pathogen among individuals within a population (Coleman et al. 2001).
In the case of HD, acquired and innate immune responses exist and are typically life-long for
white-tailed deer (Stallknecht et al. 1991, Gaydos et al. 2002a, Gaydos et al. 2002b).
Seroprevalence data and long-term HD datasets have demonstrated endemic stability for HD in
the southern regions of the US, with a reduction in seroprevalence and an increase in HD
mortality along a north-south latitudinal gradient (Stallknecht et al. 1996, Stallknecht et al.

80

2015). Using seroprevalence data from 16 states, research by Park and colleagues (2013) found
that intermediate levels of infection in a population (~50%) resulted in the greatest number of
symptomatic disease cases and mortality reports, which is contrary to commonly predicted
positive linear correlations between pathogen transmission and disease incidence (Park et al.
2013). Over time, naĂŻve deer populations will recover as deer are born or immigrate into an
affected area, but population recovery and risk for future mortality will be associated with the
maintenance of seroprevalence within the population (Park et al. 2013).
Given the vector-borne nature of HD and its increased frequency and intensity in
emergent zones, there is a research need to investigate the role of drought as a predictor of HD
across the gradient of epidemiological zones of known disease occurrence. Our understanding of
the epidemiology of HD is incomplete especially related to identifying risk factors associated
with emerging HD outbreaks. I sought to evaluate the role of drought severity on changes in HD
reports across the eastern US for a 15-year period when the disease has increased its northern
range. My objectives were to 1) develop a spatiotemporal model to evaluate if drought severity
and related environmental variables explain changing patterns of HD presence; and 2) to
determine if this potential risk factor varies in importance over the present range of HD in the
eastern United States. I did this by using spatially-explicit statistical models to evaluate drought
severity and related environmental factors (wetland cover and physiographic region) on HD
mortality patterns from a longitudinal dataset of presence-absence county-based HD surveillance
data. I use remotely sensed data and drought monitor data from a period of known HD reporting
increase (2000-2014) to characterize explanatory variables for a spatiotemporal model to address
variation in disease patterns over space and time. This approach provides a greatly improved

81

understanding of the role of drought as a predictor of HD outbreaks across the epidemiological
range of the disease.
STUDY AREA
My research focused on counties within 23 states spanning the majority of the eastern
US, excluding New England, New York, and Minnesota where HD has not been found prior
during the time period of interest. These states included Alabama, Arkansas, Delaware, Florida,
Georgia, Illinois, Indiana, Iowa, Kentucky, Louisiana, Maryland, Michigan, Mississippi,
Missouri, New Jersey, North Carolina, Ohio, Pennsylvania, South Carolina, Tennessee, Virginia,
West Virginia, and Wisconsin. All states had previously reported HD at least once between the
years 2000-2014 and represent the epidemiological gradient from enzootic in the southeastern
US to emergent and incursive in the upper Midwest and Northeast (Stallknecht et al. 2015). This
was demonstrated in the variation of HD occurrence and frequency of reports among the states
selected for study (Stallknecht et al. 2015). All states examined in our study were within the
known range for white-tailed deer (Hewitt 2011).
METHODS
Data Collection
I acquired data from a longitudinal HD surveillance survey conducted by the
Southeastern Cooperative Wildlife Disease Study (SCWDS), at the College of Veterinary
Medicine, University of Georgia (Nettles et al. 1992, Stallknecht and Howerth 2004). Surveys
have been conducted annually since 1982 for HD status within each county for the entire US by
state agency deer biologists and reported to SCWDS. I chose to evaluate 15 years, spanning from
2000-2014, when HD outbreaks precipitously increased in frequency and range across the
82

eastern US (Stallknecht et al. 2015). Four categories of HD reporting criteria are included in the
long-term survey, although I confined the study to the three categories associated with direct
deer mortality: 1) sudden and unexplained deer mortality that occurs in late summer and early
fall; 2) necropsy diagnosis of HD as determined by a trained professional; and 3) isolation or
molecular-based detection of EHDV or BTV. If a county reported any HD mortality in a given
year, that county was identified as having HD present (1) in that year, rather than absent (0).
Despite existing count data for deer mortalities in many HD reporting counties, I chose a
conservative presence-absence approach for disease reporting to reduce uncertainty in reporting
bias across states in the annual surveys. In total, 23 states were included, all of which had
reported cases of HD within the 15-year time period at the state level. All counties (1,830) within
each state were subsequently included in the analysis and all covariate data were estimated at the
county level. To account for spatial autocorrelation at the county level, I calculated the spatial
centroid for each county using a geographic information system (GIS) and used the resulting
latitude and longitude as explanatory covariates (Berry et al. 2013; ArcMap, Environmental
Systems Research Institute. Redlands, California, USA).
I calculated drought values for each year and county using the DM tabular data (Svoboda
et al. 2002). Archived DM data indicate the percent area of each DM intensity category that is
present in a given county, at a given weekly interval throughout the year
(http://droughtmonitor.unl.edu/). Drought intensity categories correspond with well-known
objective indicators of drought (e.g., Palmer drought severity, standard precipitation, etc.) and
include D0, abnormally dry; D1, moderate drought; D2, severe drought; D3, extreme drought,
and D4, exceptional drought (Svoboda et al. 2002). These data ranged from 0-100 for any given
drought category and year, and were extracted and summed across drought categories for each

83

county, resulting in a cumulative drought severity and coverage index for each year. The range
of cumulative drought severity across all counties and years was 0-500, with a mean of 82.88.
The time period selected as the predictive index was the last week in August of each year, as HD
reporting for deer is seasonal and peaks during August and September (Howerth et al. 2001,
Stallknecht et al. 2015).
Because wetland cover is closely associated with Culicoides breeding habitat and reports
of HD related deer mortality (Purse et al. 2005, Gaydos et al. 2004, Berry et al. 2013), I created a
wetland classification in GIS from the 2006 National Land Cover Database (NLCD) across the
23 state area. I chose 2006 as an index of average available land cover data representing the
condition and least potential land cover change for a county within our 2000-2014 year dataset
(NLCD data are currently available for 2001, 2006, and 2011 only). I calculated the percentage
of wetland within each county based on 30 m pixels of wetland cover, with values ranging from
0% - 93% and a mean of 9%, and included this value as a covariate in the model (Fig. 3.1). I
predicted that the amount of wetlands interacts with drought conditions and then affects vector
capacity for HD transmission and subsequent HD occurrence. Given this association, I predicted
greater percentages of wetland cover would positively correlate with HD occurrence.
Large-scale HD reporting patterns have historically reflected physiographic variation in
addition to latitudinal variation (Stallknecht and Howerth 2004, Ruder et al. 2015a). To evaluate
the role of this environmental variable on HD presence, I used broad physiographic regions
across the 23 state study area including the Appalachian highlands, the Atlantic plain, the interior
highlands, the interior plain, and the Laurentian uplands (Fig. 3.2). These regions have broad
environmental characteristics that may interact with drought at regional scales and explain HD
transmission patterns. I assigned each county a categorical value, based on the physiographic
84

region that comprised the majority of the area within the county boundaries. No counties in the
Laurentian uplands were positive for HD in the study period, so were categorized as the most
similar physiographic region.
I predicted that a correlative temporal relationship existed in HD presence by county,
dependent on HD presence the previous year. Herd immunity would likely increase in the year
following an HD outbreak within a county due to an increase in acquired immunity for deer
surviving exposure to BTV or EHDV, potentially reducing HD reporting in subsequent years
(Couvillion et al. 1982, Gaydos et al. 2004). However, in a study by Xu and colleagues (2012), a
temporal autocorrelation term resulted in the odds of HD occurrence in a county being 0.52 times
higher if the county reported HD presence in the previous year. I calculated a temporal
autocorrelation term (T), similar to Xu and colleagues (2012) to evaluate the role of a one-year
time lag on the probability of HD presence in our study.

85

Figure 3.1 Wetland cover from 2006 NLCD data across the 1,830 county study area.

86

Figure 3.2 Physiographic regions of the Eastern US.

87

Data Analysis
To estimate the role of drought on the probability of HD mortality across the study area, I
used a generalized linear mixed model (GLMM) via the “lme4” package in R 3.4 (R
development team 2010). I developed a priori models based on predicted ecological relevance
for the spatial and temporal relationship between drought, environment, and HD occurrence,
where the response variable was always HD mortality (presence or absence). I fit models with
predictors including state and year as random effects, and drought severity, latitude, longitude,
physiographic region, and percentage wetland cover as fixed effects. Two random intercepts
were allowed to vary by state and year to account for variation among state disease reporting and
annual stochasticity respectively. Interaction terms were included between drought and latitude,
and drought and wetland (Table 3.1). A temporal autocorrelation term, 𝑇𝑖𝑗 , was included as a
fixed effect to account for HD presence in one year affecting the probability of HD the following
year for the same county, where 𝑇𝑖𝑗 , = 𝑌𝑖,𝑗+1 , with 𝑌𝑖,𝑗+1 being the HD response in county 𝑖 if
that county had HD present in the previous year. Following a binary GLMM for presenceabsence data, I used a logit-link format to estimate the probability of HD presence (𝑃), in a given
county (𝑖), and a given year (𝑗). The following general model structure was:

𝑙𝑜𝑔𝑖𝑡 (𝑃𝑖𝑗 ) = log

(𝑃𝑖𝑗 )
=
1 − 𝑃𝑖𝑗

𝛼 + 𝛽1 ∗ 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 + 𝛽2…5 ∗ 𝑅𝑒𝑔𝑖𝑜𝑛𝑖𝑗 + 𝛽6 𝑊𝑒𝑡𝑙𝑎𝑛𝑑𝑖 + 𝛽7 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒𝑖𝑗 + 𝛽8 𝐿𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑒𝑖𝑗
+ 𝛽9 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 ∗ 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒𝑖𝑗 + 𝛽10 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 ∗ 𝑊𝑒𝑡𝑙𝑎𝑛𝑑𝑖 + 𝛾𝑇𝑖𝑗 + 𝑌𝑒𝑎𝑟𝑗
+ 𝑆𝑡𝑎𝑡𝑒𝑖 ,
where 𝑌𝑒𝑎𝑟𝑗 ~𝑁(0, 𝜎𝑎2 ) and 𝑆𝑡𝑎𝑡𝑒𝑖 ~ 𝑁(0, 𝜎𝛿2 ).
88

Model component 𝛼 is the model intercept, 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 is the cumulative drought severity for the
last week of August in each year and county, 𝑅𝑒𝑔𝑖𝑜𝑛𝑖𝑗 are categorical variables for the 4
available physiographic regions for each county and year, 𝑊𝑒𝑡𝑙𝑎𝑛𝑑 is the percent wetland
coverage in each county, 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒 is the latitude of each county, 𝐿𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑒 is the longitude of
each county, 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 ∗ 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒𝑖 is the interaction between drought and latitude,
𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 ∗ 𝑊𝑒𝑡𝑙𝑎𝑛𝑑𝑖𝑗 is the interaction between drought and wetland, 𝛵𝑖𝑚𝑒𝑙𝑎𝑔 is the
temporal autocorrelation term, 𝑌𝑒𝑎𝑟𝑗 is the random intercept for year, and 𝑆𝑡𝑎𝑡𝑒𝑖 is a random
effect term for each state. The terms 𝛽1…𝑘 and 𝛾 are the fitted parameter coefficients for the
predictive variables. I used Akaike’s Information Criterion adjusted for small sample size
(AICc), to select the most parsimonious model which best explained the probability of HD
presence for all preliminary candidate models (Table 3.1, Burnham and Anderson 2002).
Univariate models based on the full candidate model were evaluated to assess collinearity and
relative variable importance (see Appendix F). All models that were evaluated were conditional
on the random effects structure selected initially using AIC. The area under the receiving
operating characteristic (ROC) curve was used to assess model fit (Hosmer and Lemeshow
2000).

89

Table 3.1 Notation, description, and number of parameters providing the base structure for a
priori models for estimating the probability of HD presence in a county and year during 20002014. All models are conditional on the random effects structure selected initially using AIC.
Model notation

Description

𝛼 + 𝛽1 ∗ 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 + 𝛽2…5 ∗ 𝑅𝑒𝑔𝑖𝑜𝑛𝑖𝑗
+ 𝛽6 𝑊𝑒𝑡𝑙𝑎𝑛𝑑𝑖 + 𝛽7 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒𝑖𝑗
+ 𝛽8 𝐿𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑒𝑖𝑗
+ 𝛽9 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 ∗ 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒𝑖𝑗

Full model; Drought, physiographic region, %
wetland, latitude, longitude, an interaction between
drought and latitude, an interaction between drought
and wetland, a 1 year time lag effect, a random
effect (intercept) for year, and a random effect
(intercept) for state

+ 𝛽10 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 ∗ 𝑊𝑒𝑡𝑙𝑎𝑛𝑑𝑖𝑗
+ 𝛾𝑇𝑖𝑚𝑒𝑙𝑎𝑔𝑖𝑗 + 𝑌𝑒𝑎𝑟𝑗 + 𝑆𝑡𝑎𝑡𝑒𝑖

𝛼 + 𝛽1 ∗ 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 + 𝛽6 𝑊𝑒𝑡𝑙𝑎𝑛𝑑𝑖
+ 𝛽7 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒𝑖𝑗
+ 𝛽8 𝐿𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑒𝑖𝑗
+ 𝛽9 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 ∗ 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒𝑖𝑗

Drought, % wetland, latitude, longitude, an
interaction between drought and latitude, an
interaction between drought and wetland, a 1 year
time lag effect, a random effect (intercept) for year,
and a random effect (intercept) for state

+ 𝛽10 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 ∗ 𝑊𝑒𝑡𝑙𝑎𝑛𝑑𝑖𝑗
+ 𝛾𝑇𝑖𝑚𝑒𝑙𝑎𝑔𝑖𝑗 + 𝑌𝑒𝑎𝑟𝑗 + 𝑆𝑡𝑎𝑡𝑒𝑖

𝛼 + 𝛽1 ∗ 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 + 𝛽2…5 ∗ 𝑅𝑒𝑔𝑖𝑜𝑛𝑖𝑗
+ 𝛽7 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒𝑖𝑗
+ 𝛽8 𝐿𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑒𝑖𝑗

Drought, physiographic region, latitude, longitude,
an interaction between drought and latitude, a 1 year
time lag effect, a random effect (intercept) for year,
and a random effect (intercept) for state

+ 𝛽9 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 ∗ 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒𝑖𝑗
+ 𝛾𝑇𝑖𝑚𝑒𝑙𝑎𝑔𝑖𝑗 + 𝑌𝑒𝑎𝑟𝑗 + 𝑆𝑡𝑎𝑡𝑒𝑖

𝛼 + 𝛽1 ∗ 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 + 𝛽7 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒𝑖𝑗
+ 𝛽8 𝐿𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑒𝑖𝑗
+ 𝛽9 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 ∗ 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒𝑖𝑗

Drought, latitude, longitude, an interaction between
drought and latitude, a 1 year time lag effect, an
random effect (intercept) for year, and a random
effect (intercept) for state

+ 𝛾𝑇𝑖𝑚𝑒𝑙𝑎𝑔𝑖𝑗 + 𝑌𝑒𝑎𝑟𝑗 + 𝑆𝑡𝑎𝑡𝑒𝑖

𝛼 + 𝛽1 ∗ 𝐷𝑟𝑜𝑢𝑔ℎ𝑡𝑖𝑗 + 𝛽7 𝐿𝑎𝑡𝑖𝑡𝑢𝑑𝑒𝑖𝑗
+ 𝛽8 𝐿𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑒𝑖𝑗 + 𝛾𝑇𝑖𝑚𝑒𝑙𝑎𝑔𝑖𝑗
+ 𝑌𝑒𝑎𝑟𝑗 + 𝑆𝑡𝑎𝑡𝑒𝑖

Drought, latitude, longitude, a 1 year time lag
effect, an random effect (intercept) for year, and a
random effect (intercept) for state

90

RESULTS
I found evidence that drought severity had a significant interaction with large-scale
patterns of herd immunity and wetland cover to explain changing patterns of HD presence across
1,830 counties from 2000-2014 (Fig. 3.3, Fig. 3.5). A total of 2,795 counties during the 15 year
time period reported presence of HD mortality. The best model using model selection criteria
included the full suite of random and fixed effects, except for longitude (See Appendix F). This
included 2 interaction terms (drought severity*latitude, drought severity*wetland), the associated
fixed effect terms (drought severity, latitude, wetland), categorical fixed effects for
physiographic regions and temporal autocorrelation, and random intercepts for state and year, all
which contributed to explaining the probability of disease presence in a given year and county
(Table 3.2). My statistical model shows clear evidence for the role of drought severity as a
predictor for HD presence when modeled as a function of either latitude or wetland cover.
The interaction between drought and latitude proved to be a positive predictor of HD
presence in counties that were greater than 33 degrees in latitude. In these counties, as drought
severity increased, the probability of HD increased (Fig. 3.3). The magnitude of predictive
strength increased significantly with increasing latitude, with a probability of HD presence
reaching nearly 0.40 in the most severe drought conditions at 45 degrees latitude (counties in the
states of Michigan and Wisconsin, Fig. 3.3, Fig. 3.4). Alternatively, drought severity was not
predictive in southern latitudes (latitudes of 30-35 degrees), with a probability of HD presence
consistently near 0.15 across drought conditions. High levels of population-level HD
seroprevalence have been demonstrated for this region and fewer susceptible animals are likely
(Howerth et al. 2001, Stallknecht et al. 1996, Park et al. 2013, Stallknecht et al. 2015). Further,
drought was an important predictor when interacted with wetland in the model. The total percent
91

of wetland cover for each county interacted with drought severity across the study area, reducing
the probability of HD with increasing wetland values. So, as wetland cover and drought severity
increased, the probability of HD presence decreased. However, the effect of this relationship was
minimal, becoming less significant with increasing percent wetland.
The physiographic region for each county was an important predictor in the best model as
a categorical variable, although the significance varied across the four regions. The Appalachian
Highlands was the intercept reference category and had a positive relationship with the
probability of HD presence. The odds of HD presence in the Interior Highlands was slightly
greater than the Appalachian highlands (9%), while the odds of HD presence in the Atlantic
plains or the Interior plains was lower than the Application highlands (37% and 12%,
respectively, Table 3.3).
Both random effects were included in all models based on AIC. For the best model, the
random effect average State variance was 0.82 (SD = 0.91) and the average Year variance was
0.81 (SD = 0.90). I included both latitude and longitude centroids for each county in my
predictive model to incorporate spatial covariates which would address my prediction of spatial
autocorrelation of HD presence among neighboring counties (Barry et al. 2013). However, the
longitude term was not a significant variable in any model and there was little ecological
justification for including this term. Thus, despite the global model having a competitive ΔAICc
value (within 2 units), I chose not to include this model in my final estimates as it was not more
informative and likely an over parameterization to improve model fit. Alternatively, latitude had
strong justification for inclusion as a spatially explicit predictive term and as an ecological
prediction. Latitude was a surrogate term for herd immunity in this model, given previous
evidence for endemic stability in southern states within the US (Stallkencht et al. 1996, Park et
92

al. 2013, Stallknecht et al. 2015). The interaction between drought severity and latitude, as
mentioned earlier, greatly varied as a risk factor for HD presence in a north-south gradient across
the eastern US. This tells us that the coefficient of drought severity depends on the value
of latitude and vice versa.
Finally, I found the temporal autocorrelation term was a strong positive predictor of HD
presence in a given county and year, similar to previous research using a similar temporal
covariate (Xu et al. 2012). If a given county had HD present in a given year, the following year
was twice as likely to have HD present in that county (Table 3.2). For the best model, I used the
area under the ROC curve to evaluate the sensitivity and specificity of the model and found the
best model had acceptable discrimination with a value of 0.726 (See Appendix D, Hosmer and
Lemeshow 2000).

93

Figure 3.3 The probability of HD mortality presence with associated standard error (shaded area)
in logit scale on the y-axis and cumulative drought severity on the x-axis. The relationships are
then grouped by latitude. As latitude increases the slope of the line increases.

94

Figure 3.4 The change in the marginal effect of latitude on the coefficient of drought severity
with associated standard error (shaded area). This demonstrates the statistical significance of
these conditional coefficients and the relationship between the interacting model terms (drought
severity and latitude). With increasing latitude the magnitude increases for the coefficient of
drought severity on HD mortality. A superimposed histogram is included at the bottom of the
marginal effect plot to view the distribution of the data.

95

Figure 3.5 The probability of HD mortality presence with associated standard error (shaded area) in logit scale given the
percent wetland cover by county and the interaction with drought severity

96

Table 3.2 Average coefficient estimates for the fixed effect parameters for HD mortality
presence by county across 23 states in the eastern US, for 2000-2014.
Fixed effects:
(Intercept)Region1
Drought Severity
Latitude
Wetland
Region2
Region3
Region4
TAuto
Drought: Latitude
Drought: Wetland

Estimate

Std.
Error

z value

Pr(>|z|)

Significance

Odds
Ratio

3.350
-0.016
-0.167
0.106
-0.455
0.085
-0.126
0.698
0.001

0.969
0.001
0.024
0.043
0.102
0.145
0.106
0.072
0.000

3.458
-21.478
-6.835
2.482
-4.464
0.589
-1.191
9.748
25.728

0.001
< 2e-16
0.000
0.013
0.000
0.556
0.234
< 2e-16
< 2e-16

***
***
***
*
***

***
***

28.503
0.984
0.846
1.112
0.635
1.089
0.881
2.010
1.001

-0.001

0.000

-3.765

0.000

***

0.999

DISCUSSION
I found that extreme drought, when interacting with latitude, is a strong predictor of the
occurrence of HD outbreaks in regions with emerging and incursive epizootic HD events. The
influence of a single climatic predictor, drought severity, has never been assessed at a large
spatial-temporal scale for HD outbreaks. Interestingly, this climatic force was most informative
when allowed to vary across a latitudinal gradient and with varying wetland cover at the county
level. Further, models were significantly improved when HD presence in a previous year was
accounted for, proving the importance of temporal effects between years. My research provides
spatially explicit evidence for the link between climate forces and emerging disease patterns
across the range of this disease using a widely available drought severity index.
I demonstrated the role of a specific climate variable, drought severity, as a predictor of
HD occurrence across the eastern US. My resulting spatio-temporal model confirmed the
importance of the interaction of drought severity with environmental characteristics to explain
changing patterns of HD presence. The interaction between drought and wetland cover, in

97

particular, was indicative of the vector ecology associated with HD. Culicoides species have
specific breeding habitat requirements which are typically associated with soils bordering
wetland systems and riparian habitats (Wilkening et al. 1985). Because these water-oriented
systems are inherently affected by changing water dynamics and drought conditions, the
interaction of drought and wetland in the model was important for understanding one aspect of
vector-driven spatial heterogeneity across our study. Increased wetland area in a county was
initially believed to equate to an increased potential for midge population abundance in late
summer. While I did not have information to inform this mechanistic relationship, I did find that
the predictive effect of increased percent wetland cover at a county scale varied with increased
drought severity (Fig. 3.5). At low levels of wetland cover, the relationship with drought was
slightly positive, increasing the probability of HD with increased drought severity (Fig. 3.5).
However, when wetland cover was within the top 20%, the relationship between wetlands and
drought shifted to that of a protective effect, where increased wetland cover reduced the
probability of HD presence (Fig. 3.5). This relationship may be explained by counties that are
primarily wetland cover having less available high-quality habitat for deer, and thus reducing the
overall likelihood of the vector interacting with the host. Alternatively, extremely high wetland
cover in an area may not be ideal habitat for midge species responsible for HD transmission.
The forces driving recent changes in HD patterns have not been well studied, particularly for
emergent regions of disease where mortality of white-tailed deer is most severe. This is due, in
part, to the complex epidemiology of HD. The transmission of this disease varies spatially in a
north-south gradient, mirroring a pattern of endemic stability across the known range of the
disease (Stallknecht et al. 1996, Park et al. 2013, Stallknecht et al. 2015). Endemic stability has
been demonstrated for HD, and results in heterogeneous disease transmission and herd immunity

98

across large spatial scales (Park et al. 2013, Howerth et al. 2001). Accordingly, in the most
southern counties in the study where HD is enzootic, drought severity was not an important
predictor of disease, lending to the theory of herd immunity as a driving epidemiological force in
endemic regions (Park et al. 2013). With precipitous increases in latitude, however, I
demonstrated the increasing strength of drought severity as a predictor of disease presence (Fig.
3.3, Fig. 3.4). Thus, I addressed my initial objective to determine the importance of variation in
drought severity as a risk factor over the present range of HD in the eastern United States. My
finding was similar to a recent study of West Nile virus (WNV) patterns in humans, where the
primary forces explaining interannual variation for this disease included drought and immunity,
where regions with low immunity were predicted to have the greatest WNV epidemics (Paull et
al. 2017). Vector-borne diseases are susceptible to climate variation and, given these findings,
one can assume that if extreme drought severity increases, increase in HD occurrence in
epizootic and incursive regions may increase.
Temporal autocorrelation was predictive of HD presence in white-tailed deer for a given
county and year, but with a different relationship to previous-year HD outbreaks than expected. I
predicted that if HD was present in a county in one year, the following year would have a
negative correlation in HD presence due to fewer susceptible deer being available for infection,
and thus mortality reporting. Rather, I found a positive relationship with previous-year HD
status, with the odds of disease reporting to be twice as likely if HD was reported the year before
in the same county (Table 3.2). While this relationship was initially counterintuitive, the
transmission and accumulation of herd immunity for white-tailed deer is locally patchy and may
not be homogenous within a county (Ruder et al. 2015a). Further, the relationship with increased
HD seroprevalence for a population is non-linear with respect to HD reporting (Park et al. 2013),

99

which challenges the assumption that increased HD reporting always results in increased herd
immunity. Regardless, the correlation between adjacent reporting years was an important
predictor of HD presence over the 15-year study period.
Vector-borne diseases are considered to be most at risk of increasing their distribution
and frequency with predicted climate change scenarios (Harvell et al. 2002). The population
ecology and distribution of Culicoides species plays a significant role in the epidemiology of
hemorrhagic diseases (Mellor et al. 2000, Wilson and Mellor 2008, Savini et al. 2011). Midges,
in particular, have recently expanded their range northward in many regions, which has been
largely attributed to milder winter temperatures (Harvell et al. 2002). Due to the vector-borne
nature of HD, the prevalence of the disease is frequency-dependent and influenced by climatic
variables (Gaydos et al. 2004, Sleeman et al. 2009, Xu et al. 2012, Ruder et al. 2015). In North
America, adult-stage midges associated with HD die shortly after the first frost in autumn
(Nettles and Stallknecht 1992). If predictions of warmer average temperature are correct, the
basic reproduction number (R0) of pathogens is expected to increase and the pathogen will
persist longer on the landscape (Harvell et al. 2002). This spatial-temporal drought severity
model will provide insight for predicting areas at risk for future disease outbreaks.
MANAGEMENT IMPLICATIONS
My research was intended to better understand the drivers of HD across a large spatialtemporal gradient experiencing high variation in HD mortality events. Outbreaks of HD often
result in a devastating impact on local deer populations in epizootic and incursive disease regions
(Fischer et al. 1995, Gaydos et al. 2004, Stevens et al. 2015). However, managing these local
mortality events is extremely difficult due to multiple factors affecting the high annual variation
of the disease. Host dynamics can be complex and affect disease transmission dynamics. Based
100

on this study, the way in which drought severity and herd immunity interact may have important
consequences for local and regional management of deer populations.
The population impacts of increasing HD mortality and the subsequent management
implications are complicated by herd immunity dynamics, which is demonstrated by variation of
HD presence across latitude in the study. Given the apparent high mortality rates for naĂŻve deer
populations, HD can be considered to have top-down effects on deer abundance, producing
direct numeric reductions in the population (Gaydos et al. 2004, Holdo et al. 2009). Recent work
by Park et al. (2013) demonstrated that across a 16 state region, a greater probability of reporting
deer with HD symptoms occurred when deer population seroprevalence was intermediate
(around 50%). However, this study did not include recent disease reports in the upper Midwest
or Northeastern US where severe clinical disease and mortality have been observed (Stallknecht
et al. 2015). Typically, a disease outbreak that occurs within the lifespan of an animal will result
in antibody development or death upon infection, leaving a percentage of the population with
immunity to the pathogen. If a local deer population is not able to accumulate sufficient
immunity from pathogen exposure among individuals, deer populations may be prone to
continued die-offs given optimal environmental conditions (including severe drought) for disease
transmission. The effects of disease on the population dynamics of wildlife species depend on
the epidemiological zone the population is within and the associated epizootiology for the given
species (Gulland 1995, Osnas et al. 2009, Ruder et al. 2015). Epizootic disease events may result
in mass mortalities of animals which, coupled with stochastic events, can lead to extirpation of
local populations or species extinction (Harwood and Hall 1990). White-tailed deer are a highly
valued native wildlife species across the US and may be at risk of increased disease events in
regions of emergent HD.

101

Unfortunately, HD will remain difficult to manage in free-ranging deer populations. HD
mortality is not likely to increase with greater deer density and is not assumed to be sex- or ageclass specific (Fischer et al. 1995, Begon et al. 2002). However, maternal antibodies are
transmitted to fawns and persist for up to 18 weeks, resulting in protection of young deer from
late summer EHD outbreaks, regardless of drought conditions (Gaydos et al. 2002c). Over time,
deer populations will recolonize as reproduction continues and new deer immigrate into an area,
but this may be influenced by social structuring and movement rates (Porter et al. 1991). Midge
control options may be viable for isolated, localized wetlands, but is not a practical management
option at larger scales. Further, the relationship between wetlands and HD transmission appears
to vary at differing spatial scales, indicating more research is needed to understand the ecological
relationship between the disease vector and wetland dynamics.
HD is an emerging and reemerging issue globally. Similar to increasing HD presence in
northern US states, several other countries have observed increases in morbidity and mortality of
domestic ungulates infected with bluetongue virus. Emergent bluetongue outbreaks in Europe,
Northern Africa, and the Middle East have resulted in large mortality events and economic
losses. Researchers have found that wind dispersal and warmer temperatures have resulted in the
successful invasion of new vector species and subsequent disease occurrence (Purse et al. 2005).
While North America has not suffered the same type of loses from bluetongue virus in domestic
animals, there is concern that newly emergent HD serotypes may result in increased mortality for
both wild and domestic animals in the future. Continued efforts to understand the role of climate
on vector-borne disease patterns will be of increasing importance as climate change occurs and
diseases continue to spread in to naĂŻve populations.

102

APPENDICES

103

APPENDIX A: MODEL VALIDATION RESULTS FOR CHAPTER 3

Area under the ROC curve for best model: 0.7252321
Associated R code for AUC value: auc(roc(prob2,as.factor(FULLFE2$model$HDMORT)))

104

APPENDIX B: MODEL OUTPUT FOR BEST MODEL FOR CHAPTER 3

HD32 output
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation)
['glmerMod']
Family: binomial ( logit )
Formula: HDMORT ~ SUM * CENTROID_Y.y + SUM * Wetland + SUM +
CENTROID_Y.y +
Wetland + physio + (1 | fYear) + (1 | fState.y) + TAuto
Data: Alldatacombo
AIC BIC logLik deviance df.resid
13128.7 13226.5 -6552.4 13104.7 25608
Scaled residuals:
Min 1Q Median 3Q Max
-2.3936 -0.2974 -0.1912 -0.1277 14.7954
Random effects:
Groups Name
Variance Std.Dev.
fState.y (Intercept) 0.8141 0.9023
fYear (Intercept) 0.8059 0.8977
Number of obs: 25620, groups: fState.y, 23; fYear, 14
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)
3.35E+00
9.69E-01
3.458 0.000543
***
SUM -1.60E-02
7.46E-04
-21.478
< 2e-16
***
CENTROID_Y.y
-1.67E-01
2.44E-02
-6.835 8.22E-12
***
Wetland
1.06E-01
4.28E-02
2.482 0.013083
*
physio2
-4.55E-01
1.02E-01
-4.464 8.04E-06
***
physio3
8.55E-02
1.45E-01
0.589 0.555867
physio4
-1.26E-01
1.06E-01
-1.191 0.233563
TAuto 6.98E-01
7.16E-02
9.748 < 2e-16
***
SUM:CENTROID_Y. y 5.220e-04 2.03E-05
25.728 < 2e-16
***
SUM:Wetland -8.67E-04
2.30E-04
-3.765 0.000166
***
--Signif. codes: 0 ‘***’ 0.0
01 ‘**’ 0.0
1 ‘*’ 0. 05 ‘.’ 0.
1‘’1
Correlation of Fixed Effects:
(Intr) SUM CENTRO Wetlnd physi2 physi3 physi4 TAuto SUM:CE
SUM
-0.064
CENTROID_Y. -0.946 0.066
Wetland -0.097 0.076 0.120
105

physio2
physio3
physio4

-0.277 -0.007 0.245 -0.403
-0.175 0.011 0.144 -0.113 0.394
0.092 0.011 -0.148 -0.170 0.316 0.540

TAuto

-0.050 -0.004 0.039 -0.007 0.037 0.010 0.001

SUM:CENTROID 0.078 -0.967 -0.085 -0.080 -0.007 -0.009 -0.012 0.010
SUM:Wetland 0.147 -0.028 -0.152 -0.321 -0.199 -0.075 -0.014 0.001 0.097

106

APPENDIX C: SUPPLEMENTARY MODEL SELECTION DATA FOR CHAPTER 2
Table C.1 Full model selection table using AIC criteria for all model variations predicting the
probability of HD presence in 1830 counties within 23 states during 2000-2014.

Model Name

K

AICc

Delta_AICc

AICcWt

Cum.Wt

LL

Drought*Latitude, Drought*Wetland +

12

13129.7

0

0.71

0.71

-6552.84

Global Model

13

13131.55

1.85

0.28

0.99

-6552.77

Drought*Latitude, + physio +T

10

13139.51

9.81

0.01

0.99

-6559.75

Drought*Latitude, + wetland + Physio +T

11

13140.24

10.54

0

1

-6559.11

Drought*Latitude + Wetland + Physio +

12

13141.93

12.23

0

1

-6558.96

Drought*Latitude, Drought *Wetland +T

9

13149.38

19.68

0

1

-6565.68

Drought*Latitude, + Wetland + T

8

13168.23

38.53

0

1

-6576.11

Drought*Latitude, + Long +T

8

13174.34

44.64

0

1

-6579.17

Drought *Wetland + Latitude + Long +

12

13216.18

86.48

0

1

-6596.09

Drought *Wetland +T

8

13233.75

104.05

0

1

-6608.87

Drought + Latitude + physio + Wetland

10

13254.65

124.95

0

1

-6617.32

11

13256.43

126.73

0

1

-6617.21

Drought + Physio + Latitude +T

9

13260.06

130.36

0

1

-6621.03

Drought + Latitude + Long + physio +T

10

13260.92

131.22

0

1

-6620.46

Drought + Latitude + T

6

13291.16

161.46

0

1

-6639.58

Drought + Latitude + Wetland + T

7

13292.08

162.38

0

1

-6639.04

Drought + Physio

8

13293.72

164.02

0

1

-6638.86

Drought + Latitude + Long + Wetland

8

13293.98

164.28

0

1

-6638.99

Drought + Long + T

6

13316.88

187.18

0

1

-6652.44

Drought + Wetland

6

13316.99

187.29

0

1

-6652.5

wetland + physio

8

13662.93

533.23

0

1

-6823.46

Latitude + T

5

13664.22

534.52

0

1

-6827.11

Temporal autocorrelation (T)

4

13686.02

556.32

0

1

-6839.01

Drought Latitude interaction

6

13839.47

709.77

0

1

-6913.73

Drought Physio Inter

10

13848.59

718.89

0

1

-6914.29

Physio + T

long +T

Physio +T

+T
Drought + Latitude + Long + Wetland +
Physio +T

107

Table C.1 (cont’d)
Drought Wetland Inter + Latitude

7

13896.97

767.27

0

1

-6941.48

Drought Wetland Interaction

6

13918.24

788.54

0

1

-6953.12

Latitude + Drought Categories

9

13959.1

829.4

0

1

-6970.55

Longitude Drought interaction

6

13977.4

847.7

0

1

-6982.7

Drought categories

8

13980.75

851.05

0

1

-6982.37

Drought

4

13998.9

869.2

0

1

-6995.45

Physio

6

14345.75

1216.05

0

1

-7166.87

Latitude + Wetland + T

5

14354.81

1225.11

0

1

-7172.4

Latitude

4

14368.94

1239.24

0

1

-7180.47

Wetland

4

14382.19

1252.49

0

1

-7187.09

Longitude

4

14386.59

1256.89

0

1

-7189.3

Intercepts only null model

3

14388.51

1258.81

0

1

-7191.26

108

LITERATURE CITED

109

LITERATURE CITED

Altizer, S., R. S. Ostfeld, P. T. Johnson, S. Kutz, and C. D. Harvell. 2013. Climate change and
infectious diseases: from evidence to a predictive framework. Science 341: 514-519.
Berry, B. S., K. Magori, A. C. Perofsky, D. E. Stallknecht, and A. W. Park. 2013. Wetland cover
dynamics drive hemorrhagic disease patterns in white-tailed deer in the United States.
Journal of Wildlife Diseases 49: 501-509.
Begon M., M. Bennett, R. G. Bowers, N. P. French, S. M. Hazel, and J. Turner. 2002. A
clarification of transmission terms in host-microparasite models: numbers, densities and
areas. Epidemiology and Infection 129: 147-153.
Burnham, K. P., and D. R. Anderson. 2002. Model selection and inference: a practical
information-theoretic approach. Springer Publishing, Second edition, New York, New
York, USA.
Coleman, P. G., B. D. Perry, and M. E. J. Woolhouse. 2001. Endemic stability—a veterinary
idea applied to human public health. The Lancet 357: 1284-1286.
Couvillion, C. E., V. F. Nettles, and W. R. Davidson. 1981. Hemorrhagic disease among whitetailed deer in the southeast from 1971through 1980. In Proceedings of the 85th annual
meeting of the United States Animal Health Association, St. Louis, Missouri 522–537.
ESRI [Environmental Systems Research Institute.]. 2013. ARC/INFO. Version 10.0
Environmental Systems Research Institute, Redlands, California, USA.
Fischer, J. R., L. P. Hansen, J. R. Turk, M. A. Miller, W. H. Fales, and H. S. Gosser. 1995. An
epizootic of hemorrhagic disease in white-tailed deer (Odocoileus virginianus) in
Missouri: necropsy findings and population impact. Journal of Wildlife Diseases 31: 3036.
Gaydos, J. K., W. R. Davidson, E.W. Howerth, M. Murphy, F. Elvinger, D. E. Stallknecht.
2002a. Cross-protection between epizootic hemorrhagic disease virus serotypes 1 and 2
in white-tailed deer. Journal of Wildlife Diseases 38: 720-728.
Gaydos, J. K., W. R. Davidson, F. Elvinger, D. G. Mead, E. W. Howerth, and D. E. Stallknecht,
2002b. Innate resistance to epizootic hemorrhagic disease in white-tailed deer. Journal of
Wildlife Diseases 38 4:713-719.
Gaydos, J. K., D. E. Stallknecht, D. Kavanaugh, R. J. Olson, and E. R. Fuchs. 2002c. Dynamics
of maternal antibodies to hemorrhagic disease viruses (Reoviridae: Orbivirus) in whitetailed deer. Journal of Wildlife Diseases 38: 253-257.

110

Gaydos, J. K., J. M. Crum, W. R. Davidson, S. S. Cross, S. F. Owen, and D. E. Stallknecht.
2004. Epizootiology of an epizootic hemorrhagic disease outbreak in West Virginia.
Journal of Wildlife Diseases 40: 383-393.
Gibbs E. P. J., and E. C. Greiner. 1989. Bluetongue and epizootic hemorrhagic disease. In:
Monath T, ed. Arboviruses: Epidemiology and Ecology Boca Raton, FL: CRC Press Inc.
2: 39-70.
Guis, H., C. Caminade, C. Calvete, A. P. Morse, A. Tran, and M. Baylis. 2012. Modelling the
effects of past and future climate on the risk of bluetongue emergence in Europe. Journal
of the Royal Society Interface 9: 339-350.
Gulland, F. M. D. 1995. The impact of infectious diseases on wild animal populations: a
review. Ecology of infectious diseases in natural populations. Cambridge University
Press, Cambridge 20-51.
Hartemink, N. A., B. V. Purse, R. Meiswinkel, H. E. Brown, A. De Koeijer, A. R. W. Elbers, and
J. A. P. Heesterbeek. 2009. Mapping the basic reproduction number (Ro) for vectorborne diseases: A case study on bluetongue virus. Epidemics 1: 53-161.
Harvell, C. D., C. E. Mitchell, J. R. Ward, S. Altizer, A. P. Dobson, R. S. Ostfeld, and M. D.
Samuel. 2002. Climate warming and disease risks for terrestrial and marine
biota. Science 296: 2158-2162.
Harwood, J., and A. Hall. 1990. Mass mortality in marine mammals: its implications for
population dynamics and genetics. Trends in Ecology and Evolution 8: 254-257.
Hewitt, D. G. ed. 2011. Biology and management of white-tailed deer. CRC Press, 2011.
Holdo, R. M., A. R. Sinclair, A. P. Dobson, K. L. Metzger, B. M. Bolker, M. E. Ritchie, and R.
D. Holt. 2009. A disease-mediated trophic cascade in the Serengeti and its implications
for ecosystem. PLoS biology, 7:9
Hosmer, D. W., and S. Lemeshow. 2000. Applied logistic regression. John Wiley and Sons, New
York, New York, USA.
Howerth, E. W., D. E. Stallknecht, and P. D. Kirkland. 2001. Bluetongue, epizootic hemorrhagic
disease, and other orbivirus-related diseases. In: Infectious Diseases of Wild Mammals.
3rd edition. Edited by E.S. Williams and I.K. Barker. Ames, Iowa: Iowa State Press;
2001:77–97.
Jones, K. E., N. G. Patel, M. A. Levy, A. Storeygard, D. Balk, J. L. Gittleman, and P. Daszak.
2008. Global trends in emerging infectious diseases. Nature 451: 990-993.
Koelle, K., X. RodĂł, M. Pascual, M. Yunus, and G. Mostafa. 2005. Refractory periods and
climate forcing in cholera dynamics. Nature 436: 696-700.
111

Mellor, P.S., and J. Boorman. 1995. The transmission and geographical spread of African horse
sickness and bluetongue viruses. Annals of Tropical Medicine & Parasitology 89: 1-15.
Mellor, P.S., J. Boorman, and M. Baylis. 2000. Culicoides biting midges: their role as arbovirus
vectors. Annual Review of Entomology 45: 307-340.
Mullens, B. A., A. C. Gerry, T. J. Lysyk, and E. T. Schmidtmann. 2004. Environmental effects
on vector competence and virogenesis of bluetongue virus in Culicoides: interpreting
laboratory data in the field. Veterinaria Italiana 40: 160-163.
Nettles V. F., W. R. Davidson, and D. E. Stallknecht. 1992. Surveillance for hemorrhagic disease
in white-tailed deer and other wild ruminants. Proceedings of the Southeast Association
of Fish Wildlife Agencies 46: 138-146.
Osnas, E. E., D. M. Heisey, R. E. Rolley, and M. D. Samuel. 2009. Spatial and temporal patterns
of chronic wasting disease: fine‐scale mapping of a wildlife epidemic in
Wisconsin. Ecological Applications 19: 1311-1322.
Park, A. W., K. Magori, B. A. White, and D. E. Stallknecht. 2013. When more transmission
equals less disease: reconciling the disconnect between disease hotspots and parasite
transmission. PloS one 8: 4.
Paull, S. H., D. E. Horton, M. Ashfaq, D. Rastogi, L. D. Kramer, N. S. Diffenbaugh, and A. M.
Kilpatrick. 2017. Drought and immunity determine the intensity of West Nile virus
epidemics and climate change impacts. Proceedings of the Royal Society B 284: 20162078.
Purse, B. V., P. S. Mellor, D. J. Rogers, A. R. Samuel, P. P. Mertens, and M. Baylis. 2005.
Climate change and the recent emergence of bluetongue in Europe. Nature Reviews
Microbiology 3: 171-181.
Porter, W. F., N. E. Mathews, H. B. Underwood, R. W. Sage, and D. F. Behrend. 1991. Social
organization in deer: implications for localized management. Environmental
Management 15: 809-814.
Rogers, D. J. and S. E. Randolph. 2006. Climate change and vector-borne diseases. Advances in
Parasitology 62: 345-381.
Ruder, M. G., T. J. Lysyk, D. E. Stallknecht, L. D. Foil, D. J. Johnson, C.C. Chase, D. A.
Dargatz, and E. P. J. Gibbs. 2015a. Transmission and epidemiology of bluetongue and
epizootic hemorrhagic disease in North America: Current perspectives, research gaps,
and future directions. Vector-Borne and Zoonotic Diseases 15: 348-363.
Ruder, M. G., D. E. Stallknecht, E. W. Howerth, D. L. Carter, R. S. Pfannenstiel, A. B. Allison,
and D. G. Mead. 2015b. Effect of temperature on replication of epizootic hemorrhagic
112

disease viruses in Culicoides sonorensis (Diptera: Ceratopogonidae). Journal of Medical
Entomology 52: 1050-1059.
R Development Core Team. 2011. R: A language and environment for statistical computing. R
Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL
http://www.R-project.org.
Savini, G., A. Afonso, P. Mellor, I. A. O. Aradaib, H. Yadin, M. Sanaa, W. Wilson, F. Monaco,
and M. Domingo. 2011. Epizootic heamorragic disease. Research in Veterinary
Science 91:1-17.
Shope, R. E., L. G. MacNamara, and R. Mangold. 1960. A virus-induced epizootic hemorrhagic
disease of the Virginia white-tailed deer (Odocoileus virginianus). The Journal of
Experimental Medicine 111:155-170.
Sleeman J. M., J. E. Howell, W. M. Knox, P. J. Stenger. 2009. Incidence of hemorrhagic disease
in white-tailed deer is associated with winter and summer climatic conditions. Eco Health
6: 11-15.
Stallknecht, D. E., M. L. Kellogg, J. L. Blue, and J. E. Pearson, 1991. Antibodies to bluetongue
and epizootic hemorrhagic disease viruses in a barrier island white-tailed deer
population. Journal of Wildlife Diseases 27: 668-674.
Stallknecht, D. E., M. P. Luttrell, K. E. Smith, and V. F. Nettles. 1996. Hemorrhagic disease in
white-tailed deer in Texas: a case for enzootic stability. Journal of Wildlife Diseases 32:
695-700.
Stallknecht, D. E. and E. W. Howerth. 2004. Epidemiology of bluetongue and epizootic
hemorrhagic disease in wildlife: surveillance methods. Veterinaria Italiana 40: 203-207.
Stallknecht, D. E., A. B. Allison, A.W. Park, J. E. Phillips, V. H. Goekjian, V. F. Nettles, and J.
R. Fischer. 2015. Apparent increase of reported hemorrhagic disease in the midwestern
and northeastern USA. Journal of Wildlife Diseases 51: 348-361.
Svoboda, M., D. LeComte, M. Hayes, R. Heim, K. Gleason, J. Angel, B.
Rippey, R. Tinker, M. Palecki, D. Stooksbury, D. Miskus, and S. Stephens. 2002. The
drought monitor. Bulletin of the American Meteorological Society 83: 1181-1192
Stevens, G., B. McCluskey, A. King, E. O’Hearn, and G. Mayr. 2015. Review of the 2012
epizootic hemorrhagic disease outbreak in domestic ruminants in the United States. PloS
one 10: e0133359
Thomas, F. C., N. Willis, and G. Ruckerbauer. 1974. Identification of viruses involved in the
1971 outbreak of hemorrhagic disease in southeastern United States white-tailed deer.
Journal of Wildlife Diseases 10:187-189.

113

Wilkening, A.J., D. L. Kline, and W. W. Wirth. 1985. An annotated checklist of the
Ceratopogonidae (Diptera) of Florida with a new synonymy. Florida Entomologist 511537.
Wittmann, E. J., and M. Baylis. 2000. Climate change: effects on Culicoides-transmitted viruses
and implications for the United Kingdom. The Veterinary Journal 160: 107-117.
Wilson, A., and P. Mellor. 2008. Bluetongue in Europe: vectors, epidemiology and climate
change. Parasitology Research 103: S69-77.
Xu, B., M. Madden, D. E. Stallknecht, T. W. Hodler, and K. C. Parker. 2012. Spatial-temporal
model of hemorrhagic disease in white-tailed deer in south-east USA, 1983 to 2000.
Veterinary Record 170: 288.

114

POSTFACE

By investigating the local distribution and population-level impacts on deer in areas
affected by EHD, developing and applying population estimation techniques appropriate for
localized population change assessment, and evaluating landscape-level drivers of hemorrhagic
disease across an epidemiological gradient, I was able to address uncertainties on the population
impacts and risks of emergent EHD for white-tailed deer. I briefly review my findings in this
section and discuss relevant research needs for future study.
Chapter 1 sought to determine whether localized deer population impacts occurred after a
severe EHD outbreak in 2012 in Michigan. The experimental study design implemented allowed
me to relate population impacts to EHD and demonstrate the local spatial distribution of
population impacts relative to a high-risk riparian habitat. Analyses confirmed that population
impacts did occur and were attributable to EHD mortality. These impacts were most severe
within 1 km of a river corridor with known EHD reports and dissipated within 5 km of an EHD
focal area. These findings were substantiated by comparison to a non-disease site within
Michigan, where deer populations were stable across habitat types and over time. Further, where
population impacts were observed, we consistently detected annual population increases over the
4-year monitoring period following the disease outbreak. This increase over time provided
evidence of deer population recovery after a sudden and localized population depression.
Chapter 2 sought to develop and apply population estimation methods for localized
assessment of deer populations. This was accomplished using aerial surveys and implementing
both distance sampling (Buckland et al. 2001) and N-mixture models (Royle 2004) for count
data collected from line-transects. Further, N-mixture models were evaluated for sensitivity to
115

spatial unit size selection. Analyses of this chapter demonstrated similar population estimates
beween analytical frameworks, indicating that N-mixture models were successfully applied to
aerial survey data for deer. Increased precision of abundance estimates for the N-mixture models
compared to distance sampling estimates resulted in the clear demonstration of a population
increase over time (abundance was greater in 2016 than 2014). These results mirror findings
from ground surveys conducted annually for the same study area during summer months.
Population impacts from EHD mortality in 2012 were suspected in this area and population
increases align with these predictions. Finally, we found that abundance and detection
probability estimates for N-mixture models was highly sensitive to the selection of spatial unit
size. Spatial units most similar to the movement ecology of deer (2000 m x 900 m) provided
realistic values and acceptable variance for abundance and detection parameters.
Chapter 3 sought to evaluate the role of large-scale drought severity across the
epidemiological range of hemorrhagic disease (HD) in the eastern United States. A 15-year data
set of HD occurrence by county from an annual survey conducted by the Southeastern
cooperative wildlife disease study was analyzed in this chapter. I developed a spatial temporal
model to evaluate drought severity and related environmental variables that might explain
changing patterns of hemorrhagic disease presence. My results indicated that drought severity
had a significant interaction with herd immunity, modeled as on index of latitude. As latitude
increased into northern regions were HD outbreaks are epidemic or incursive, the predictive
value for drought severity on HD presence increased significantly. This relationship was not
important in the southeastern US, where high seroprevelence and frequent HD events occur.
Further, the presence of HD was positively related to whether disease had occurred in the county
in the previous year and if the county existed in the Interior Highlands physiographic region.
116

These findings provide clear evidence of climate-based drivers of HD across the range of the
disease.
Areas of future research to inform EHD epidemiology and wildlife impacts should
include proactive herd immunity surveillance, regional vector confirmation, and refined largescale modeling efforts. These areas of focus will fill gaps in our understanding of the
epidemiology of the disease and improve our ability to predict future mortality events.
Systematic seroprevalence information for Michigan and other states experiencing epidemic
EHD outbreaks should be collected to create a clear spatial distribution of enzootic, epidemic,
and incursive regions of the disease. Seroprevalence within a deer population will inform the
proportion of susceptible individuals within the population and subsequently the risk associated
with future HD outbreaks. Culicoides surveys using standard protocols should be conducted in
conjunction with outbreaks to help determine regionally specific vectors and life-stage
requirements. Continued landscape-level modeling efforts of characteristics affecting midge
presence and disease occurrence, and forecasting models predicting EHD outbreaks under
varying immunity, climate, and land use scenarios would be informative for understanding
current and future conditions conducive to HD presence. Finally, state agencies should develop
systematic and consistent surveillance plans for assessing outbreak impacts of EHD. If deer
mortality information were collected using similar methods and survey designs across regions,
outbreak impacts could be mapped and compared for disease management and epidemiological
characterization.

117

LITERATURE CITED

118

LITERATURE CITED

Buckland, S. T., D. R. Anderson, K. P. Burnham, J. L. Laake, D. L. Borchers, and L. Thomas.
2001. Introduction to distance sampling: estimating abundance of biological populations.
Oxford University Press, Oxford, United Kingdom.
Royle, J. A. 2004. N‐mixture models for estimating population size from spatially replicated
counts. Biometrics 60: 108-115.

119