GENETIC RELATEDNESS OF EPIZOOTIC HEMORRHAGIC DISEASE VIRUS SEROTYPE 2 FROM 2012 OUTBREAK IN THE USA Open Access

Two closely related orbiviruses (family Reoviridae), epizootic hemorrhagic disease virus (EHDV) and bluetongue virus (BTV) cause hemorrhagic disease (HD), a leading cause of morbidity and mortality in white-tailed deer (WTD; Odocoileus virginianus) throughout North America. These segmented, double-stranded RNA viruses are transmitted by biting midges from the genus Culicoides (Howerth et al. 2001). While both EHDV and BTV affect WTD, mortality events in the eastern US are most often associated with EHDV (Stallknecht et al. 2015).

The EHDV serotypes 1 and 2 (EHDV-1 and -2) are distributed broadly across the US; a third serotype, EHDV-6, emerged as a reassortant between an exotic EHDV-6 (CSIRO-753–like) and North American EHDV-2 (Alberta) strains (Allison et al. 2012). Since its initial detection in 2006, EHDV-6 has become endemic in the US (Allison et al. 2010; Stallknecht et al. 2015). Serotype 2 is found in Asia, Australia, and North America. A closely related EHDV-2, called strain Ibaraki, circulates in the eastern hemisphere throughout Japan and Australia. Strain Ibaraki causes disease in cattle, whereas North American EHDV-2 viruses primarily affect deer (Inaba 1975; Howerth et al. 2001; Anthony et al. 2009b). North American EHDV-2 viruses dominated recent large-scale outbreaks in 1988, 1996, 2007, and 2012 (Stallknecht et al. 2015). During the 2012 outbreak, 66% (135/205) of the EHDV and BTV viruses isolated from animals in 21 states were identified as EHDV-2. During this outbreak, previously undocumented levels of clinical disease were also observed in domestic ruminants from confirmed EHDV infection (Stevens et al. 2015).

Viral overwintering mechanisms and outbreak dynamics of HD are poorly understood. Hemorrhagic disease is seasonal, with disease typically documented during summer through late fall, subsiding shortly after periods of colder weather that negatively affect vector abundance and activity (Nettles and Stallknecht 1992). There are two predominant theories as to how these viruses are maintained during winter and spread during regional outbreaks. One theory is that a single virus strain is introduced from areas where overwintering conditions are favorable and subsequently spreads during an outbreak. Another theory is that multiple viruses expand from local reservoirs that sustain virus transmission or viability, either in the ruminant or insect host, so that multiple strains of virus comprise the outbreak. If the former is correct, one single genotype of virus would be found throughout the outbreak. If the latter, subtle genetic diversity, such as a collection of single-nucleotide polymorphisms (SNPs) forming distinct genotypes, would be detected. The involvement of multiple viral strains within a single outbreak could result from the movement of infected vectors or animals from enzootic-subclinical regions in the southeastern US (Sellers and Mellor 1993; Koenraadt et al. 2014). Alternatively, multiple viral strains could expand from cryptic refuges within the outbreak area (Murphy et al. 2005; Maclachlan 2011).

Because of the widespread nature of EHDV-2 during the 2012 outbreak, an unprecedented opportunity emerged to examine the relatedness of geographically distinct isolates representing a single serotype during a single epidemic. Because many of the genes composing EHDV show very little variation, we examined two genes that encode viral proteins under the most selection pressure in the vector-virus-host system. The outer capsid protein, VP2, is both the mammalian receptor binding protein and the antigen determining protein; this gene experiences selection pressure from the ruminant host (Oldfield et al. 1991; Maan et al. 2007). The outer protein of the inner capsid, VP7, is the insect receptor binding protein (Xu et al. 1997; Tan et al. 2001). Although VP7 is functionally constrained by protein-protein interactions within the virion, this gene potentially receives selection pressure from the vector (Tan et al. 2001; Roy 2008).

Existing studies utilizing sequence analyses have failed to shed light on EHDV overwintering mechanisms. The reasons for this are varied, but generally relate to geographic and temporal scale, as well as small sample sizes. Most studies have focused on questions related to virus evolution or the detection of topotypes and include multiple EHDV serotypes over multiple years or examine just one protein or different proteins of interest (Cheney et al. 1995; Mecham et al. 2003; Anthony et al. 2009a, b). Additionally, full-length gene or segment sequences are not always evaluated. A large number of EHDV-2 isolates from the 2012 HD outbreak provided a unique opportunity to focus on a single EHDV serotype during a single year over a wide geographic area in North America. The objective of our study was to determine if a single strain or multiple EHDV-2 strains comprised the outbreak. We hypothesized that genetic diversity would reflect source virus diversity, with minimal diversity associated with the expansion of an introduced strain and increased diversity and geographic clustering associated with local expansion of endemic foci.

