Patterns of Hemorrhagic Disease in White-Tailed Deer (Odocoileus virginianus) in the Great Plains of the USA, 1982–2020

Hemorrhagic disease (HD) of white-tailed deer (WTD, Odocoileus virginianus) is caused by epizootic hemorrhagic disease virus (EHDV) and bluetongue virus (BTV), which may also cause disease in some other ruminant species (Nettles et al. 1992). Both EHDV and BTV are in the genus Orbivirus (Reoviridae family) and are vectored by Culicoides biting midges (Stallknecht et al. 2002). Currently, Culicoides sonorensis and Culicoides insignis are the only confirmed vector species for EHDV and BTV, but other Culicoides species are suspected vectors (McGregor et al. 2022). Hemorrhagic disease is considered one of the most important viral diseases of WTD in North America (Stallknecht and Howerth 2004). Hemorrhagic disease in deer ranges from inapparent infection to acute and fatal febrile disease typically characterized by edema, effusion, and hemorrhage in multiple organs and tissues (Howerth et al. 2001). Some deer may survive acute HD infection to develop chronic sequelae, such as cracked hoof walls or scarred ruminal papillae (Couvillion et al. 1981).

Reports of cyclical WTD mortality events consistent with HD have been documented in the US since the late 1800s (Nettles and Stallknecht 1992), but it was not until the 1950s and 1960s that the viral causes of HD were confirmed (Shope et al. 1960; Stair et al. 1968). Following this, classical and molecular viral diagnostics enabled long-term monitoring of HD in the US and helped to characterize the diversity of viruses responsible for outbreaks (Stallknecht and Howerth 2004; Ruder et al. 2016). Annual HD surveillance indicates that EHDV-2 is the predominant virus causing HD outbreaks among WTD in the US, but other EHDV and BTV serotypes are commonly involved with HD outbreaks in wild ruminants (Ruder et al. 2016; Rivera et al. 2021). Although severe disease associated with EHDV and BTV infection is most commonly reported in WTD, other wild ruminant species commonly affected during outbreaks include mule deer (Odocoileus hemionus), bighorn sheep (Ovis canadensis), and pronghorn (Antilocapra americana; Howerth et al. 2001).

Clinical disease associated with HD in individual deer is variable, but general patterns of infection and disease in WTD populations vary predictably depending on the geographic region (Stallknecht et al. 2002; Ruder et al. 2015). Historically, EHDV and BTV circulated broadly between the latitudes of 35° S and 40–45° N in areas with habitats and climatic conditions that support Culicoides spp. vector populations (Purse et al. 2005). In the US, geographic areas may be described as having either epizootic or enzootic disease patterns (Stallknecht et al. 2002). Regions with enzootic patterns are generally at more southern latitudes where EHDV and BTV frequently circulate but reports of HD in WTD are less common (e.g., 2–3-yr cycles), because inapparent infections or mild disease without mortality predominate (Couvillion et al. 1981; Nettles et al. 1992). Regions with epizootic patterns generally occur at more northern latitudes where virus circulation is intermittent, herd immunity is low, and intense HD epizootics with high mortality occur infrequently (e.g., 8–10-yr cycle historically) (Stallknecht et al. 1996; Ruder et al. 2015; Stilwell et al. 2021). These epidemiologic patterns have been best characterized for states in the eastern and southeastern US, ranging from enzootic patterns throughout much of the coastal plain physiographic region in the southeast, to epizootic patterns in the upper midwest, northeast, and Piedmont and Appalachian Mountain areas in the southeast (Stallknecht et al. 1991; Nettles and Stallknecht 1992; Stallknecht et al. 2002). Over the last two decades northward expansion of HD has occurred in parts of the northeastern and midwestern US, highlighting the dynamic nature of this disease system (Stallknecht et al. 2015; Casey et al. 2021). In addition to these changing regional patterns, northward shifts in latitude of reported HD have occurred within individual states, such as Illinois (Boyer et al. 2007; Dorak et al. 2022). Much of the observed regional variation in infection outcome is related to differences in population immunity, including previous exposure to the viruses, passive antibody transfer to fawns, as well as genetic factors influencing innate immunity (Stallknecht et al. 1996; Gaydos et al. 2002a, b).

Sustained and long-term monitoring of HD in wild ruminant populations by wildlife management agencies is important for understanding epidemiological patterns in the US and how cyclical variation in outbreaks might affect population dynamics and management goals. Given the cyclical nature of HD in areas where deer are immunologically naïve to EHDV or BTV infection (e.g., northern latitudes), it is difficult to predict when outbreaks might occur. However, passive surveillance and outbreak investigations have helped to identify new virus serotypes (Ruder et al. 2017), novel patterns of disease (Stallknecht et al. 2015; Dorak et al. 2022), and potential risk factors such as severe drought (Christensen et al. 2020). Although investigation of an EHDV-2 outbreak in West Virginia revealed an estimated 20% mortality rate in the affected WTD population (Gaydos et al. 2004), questions remain regarding the impact of more frequent HD outbreaks on a population. Sustained monitoring of HD outbreaks is fundamentally important for understanding the potential population impacts of HD, especially in areas where deer populations are immunologically naïve and there is evidence of northern expansion and increased frequency and intensity of HD outbreaks (Stallknecht et al. 2015; Dorak et al. 2022).

Historically, severe HD outbreaks have been reported in parts of the Great Plains of the US, including the states of North Dakota, South Dakota, Nebraska, and Kansas (Stevens et al. 2015). However, HD in the Great Plains remains poorly characterized. Describing historic and current patterns of reported HD is the first step in recognizing changes in this disease system and how HD patterns might change in the future. Anecdotally, patterns of HD in this region may be changing, similar to the documented northern expansion and increased frequency and intensity of HD in parts of the eastern US (Stallknecht et al. 2015). The objectives of our study were 1) to broadly describe historic patterns of HD in the Great Plains through the use of annual survey data between 1982 and 2020; and 2) to determine if patterns of reported HD changed in frequency, intensity, or geographic distribution during this time period.

The Great Plains encompasses approximately 2.17 million square kilometers across the US states of Montana, Wyoming, Colorado, New Mexico, Texas, Oklahoma, Missouri, Iowa, Minnesota, North Dakota, South Dakota, Nebraska, Kansas, and Oklahoma (Zhu et al. 2011). The Great Plains encompasses highly diverse habitats and spans the latitudinal gradient that has been significant in determining risk of HD outbreaks in other states (Omernik 1987; Christensen et al. 2020; Casey et al. 2021). Our study focused on the interior portion of the region that includes the five states of North Dakota, South Dakota, Nebraska, Kansas, and Oklahoma (hereafter referred to as the Great Plains), because HD patterns have yet to be described in this region (Supplementary Material Fig. S1). Our study area ranged from 33–49° N latitude and 104–96° W longitude, including approximately 958,300 km².

The reports of HD between 1982 and 2020 were obtained from a long-term (1982–present) survey of all state wildlife management agencies conducted by the Southeastern Cooperative Wildlife Disease Study (SCWDS; Nettles et al. 1992; Ruder et al. 2016). This survey is completed annually by state wildlife agency personnel in all 50 states who are responsible for the conservation and management of deer populations in the US (Ruder et al. 2016). The survey uses syndromic surveillance to report HD based on characteristic field signs and postmortem lesions associated with HD in WTD populations, as well as diagnostic testing of select mortality events. Data are collected at the county level in a binary format of 1 or 0 to represent presence or absence of HD, respectively, based on one or more of the four criteria established at the inception of the survey in 1982. These criteria have remained consistent throughout the entire study and include: 1) sudden and unexplained deer mortality that occurs during late summer and early autumn; 2) necropsy-based diagnosis of HD based on lesions; 3) isolation or molecular detection of EHDV or BTV from postmortem tissues; and 4) observation of deer with sloughing hooves (Nettles et al. 1992; Ruder et al. 2016). Depending on the outbreak, a report of HD in a county might involve a small cluster of dead deer or potentially dozens to hundreds of dead deer. For our project, the presence of HD in counties was determined using the first three reporting criteria, to focus on HD mortality. The fourth criterion was excluded from the data set because it focuses more on HD morbidity and was less consistently used across all five states over the 38-yr study period.