Viruses utilized for this study were isolated from diagnostic samples submitted for EHDV and BTV testing during a large-scale 2012 HD outbreak in the US (Fig. 1; Abdy et al. 1999; Allison et al. 2010). From July to November 2012, 135 EHDV-2 viruses were isolated from WTD (n=124), cattle (Bos taurus; n=9), mule deer (Odocoileus hemionus; n=1), and an alpaca (Vicugna pacos; n=1). Most of the samples were from wild populations; 15% (n=19) of WTD were captive individuals. A geographically representative subset of these viruses was selected for sequencing, including at least one isolate from each of the 21 states where EHDV-2 was isolated (Table 1 and Fig. 1). To eliminate the variable of host difference, isolates were chosen from genus Odocoileus. To provide an internal control to confirm that viruses circulating in close proximity are genetically conserved, we compared sequences from two isolates (CC12-351, CC12-352) from a single county (Greene County, Missouri).

Viral RNA was extracted from low passage cultures using QIAamp^® Viral RNA Mini kit (Qiagen, Hilden, Germany) and frozen at –20 C until use. Viral RNA for VP2 and VP7 were amplified using sequencing primers and SuperScript^® III One-Step RT-PCR Platinum^® Taq High Fidelity (Invitrogen, Life Technologies, Carlsbad, California, USA). Sequencing primers were designed using alignments of full-length gene sequences and are available upon request. Each 50-lL single-tube reverse transcription-PCR was assembled according to SuperScript protocol using the following cycling parameters: complementary DNA generation at 45 C for 30 min and 2 min at 94 C followed immediately by 40 PCR cycles of 94 C for 15 s, 47 C for 30 s, and 68 C for 1 min, with a final extension cycle of 68 C for 5 min.

Amplicons were separated on a 1.5% agarose gel stained with ethidium bromide, excised, and gel purified using QIAquick^® Spin kit (Qiagen, Germantown, Maryland, USA). Purified complementary DNA was stored at –20 C until sequenced. BigDye^® Terminator v3.1 Sequencing kit (Life Technologies) reactions were set up according to protocol and purified by Sephadex column. Sanger sequencing was performed by either Georgia Genomics Facility (Athens, Georgia, USA) or Macrogen USA (Rockville, Maryland, USA).

All sequencing reads were trimmed and aligned to references HM636898 and HM636903 for VP2 and VP7, respectively, using Geneious version 8.1.8 (Kearse et al. 2012). Consensus sequences were trimmed to the open reading frame (ORF) AUG to TAG for VP2 and AUG to TAC for VP7. Gene sequences from this study are listed in GenBank under accession numbers MH347379–MH347412 and MH347413–MH347447 for VP2 and VP7 genes, respectively. Gene sequences from 2012 isolates were aligned using ClustalW, MUSCLE, and MAFFT plugins for Geneious v8.1.8 (Edgar 2004; Larkin et al. 2007; Katoh and Standley 2013). Maximum likelihood (ML) trees for VP2 were built using 1,000 bootstraps and the TN93+I model from the representative MUSCLE alignment, as determined by best model fit test in MEGA v 7.0.16 (Kumar et al. 2016). Additional full-length gene sequences for VP2 from GenBank were included for analysis. Maximum likelihood trees of VP7 gene sequences were built in MEGA using the T92+G+I algorithm, as determined by best model fit test in MEGA, and performed using 1,000 bootstraps based on the MUSCLE alignment (Nei and Kumar 2000; Kumar et al. 2016). As with the large VP2 trees, other full-length VP7 sequences from EHDV serotypes were included in the ML tree with 1,000 bootstraps. Each protein was analyzed using Tajima's neutrality test, a test of homogeneity of substitution patterns between sequences, and the ML estimate of transition/transversion bias in MEGA (Tajima 1989; Nei and Kumar 2000; Kumar and Gadagkar 2001). Additionally, a Mantel matrix comparison test was performed in XLSTAT using the VP2 pairwise percent nucleotide identity matrix and a matching distance matrix created from the central point from each county of isolation because of a lack of geolocation data for each sample (Mantel 1967; XLSTAT 2016).

Complete ORF sequences of two genes were obtained from 34 different EHDV-2 viruses isolated during the 2012 outbreak. Of the 34 viruses examined by pairwise comparison, 29 were unique at the nucleotide level for VP2 genes and 19 isolates were unique at the nucleotide level for VP7.