Virus detection data (criterion 3) was also compiled for other wild ruminant species in the study area. These data, combined with criterion 3 for WTD, help to corroborate the spatial and temporal patterns of suspected HD (criteria 1 and 2) in WTD for the Great Plains between 1982 and 2020. Use of confirmatory diagnostic testing on all dead wild ruminants during outbreaks is not possible, but their use on a subset provides confidence in HD reports based on criteria 1 and 2 during outbreaks. Virus detections from postmortem tissues (e.g., spleen, lung, lymph node) were performed at SCWDS, as described (Kienzle et al. 2017), or at state or federal veterinary diagnostic laboratories, and reported in the annual HD survey. Reported virus detection data were organized by state, species, serogroup, and serotype, along with the number of years the specific serotypes were reported.

We estimated the potential change in reported HD patterns through data visualization, evaluation of intensity, and frequency of HD reports. To visually observe potential changes in reported frequency of HD trends in the Great Plains, we created county-level maps split into four decades: 1982–1991, 1992–2001, 2002–2011, and 2012–2020. We define a decade as a set of 9 yr, with the last set resulting in 8 yr, 2012–2020. Reporting frequency was defined as the number of years a county or state reported HD within each of these four decades.

For our temporal analysis, we used linear regression to estimate the change in HD reporting intensity from 1982-2020. We defined reporting intensity as the percentage of counties in a state with reported HD each year (i.e., number of counties with reported HD divided by total number of counties in the state × 100). Representing reporting intensity as a percentage provided a scale for all counties across the Great Plains states. Scaling the counties allowed us to compare the temporal change in HD reporting intensity from 1982 to 2020 for all five states without bias towards states that have more counties than others. For simplicity, we excluded virus detection data (criterion 3) from wild ruminant species other than WTD from statistical analyses, because the data did not provide new county detections in addition to counties where HD was reported in WTD populations. Therefore, all statistical analyses were conducted on WTD (deer) reports only. Finally, we used linear regression over all counties, regardless of state, to determine a statistical significance in change in HD reporting frequency from 1982 to 2020 in the Great Plains.

For spatial analysis, we used generalized linear mixed models (GLMM) with a binomial distribution to analyze changes in the probability of detecting HD for the Great Plains (Johnson et al. 2006). We developed univariate and multivariate models with independent variables of latitude, longitude, and year to assess how the dependent variable (HD report probability for a county each year) changed over time. We included interaction terms in the multivariate models to determine how the probability of detecting HD changed relative to other independent variables. To account for the numerical differences between the ranges of the independent variables of time and location, we centered and scaled latitude, longitude, and year by standard deviation prior to running models.

Model selection was not used for the temporal linear regression for each state, but we used Akaike information criterion corrected for small sample size (AIC_c) for model selection to determine the top model for spatial changes in HD patterns (Burnham and Anderson 2002). We used our top spatial model to determine how patterns of HD reports changed both for counties in the Great Plains (excluding state boundaries) and for each state. We ran the top model by each state to account for reporting differences that may have occurred between states between 1982 and 2020. In the model used for all counties in the Great Plains, state was used as a random effect to account for the variation in annual HD reports from all states (Gillies et al. 2006). All linear regression and GLMM models were run using the lme4 package in R version 4.0.3, and statistical significance was determined with P≤0.05 (Bates et al. 2015; R Core Team 2022).

Over the 38 yr of the study, all states submitted HD reports for each county regardless of HD status. In total, 1,262 county-level HD reports were submitted between 1982 and 2020 across all five states. Counties with the majority of HD reports followed a diagonal pattern from southeast to the northwest corners of the region, similar to larger HD patterns previously described by Stallknecht et al. (2015; Fig. 1A). The counties where EHDV and/or BTV was confirmed through diagnostic testing (criterion 3) during the entire study period are shown in Figure 1B, revealing those counties in the study area where virus confirmation was achieved throughout the distribution of suspected HD (criteria 1 and 2). Virus detections from WTD and other free-ranging ruminant species including mule deer, elk, bighorn sheep, and pronghorn, are shown in Table 1.

Changes in report frequency and report intensity varied for each state. Report frequency increased for all states throughout the decades, but Nebraska experienced a slight decrease in report frequency between the third and fourth decade (Table 2). Within the five states, Oklahoma had the lowest reporting frequency throughout all decades (Table 2). South Dakota had the highest reporting frequency in counties in the fourth decade. In addition to observed increased report frequency, our spatial map of reported HD by decade revealed the observed change in geographic distribution of reported HD for all states (Fig. 2).

Reporting intensity also varied between the states over the 38 yr. At the state level, North Dakota, South Dakota, and Kansas had a significant increase (P≤0.05) in HD reporting intensity, but reporting intensity in Nebraska and Oklahoma did not significantly increase from 1982 to 2020 (Fig. 3). South Dakota experienced the largest increase in report intensity, with a probability of 0.58, whereas Oklahoma experienced the lowest increase in reporting intensity, with a probability of 0.041 (Table 3). When all counties were combined, reporting intensity did significantly increase (P≤0.05) from 1982 to 2020 for the Great Plains region (Table 3). Overall, each state experienced at least 1 yr without any counties reporting HD, but North Dakota, South Dakota, Nebraska, and Kansas experienced peaks (≥25%) in the percent of counties reporting HD in some years, suggesting geographically large-scale outbreaks (e.g., 2012 in North Dakota, South Dakota, Nebraska, and Kansas; 2020 in North Dakota and South Dakota; Fig. 3).

Lastly, our top spatial model suggested that the probability of detecting HD through the county-level reports increased in both latitude and longitude from 1982 to 2020. Our top spatial model included an interaction with all three variables of year, latitude, and longitude (Supplementary Material Table S1). Latitude and the interaction of latitude and year were not important predictors of spatial change in HD reports, but longitude, year, the interactions of latitude and longitude, and the three-way interaction of year, latitude, and longitude were significant predictors of changing patterns of HD reports (Supplementary Material Table S2). We observed an increase in the probability of detecting HD through a report for every 1.2° increase in longitude but did not observe a significant change in the probability of detecting HD through reports in latitude (Fig. 4; Supplementary Material Table S2).

Hemorrhagic disease in the Great Plains generally followed the pattern of enzootic stability in the south transitioning to patterns of epizootic mortality moving north. Consistent with our observation of very few HD reports in Oklahoma, deer populations in Oklahoma have a high population immunity, evidenced by high EHDV or BTV antibody prevalence with minimal HD mortality (Kocan et al. 1987). Moving north, the percent of counties with reported HD increased in Kansas, especially in the most recent decade (Table 2). Similar to Oklahoma, a pattern of enzootic stability has been identified in western Kansas (Flacke et al. 2004). However, Kansas appears to be an epidemiological transition state, with evidence of enzootic stability in the western counties transitioning to epizootic mortality in the east (Flacke et al. 2004). In our study, most Kansas counties reporting HD were primarily in the eastern half of the state (Figs. 1 and 2). Patterns of epizootic mortality intensified moving north from Kansas into Nebraska where HD reports have remained relatively consistent from 1982 to 2020, although Nebraska experienced peaks in mortality similar to the other four states (Fig. 3). In South Dakota over the past decade there is evidence of increasing outbreak intensity, as the percentage of counties reporting HD increased from 11.2% to 32.2% (Table 2). Furthermore, we showed an increase in the geographic distribution of reported HD in South Dakota, which was corroborated by postmortem diagnostic data (i.e., virus detections), from 1982 to 2020 (Table 1).

Although frequency of reported HD increased over the four decades, the geographic distribution of HD also changed in the Great Plains. Unlike northern expansion of HD previously observed in states east of the Great Plains (Stallknecht et al. 2015; Christensen et al. 2020; Dorak et al. 2022), our findings indicated that HD has expanded longitudinally in an eastern direction from 1982 to 2020 (Fig. 4). Hemorrhagic disease has been reported in North Dakota during all four decades of the study period, therefore we could not observe significant change in reported HD relative to latitude alone, but both visual observation and linear regression suggest that counties reporting HD have expanded eastward and northward across the Great Plains (Fig. 2).