Sequences for the mammalian binding protein (VP2) varied considerably more than VP7 sequences within this data set. Average pairwise identity of VP2 sequences from the 2012 isolates was 99.1%, ranging from 96.6–100%; average amino acid identity was 99.3%. The ORF of the VP2 gene was 2,949 bases long; of those, 2,772 bases (94%) remained fixed; 49 amino acid changes were observed in 983 residues (95% remained fixed). Gene sequences of VP2 also formed two main clades and a group of five unresolved sequences in an unrooted tree (Fig. 2), but there were considerably more polymorphic sites than seen in VP7 (see Supplementary Material Tables S1, S2). At least 11 SNPs distributed across the gene defined the mid-Atlantic group, as seen in Table S2. Within the western group (Fig. 2) was a mixture of isolates from the core of the outbreak (Missouri, Kansas) and from potential outliers not connected to the main body of the outbreak (Montana, Wyoming, Utah, Colorado). In the mid-Atlantic clade, however, the branch lengths were much shorter, indicating a higher similarity between viruses in the southeastern and mid-Atlantic states (Georgia, North Carolina, Virginia, Maryland, West Virginia, New Jersey, Pennsylvania, Kentucky). Located within the clades possessing the mid-Atlantic genotype, isolates from Indiana, Iowa, and Kansas were also found. In Figure 3, the EHDV-1 (New Jersey)–rooted ML tree of 2012 isolates, including Ibaraki and two other EHDV-2 sequences, had 487 variable bases in the 2,762 base pairs ORF alignment. The model of best fit was determined to be TN93+I+G, where the +I parameter, representing the number of invariable sites, used in the trees was 26%. In VP2, the transition/transversion ratio was 3.7. When rooted by EHDV-1 (New Jersey), much of the distinction between clades found in the unrooted tree disappeared, demonstrating the relative relatedness of the 2012 VP2 gene sequences. In the amino acid sequence tree, no clade distinction was made between any of the North American EHDV-2 isolates. The Mantel test of matrix comparison between the pairwise percent nucleotide identity matrix of VP2 clades and the corresponding geographic distance matrix yielded a statistically significant association (P<0.001) between the two matrices. Although the association was not strongly linear, the correlation value of –0.67 indicated a negative association between genetic similarity and geographic distance within the two main clades created from the 2012 analysis. Resolution between VP2 sequences from 2012 further decreased when compared with historical isolates and other serotypes with worldwide distribution (Fig. 3A). Each serotype fell out in distinct clades, visually representing relationships between global distance and temporal clustering. Tajima's neutrality test on the nucleotide alignment of VP2 produced a D statistic of –2.54, a negative result greater than two standard deviations from a null hypothesis of neutral evolution. This is indicative of recent population expansion after a bottleneck event (Tajima 1989; Nei and Kumar 2000).

Average pairwise identity of VP7 gene sequences within the group of EHDV-2 isolates from 2012 was 99.1%, ranging from 95.8% to 100%. At the amino acid level, the isolates were so conserved that only three of the 34 VP7 gene sequences were unique. In the ORF of 1,050 bases, 990 of those bases (94%) remained fixed. Of the 60 bases across the entire reading frame that differed, 58/60 (97%) were silent mutations. Only two residues were changed across the entire protein sequence. Both residue mutations were conservative amino acid substitutions (N-275-S and M-349-L). Of the 2012 isolates, two distinct clades are present within the VP7 tree (Fig. 4). The clades fell into two topotypes, based on two SNPs at nucleotides 97 and 561 (Table S2). Subclades are further distinguished by additional SNPs. The midwestern/western US clade has a genotype consisting of thymine at nucleotide 97, guanine at nucleotide 285, and cytosine at nucleotide 561, whereas isolates predominantly from the mid-Atlantic region of the US formed another clade based on a genotype consisting of cytosine at nucleotide 97 and thymine at nucleotide 561. When additional VP7 gene sequences were added, EHDV-2 still clusters, but American EHDV-6 VP7 sequences were interspersed throughout the Alberta-type North American EHDV-2, instead of clustering together (Fig. 3B). Maximum likelihood estimate of transition/transversion bias (R) was found to be 2.57, such that transitional changes that most often are silent mutations are 2.57 times more likely to occur within VP7. Approximately 62% of the gene (651/1,049 bases) was invariable (+I) within this data set. When the best model fit test (MEGA) was performed on the larger VP7 data, the model remained the same, but the proportion of invariable sites (+I) dropped to 58%.

Isolates WDL12-286, CC12-297, and CC12-463 did not fall into a resolved clade in either protein. They more closely resembled the historic isolate and reference sequence (CC126-00 from Georgia in 2000) that we used for alignment.

The internal control isolates from the same county (CC12-351 and -352) were 99.8% identical at the nucleotide level in both proteins. In VP2, there were six nucleotide differences, two of which were nonsynonymous mutations. Two nucleotide differences were present in VP7 but were both synonymous mutations causing no difference in amino acid sequence.

We expected a high level of nucleotide and amino acid identity within both gene sequences based on documented sequence conservation of EHDV viruses isolated several years apart; however, we did not expect the degree of sequence similarity we found (Mecham et al. 2003; Murphy et al. 2005). Although RNA viruses, orbiviruses are less diverse than other, better known RNA viruses such as influenza or coronaviruses, but they possess the genetic conservation typically found in other arboviruses (Weaver 2006). Orbiviruses are thought to be slightly more stable because of double-stranded genomes but are also evolutionarily limited by the vector-host transmission cycle and the double-filter effect found within arboviruses (Deyde et al. 2006). The high level of genetic similarity between gene sequences examined and the degree of robustness is consistent with a period of evolutionary stasis in which the virus is very well adapted to the vector-host system (Lauring et al. 2013). Rotaviruses (also subfamily Sedoreovirinae), although not arboviruses, have point mutations, through which lineages and sublineages within serotypes can be traced, as similarly found in the EHDV-2 results of the present study (Desselberger et al. 2001; Dennis et al. 2014).

In the alignments of both VP2 and VP7 nucleotide sequences, several point mutations were associated with the geographic origin of isolates, suggesting regional clustering. With VP2, polymorphisms at 26 nucleotide locations corresponded with division into two major clades (Table S2 and Fig. 2A). Additional nucleotide diversity within those 26 locations provided more distinct genotypes. For example, guanine instead of adenine at positions 528 and 987, and adenine instead of guanine at positions 1191 and 1512 separate the tightly clustered mid-Atlantic clade, which includes isolates from Georgia to New Jersey, from other nearby isolates in Indiana, Pennsylvania, Tennessee, and West Virginia. In the full alignment, however, the Franklin County, Kentucky; Lancaster County, Nebraska; Tama County, Iowa; Blaine County, Montana; and Barry County, Michigan, isolates have many more polymorphic sites than the rest of the 2012 isolates and share more identity with a historical isolate, CC126-00, isolated from a WTD in Georgia in 2000 (HM636898 and HM636903). In the VP2 alignment, more variable nucleotide sites were found within a distinct region; this region corresponds to the putative epitope binding region between nucleotides 480 and 1110 of the VP2 gene (Gould and Eaton 1990; Pritchard and Gould 1995; DeMaula et al. 2000; Murphy et al. 2005). In the summarized alignment of VP7 polymorphisms, only six nucleotide locations across the 1,050 base pair ORFs were associated with the tree topology (Table S2). Isolates from Colorado and Missouri stand out as having guanine and thymine rather than adenine at positions 603 and 1045, respectively. In the VP7 gene sequences, there was no distinction within the “eastern” clades, as seen in the tree, which differs from the VP2 gene sequences (Figs. 2, 4).

When examining nucleotide differences, several subtle but distinct regional genotypes emerged within both VP2 and VP7, supporting the idea that multiple strains from endemic foci expand during outbreaks (Wilson et al. 2008). These regional clusters corroborate findings from similar North American analyses using temporally and geographically distributed EHDV-2 sequences (Cheney et al. 1996; Murphy et al. 2005). Although it is possible that multiple introduction events could have provided the multiple viral strains seen during this outbreak, the slight difference found in the eastern clade, which creates an intermediate genotype and suggests evolution within the greater eastern clade, would be consistent with long-term evolution and circulation. Certain geographic and geologic factors such as landscape cover, changes in elevation, and temperature changes have been found to define boundaries to BTV spread in Europe and Australia (Bishop et al. 2000; Pioz et al. 2012). The greater eastern clade does not stop at the Appalachian Mountains; it extends west through the Appalachian Highlands to the largely agricultural land throughout the Midwest and Interior Plains, with a border somewhere near the Ohio River.