While not investigated in our study, the overall shift in the distribution, intensity, and frequency of reported HD for the Great Plains may be related to recent changes in climate conditions (i.e., increasing annual temperatures) and subsequent impact on vector populations and virus transmission (Purse et al. 2005). Temperature impacts vectorial capacity by dictating vector developmental rates, biting rates, and virogenesis (Mullens et al. 2004). As temperature increases between 27 C and 30 C, the Culicoides generation time and extrinsic incubation period decreases, resulting in more infected adult Culicoides on the landscape (Wittman et al. 2002). In regions with patterns of epizootic mortality such as North Dakota and South Dakota, changing climatic variables may enhance vectorial capacity (Sleeman et al. 2009; Salley et al. 2016). Although we acknowledge that certain climate components such as temperature and moisture influence vector-virus interactions, further research is warranted to understand how habitat in relation to climate impacts HD reports in the Great Plains.

A noteworthy change in HD reporting intensity and frequency in the Great Plains occurred between South Dakota and North Dakota. In the northeastern US, the effect of latitude was significant in predicting HD mortality relative to drought index, but reported HD generally decreased north of 40° N latitude (Stallknecht et al. 2015; Christensen et al. 2020). Similar to the trend of reduced HD reporting frequency at northern latitudes, HD reporting frequency and intensity increased at a smaller rate in North Dakota versus South Dakota, potentially related to the effect of local climate conditions on vectorial capacity. Populations of C. sonorensis (a confirmed vector of EHDV and BTV), were found during 2002 in regions of western North Dakota, western South Dakota, and Nebraska where soils are nonglaciated and annual evaporation exceeds annual precipitation (Schmidtmann et al. 2011). Whether C. sonorensis populations have since expanded into areas east of the Missouri River, or whether other resident Culicoides spp. are responsible for transmission of EHDV and BTV in eastern North Dakota and South Dakota remains unclear because of the paucity of Culicoides surveillance in the region. Other Culicoides spp. that are suspected vectors of EHDV and BTV and could play a role in the region include Culicoides stellifer, Culicoides paraenesis, Culicoides venustus, and Culicoides haematopotus, among others (McGregor et al. 2022). Regardless of the underlying factors involved, the differences in reported HD intensity and frequency between South Dakota to North Dakota might demarcate the northeastern boundary of HD in the Great Plains and changing climatic conditions might continue to shift the leading edge of HD in a northeasterly direction.

As climate continues to change, vector–virus interactions and their impact on HD morbidity and mortality in northern deer populations will also change. Continued surveillance is essential for monitoring such impacts, understanding their drivers, and anticipating implications for deer management. One notable HD outbreak occurred in 2012 across most of the Midwest, impacting both free-ranging and domestic ruminant species (Stevens et al. 2015; Ruder et al. 2015). Data from our study clearly show prominent peaks of HD reports during 2012, in Kansas, Nebraska, and South Dakota, consistent with that outbreak (Fig. 3). The 2012 Midwest HD outbreak was predominantly attributed to EHDV-2, and the severe drought conditions and high temperatures during the summer months probably contributed to the intensity of the outbreak (Stevens et al. 2015). Given that the Great Plains encompasses a gradient of different climate conditions, considerable changes in climate might affect how vector-borne diseases, such as HD, present on the landscape. Thus, it is plausible that HD may become an increasingly common mortality factor among some deer populations.