Other aspects of our data reinforce the intricate nature of the vector-virus-host system and warn against oversimplified conclusions. First, the two genes examined did not match exact geographic patterns, suggesting that they had undergone differing selection pressures. Also, there were isolates located farther away than others within geographic clusters of genetically similar, if not identical, isolates, which accounted for the imperfect relationship between genetic and geographic distance found in the Mantel distribution. This could have several causes. For example, many species of wild and domestic ruminants can develop sufficient viremia for vector transmission. Therefore, interstate movement of infected domestic animals (including captive deer) or vectors cannot be ruled out as a possible method of EHDV-2 introduction (Sellers and Pedgley 1985). Sequence analysis of paired historical isolates from the same locations should provide future insight into local maintenance of EHDV-2 strains. In BTV, reassortment plays a much larger role in intertypic diversity that historically has not been documented in EHDV (He et al. 2010; Nomikou et al. 2015; Boyce et al. 2016). Since at least 2006, EHDV serotypes 2 and 6 can, and have, undergone reassortment, thus providing another method for the two genes to have differing inheritance patterns (Allison et al. 2012). Because 33% of all of our EHDV isolations in 2012 were EHDV-6, examination of VP2 and VP7 proteins in those isolates would also help to elucidate reassortment mechanisms in North America.

Putting our 2012 data in context with sequences from all over the world and spanning all seven EHDV serotypes in Figure 3, little resolution within the various serotypes was observed. In the VP2 tree (Fig. 3A), a trend toward temporal association (year of isolation) was seen on a very small scale within EHDV-2 viruses, but a few sequences in 2012 shared higher nucleotide identity with the historic isolate CC126-00, from a WTD in Georgia during 2000 (Fig. 2). These viruses were primarily located on the outskirts of the main body of the outbreak, as reported in Figure 1, but the Tama County, Iowa, isolate CC12-463 fell well within the outbreak. Another outlier sequence from Barry County, Michigan, appeared within different clades in both VP2 and VP7 trees. Michigan experienced its first large outbreak of EHDV-6 in 2012, which could suggest another introduction of diversity within the EHDV-2 viruses. We lacked a sufficient sample size to infer a temporal relationship based on year of isolation outside of the 2012 samples. Our data was also insufficient to find any pattern relating the date of isolation and genetic similarity on a monthly or weekly scale, which would have allowed us to trace the passage of virus across the outbreak. Higher variation in VP2 gene sequences was to be expected, as was the phenotypic conservation of VP7 (Mecham et al. 2003; Maan et al. 2007; Anthony et al. 2009a, b). The strongly negative Tajima neutrality test statistic of VP2 also suggested that the differences found between clades of VP2 are from different viral origins. The negative D value (D=–2.54) obtained using Tajima's neutrality test on the nucleotide alignment corresponded with lower than expected diversity given the number of variable sites.

The VP7 gene and protein has the largest difference between nucleotide identity and amino acid identity of all viral proteins in both EHDV and BTV, such that high levels of synonymous mutations should be expected (Anthony et al. 2009a). This protein is also fairly conserved between serotypes and is commonly used as a target for EHDV and BTV differentiation. Despite expecting high conservation, there was variation in the VP7 gene sequences from the county (Greene County, Missouri) in the center of the outbreak. Given the high level of conservation between geographically distinct isolate sequences, this result was surprising and perhaps indicative of a greater degree of variation in a smaller community that is not seen in a larger scale genetic survey of an outbreak. One of the two isolates, CC12-351 from Greene County, was sequenced elsewhere using NexGen sequencing and matched our sequence, resolving any doubt about contamination or methodology (Wilson et al. 2016). While little variation was observed in VP7, topotypes are present and North American EHDV-6 VP7 genes cannot be distinguished from EHDV-2 VP7 genes in a phylogenetic tree (Fig. 3B). This suggested continuing reassortment between North American EHDV-6 viruses and the more diverse group of cocirculating EHDV-2 viruses, although EHDV reassortment is far less frequent than within BTV serotypes (He et al. 2010; Boyce et al. 2016).

Accumulation of SNPs forming distinct genotypes with regional clusters is consistent with local virus maintenance and emergence from endemic foci during advantageous vector, virus, host, and environmental conditions. Although only two genes were sequenced, complete sequences detected additional subtleties that would have been missed if partial sequences were used. The location of these two proteins in separate pieces of the viral capsid also allowed us to detect potential reassortment of the currently expanding EHDV-6 North American strain. Sequencing of additional historic EHDV-2 isolates and recent EHDV-6 isolates from the same location would elucidate those mechanisms further. In conclusion, presence of distinct regional genotypes within both genes examined is consistent with multiple strains of EHDV-2 viruses composing the 2012 outbreak.

This study was made possible by the continued support from member agencies of the Southeastern Cooperative Wildlife Disease Study. A special thanks goes to the many biologists and technicians of various agencies for collecting and submitting tissues for HD surveillance and virus isolation. Additional financial support came from a Kansas Bioscience Authority funded grant to the Kansas State University Center of Excellence for Emerging and Zoonotic Animal Diseases.

Supplementary material for this article is online at http://dx.doi.org/10.7589/2017-05-125.