One important aspect of the complex HD system is host diversity and how different hosts for EHDV or BTV affect HD distribution on the landscape. Although WTD are considered a primary host for EHDV and BTV in North America, other free-living species such as mule deer, elk, and pronghorn also are frequently involved in outbreaks in the Great Plains (Howerth et al. 2001). Furthermore, livestock hosts for EHDV and BTV, such as cattle and sheep, are common in the Great Plains, are often sympatric with wild ruminants, and are known to affect populations of certain Culicoides species (Hardy and Price 1952; Ruder et al. 2015; McGregor et al. 2022). Although wild ruminant species often suffer severe or fatal disease, cattle infected with EHDV and BTV typically develop inapparent or mild disease. Both domestic and wild ruminant species are vertebrate hosts defining this disease system, and the presence and movement of susceptible livestock species may potentially influence the occurrence and movement of EHDV and BTV on the landscape (Rivera et al. 2021). To understand the complexities of this vector-borne disease system involving multiple ruminant hosts and multiple viruses, future research and surveillance efforts should incorporate both domestic and free-living ruminant hosts, as well as Culicoides spp. vectors. Specifically, long-term studies that spatially and temporally combine investigations of free-living ruminant mortality; serological testing of harvested free-living ruminants; serologic and virologic testing of sentinel livestock; and Culicoides spp. surveillance across diverse land and animal use settings, would collectively help inform the understanding of orbivirus epidemiology in the Great Plains and elsewhere.

To our knowledge, no other HD data set exists for a similar length of time and spatial coverage. The HD survey relies on the continued participation of state wildlife management agency personnel, but we acknowledge certain levels of bias exist within the data set. Disease events in wildlife may be difficult to observe, thus potentially leading to underreporting of events such as HD. The HD survey primarily captures outbreaks with mortality, but small-scale mortality in dense vegetation or where humans are not present, for example, may not be detected. In enzootically stable regions, virus transmission events do not result in significant morbidity or mortality, and these inapparent infections are therefore not accounted for in this survey (e.g., in Oklahoma). Further, reports of HD were received across the 38-yr period from independent sources that could not be standardized between states (Stallknecht et al. 2015). Nevertheless, reliability of these data is made possible by the diligence and continued support of state agency personnel, well-characterized and consistent field signs and necropsy lesions of HD in WTD, consistent criteria for reporting HD by state agencies, the relative visibility of deer to the public, and the value of deer as a natural resource (Ruder et al. 2016). The validity of our survey is further supported by lack of outlier reports in areas with patterns of epizootic mortality (e.g., northeast US); consistency in spatial and temporal reporting of HD during outbreaks regardless of state borders; and through confirmatory diagnostics throughout the study (criterion 3; Stallknecht et al. 2015). Given the duration and consistent case criteria, this survey provides unique opportunities to study spatial and temporal changes in patterns of HD that short term projects could not capture.

The general pattern of increase in the frequency and intensity of reported HD in parts of the central and northern Great Plains, and the significant latitude and longitude interaction supports the concept that conditions have changed in favor of HD outbreaks from 1982 to 2020. Additionally, our finding that counties reporting HD expanded eastward in North Dakota and South Dakota suggests that this region may represent the leading edge of HD in the northern Great Plains. Continued monitoring of EHDV and BTV in the multihost system of the Great Plains is important to document future changes in the face of changing climate.

This study was made possible through long-term financial support from state member wildlife agencies of Southeastern Cooperative Wildlife Disease Study (SCWDS) provided by the Federal Aid to Wildlife Restoration Act (50 Stat. 917). We also thank the US Fish and Wildlife Service National Wildlife Refuge System and the US Geological Survey Ecosystems Mission Area for their long-term support of SCWDS. Furthermore, this research depends on the continuous efforts of the numerous wildlife biologists, veterinarians, and technicians from state wildlife agencies that contribute to this long-term survey and provide diagnostic samples for testing. We are especially grateful to current and former staff with the North Dakota Game and Fish Department, South Dakota Game, Fish & Parks, Nebraska Game and Parks Commission, Kansas Department of Wildlife & Parks, and Oklahoma Department of Wildlife Conservation for providing data through the annual survey. We thank D. Barber, S. Hesting, T. Nordeen, and S. Griffin for reviewing this manuscript. Additional virus detections were made by the US Department of Agriculture’s National Veterinary Services Laboratories and the South Dakota State University Animal Disease Research and Diagnostic Laboratory. Laboratory and data management support at SCWDS were provided by V. Nettles, J. Fischer, R. Poulson, C. Kienzle, N. Stilwell, A. Allison, S. Vigil, and C. McElwee.

Supplementary material for this article is online at http://dx.doi.org/10.7589/JWD-D-23-00021.