Resolving Bovine viral diarrhea virus subtypes from persistently infected U.S. beef calves with complete genome sequence

Abstract

Bovine viral diarrhea virus (BVDV) is classified into 2 genotypes, BVDV-1 and BVDV-2, each of which contains distinct subtypes with genetic and antigenic variation. To effectively control BVDV by vaccination, it is important to know which subtypes of the virus are circulating and how their prevalence is changing over time. Accordingly, the purpose of our study was to estimate the current prevalence and diversity of BVDV subtypes from persistently infected (PI) beef calves in the central United States. Phylogenetic analysis of the 5′-UTR (5′ untranslated region) for 119 virus strains revealed that a majority (82%) belonged to genotype 1b, and the remaining strains were distributed between genotypes 1a (9%) and 2 (8%); however, BVDV-2 subtypes could not be confidently resolved. Therefore, to better define the variability of U.S. BVDV isolates and further investigate the division of BVDV-2 isolates into subtypes, complete genome sequences were obtained for these isolates as well as representatives of BVDV-1a and -1b. Phylogenetic analyses of the complete coding sequence provided more conclusive genetic classification and revealed that U.S. BVDV-2 isolates belong to at least 3 distinct genetic groups that are statistically supported by both complete and individual coding gene analyses. These results show that a more complex set of BVDV-2 subtypes has been circulating in this region than was previously thought.

Keywords

molecular epidemiology persistent infection pestivirus

Introduction

Bovine viral diarrhea virus (BVDV; family Flaviviridae, genus Pestivirus)⁷ is an important pathogen of cattle that causes respiratory disease, enteritis, and immune dysfunction.² BVDV can also cross the placenta and infect the fetus, resulting in abortion, congenital malformations, or birth of persistently infected (PI) calves.² Persistent infection develops when a susceptible dam is infected with a noncytopathic strain of BVDV and becomes viremic early in gestation, infecting the fetus prior to the maturation of its immune system.²⁵ Consequently, PI calves are born BVDV antibody negative, maintain a lifelong viremia, and continuously shed virus into the environment, making them the main reservoir for the virus and the leading source of transmission to susceptible cattle.^2,9,12 As a result, herds with PI animals have higher production costs because of poor feed conversion, increased illness, and higher mortality rates.¹⁷

BVDV has a single-stranded RNA genome that is ~12.3 kb and consists of a single open reading frame flanked by 5′ and 3′ untranslated regions (UTRs). Based on sequence comparison of the highly conserved 5′-UTR, BVDV isolates have been divided into 2 genotypes, BVDV-1 and BVDV-2,^28,32 each of which contains distinct subtypes^8,41 with genetic and antigenic¹ variation. Although both genotypes are diagnosed worldwide, the prevalence of subtypes varies by geographic region. Therefore, to effectively control BVDV by vaccination, it is important to know which subtypes of the virus are circulating and how their prevalence is changing over time. To date, 3 major subtypes reportedly circulate in the United States: BVDV-1a, -1b, and -2a^9,10; however, the predominance of these subtypes has shifted from BVDV-1a to BVDV-1b over the past 20 years in BVDV strains characterized by diagnostic laboratories.³⁴ In addition, 3 U.S. isolates of BVDV were classified in 2014 as BVDV-2c based on complete genome sequencing.¹⁹ There has also been a single case of BVDV-2b infection in the United States,¹⁴ highlighting the need for continued surveillance to detect emerging viral subtypes. Given the importance of PI cattle in the transmission of BVDV, the purpose of our study was 2-fold: 1) compare the current regional prevalence of BVDV subtypes in PI beef calves with similar samples collected 10 years earlier and 2) investigate the genetic variation among these viruses with complete genome sequencing.

Materials and methods

Study population and sample collection

Calves sampled in our study were purchased by order buyers from at least 5 central states (Alabama, Mississippi, Missouri, Oklahoma, and Texas) and shipped to a southwest Kansas feedyard between August 2013 and April 2014. On feedyard arrival, all calves were tested for BVDV infection by Central States Testing (Sublette, KS) using a noncommercial antigen enzyme-linked immunosorbent assay (ELISA) on supernatants from ear notch specimens. Cattle that tested positive for BVDV were classified as PI and maintained in pens apart from the other animals. For this study, blood samples collected with EDTA as anticoagulant and ear notches were collected from 131 potential PI calves on April 30, 2014. Ear notches were sent to the University of Nebraska–Lincoln Veterinary Diagnostic Center for verification of BVDV PI status by immunohistochemistry (IHC). PI status was independently verified by reverse transcription (RT)-PCR using pestivirus-specific primers and cycling conditions listed below.

Cattle in this study were privately owned by a commercial feeding facility, and the owners and management approved the use of animals for this study. Animal welfare for the entire facility was in accordance with the National Cattlemen’s Beef Association’s Beef Quality Assurance Feedyard Welfare Assessment program (http://www.bqa.org/), and experimental protocols were approved by the University of Nebraska Institutional Animal Care and Use Committee (approval 1172).

Genetic typing of BVDV

BVDV was genotyped directly from noncultured clinical samples to avoid unnecessary handling and potential for contamination. Total RNA was extracted from 250 μL of plasma using a phenol and guanidine isothiocyanate reagent.^a Pestivirus-specific primers (F: 5′-CATGCCCATAGTAGGAC-3′ and R: 5′-CCATGTGCCATGTACAG-3′¹³) were used to amplify a 283-bp fragment of the 5’-UTR.^b Cyclic amplification reactions were carried out in a 25-μL reaction containing 2.5 mM MgCl₂, 400 μM concentration of each dNTP, 0.4 μM concentration of each primer, and 1 μL of enzyme mix.^b Cycling conditions were as follows: reverse transcription for 30 min at 50°C, inactivation of RT enzyme and activation of hot start polymerase for 15 min at 95°C followed by 40 cycles of 94°C for 30 s, 60°C for 60 s, and 72°C for 60 s. A final extension was carried out at 72°C for 7 min. Following amplification, the 283-bp fragment of the 5′-UTR was confirmed using agarose gel electrophoresis. These fragments were purified by exonuclease treatment,^c according to the manufacturer’s specifications. DNA was then precipitated with ethanol and sequenced using the above pestivirus-specific primers. Sequencing was performed using dye terminator chemistry^d and resolved on a capillary DNA analyzer.^e

A neighbor-joining tree based on Hasegawa–Kishino–Yano (HKY) distance was produced using multiple sequence comparison by long-expectation (MUSCLE) algorithm for alignment of the 5′-UTR sequences using a commercial DNA analysis software package.^f For construction of maximum likelihood trees, molecular evolutionary analysis software^g was used to determine the most appropriate nucleotide substitution model using the best-fit model tool. The appropriate model was then used for phylogenetic inference according to the maximum likelihood algorithm using PhyML methods.^f,15 Support for tree nodes was determined by 1,000 bootstrap analyses. Classification of subtypes was based on reference sequences from GenBank as reported previously.^8,19,41

Generation of 19 full-length BVDV genomes

Nineteen full-length BVDV genomes were assembled de novo by sequencing virus strains directly from the plasma of infected calves as described previously.⁴³ This culture-free method was used to avoid the potential selection of a subpopulation of the infecting virus quasi-species when grown in cultured cells. Briefly, 100 ng of total plasma RNA were used as input material for an RNA library preparation kit.^h Libraries were prepared as specified by the manufacturer’s protocol with 1 important modification: BVDV genomes lack a 3′ polyA tail, and thus the initial step of poly-A selection on oligo-dT beads was omitted. Libraries were sequenced on a desktop sequencerⁱ with a 600-cycle kit (v3) to generate 2 × 300-bp paired-end reads.

Index adapters were removed from raw sequence reads using software,^j,24 and trimmed reads were screened against a database,^f,k to remove vector sequences. Given that virus strains were sequenced directly from the plasma of infected animals, host nucleic acids accounted for ≥98% of the reads. Therefore, assembly of viral genomic sequences was accomplished using template-assisted assembly where reference genomes (strain Singer, accession L32875; New York’93, accession AF502399; or Osloss, accession M96687) were used to discriminate viral from nonviral sequences in the dataset. Reads that mapped to the reference viral genome were then assembled de novo for each plasma sample.^f

Phylogenetic and recombination analyses of full-length BVDV genomes

To compare isolates from this study with others from around the world, a neighbor-joining tree based on HKY distance was created from a MUSCLE alignment of the coding RNA sequence (CDS) of 19 full-length genomes generated in this study along with 75 strains available in GenBank.^f The best nucleotide substitution model for maximum likelihood analyses was determined as described above.^g The appropriate model was then used for phylogenetic inference according to the maximum likelihood algorithm using PhyML methods.^f,15 The robustness of each tree was tested using 1,000 bootstrap analyses. GenBank accession numbers for strains used in these analyses are listed in Supplemental Table 1 (available online at http://vdi.sagepub.com/content/by/supplemental-data).

The MUSCLE alignment was also analyzed using a software package^18,l to construct a neighbor-net phylogenetic network. The phi test for recombination⁶ was used to evaluate the sequences for recombination.^l Recombination was further evaluated using recombination detection software^m using the default settings for a full exploratory recombination scan.^{22,23,27,30,35}

Results

Genetic typing of BVDV

Of the 131 calves that tested positive for BVDV infection by antigen ELISA, 119 (91%) were verified as positive for persistent infection by IHC on skin (ear notch) specimens and by RT-PCR on plasma. The 12 animals not testing positive were removed from the study.

To determine the prevalence of BVDV genotypes in this population and assign each strain to a subgroup, phylogenetic analyses were performed on a fragment of the 5′-UTR with 2 commonly used methods, maximum likelihood (ML; Fig. 1A) and neighbor-joining (NJ; Supplemental Fig. 1, available online at http://vdi.sagepub.com/content/by/supplemental-data). These methods resulted in the identical classification of BVDV strains. Of the 119 confirmed PI calves, 110 (92%) were infected with BVDV-1. Of those, 99 (83.2%) were infected with BVDV-1b and 11 (9.2%) were infected with BVDV-1a. The remaining 9 (7.6%) were infected with BVDV-2 (Fig. 1A). No strains were classified as BVDV-2b; however, 2 strains, “83 TX” and “76 MO,” had higher pairwise sequence identity with BVDV-2c strains from Germany (98–99%) than prototypic U.S. BVDV-2a strains 890 (92%) and NewYork’93 (94%).

Figure 1.

Maximum likelihood (ML) analyses of Bovine viral diarrhea virus (BVDV) 5′-UTR (untranslated region) sequences. A. Unrooted ML phylogeny based on the 5′-UTR of 119 BVDV strains from this study as well as reference sequences from GenBank. ML phylogeny construction was carried out using the Kimura model of nucleotide substitution with a gamma shape parameter of 0.30. The log likelihood of the inferred phylogeny is −3,223. Classification of subtypes was based on reference sequences.^8,19,41 The bootstrap values from the consensus tree have been added to the topology of the identical clades on the ML tree. The scale bar represents substitutions per site. Colors indicate BVDV subtypes. B. ML phylogeny based on the 5′-UTR of 9 BVDV-2 strains from this study as well as reference sequences from around the world. ML phylogeny construction was carried out using the Kimura model of nucleotide substitution with a gamma shape parameter of 0.23. The log likelihood of the inferred phylogeny is −1,154. Samples from this study are in bold font and named by their sample collection order followed by the buyer’s state of origin, where UNK is unknown. Reference strains are listed by the isolate name followed by the country of isolation. Numbers indicate the percentage of 1,000 bootstrap replicates that support each group. The scale bar represents substitutions per site.

To better characterize BVDV-2 strains from our study, a basic local alignment search toolⁿ was used to identify the most similar sequences in GenBank, and these strains were included together with additional reference BVDV-2 strains from around the world for phylogenetic analyses. Both ML (Fig. 1B) and NJ (data not shown) trees were generated. Separation of subtype 2b from the other BVDV-2 strains was strongly supported in both NJ and ML tree–based inferences. In contrast, the separation of BVDV subtypes 2a and 2c was not statistically supported by either method as evident from the low bootstrap values and variable strain placement between trees. Phylogenetic analysis and percentage of similarity determined by pairwise distances confirmed that 2 of the BVDV-2 strains from this study were most similar to BVDV-2c strains (bootstrap value >75%); and further revealed that one BVDV-2 strain from this study was most similar to a distinct subset of BVDV-2 isolates containing U.S. CD87³⁶ and German 2a strains⁴ (bootstrap values >85%). However, given the overall low bootstrap values in the 5′-UTR NJ and ML trees, BVDV-2 strains could not be confidently subtyped.

BVDV prevalence by location

The prevalence of each genotype was also determined by the order buyer’s state of origin. The 119 PI calves were purchased at auction by order buyers from at least 5 central states; however, the precise geographic origin of each animal was not known. When analyzed by buyer’s state of origin, no significant difference was observed in the prevalence of BVDV subtypes (Table 1). Furthermore, there was no evidence to suggest geographic-related clustering of subtypes (Supplemental Fig. 1).

Table 1.

Prevalence of Bovine viral diarrhea virus (BVDV) subtypes in persistently infected cattle by order buyer location.*

State of buyer	No. of animals	BVDV subtype
State of buyer	No. of animals	1a	1b	2
Alabama	3	0	3 (100)	0
Mississippi	2	1 (50)	1 (50)	0
Missouri	35	2 (5.7)	30 (85.7)	3 (8.6)
Oklahoma	33	3 (9.1)	28 (84.8)	2 (6)
Texas	22	4 (18)	17 (77.3)	1 (5)
Unknown	24	1 (4.2)	20 (83.3)	3 (12.5)
Total	119	11 (9.2)	99 (83.2)	9 (7.6)

Subtype determined by phylogenetic analyses of the 5’-UTR. Numbers in parentheses are percentages.

Generation of 19 full-length BVDV genomes

To investigate the genetic variation among strains in our study and to increase the resolution of BVDV-2 subtyping, all 9 BVDV-2 strains, 2 BVDV-1a strains, and 8 BVDV-1b strains were selected from across the phylogenetic tree for complete genome sequencing (Table 2). Each of the 19 genomes assembled de novo into a single contig with an average genome length of 12,245 nucleotides (minimum 12,192; maximum 12,294) with mean read coverages from 33- to 1,059-fold.

Table 2.

De novo whole genome assemblies of 19 Bovine viral diarrhea virus (BVDV) strains.*

Animal ID	Strain ID	BVDV subtype	Genome size (nt)	5′-UTR (nt)	CDS (nt)	3′-UTR (nt)	GenBank accession
AL 6	USMARC-51998	1b	12,227	373	1,1697	157	KP941581
MO 113	USMARC-55926	1b	12,212	367	1,1697	148	KP941592
MO 125	USMARC-53847	1b	12,225	372	1,1694	159	KP941583
MO 18	USMARC-55478	1b	12,192	372	1,1697	123	KP941587
MO 80	USMARC-55925	1b	12,224	371	1,1697	156	KP941591
OK 104	USMARC-55924	1b	12,222	375	1,1697	150	KP941590
TX 55	USMARC-55923	1b	12,250	384	1,1697	169	KP941589
UNK 21	USMARC-55922	1b	12,243	372	1,1697	174	KP941588
OK 88	USMARC-55477	1a	12,270	370	1,1697	203	KP941586
OK 77	USMARC-53875	1a	12,294	372	1,1697	225	KP941584
MO 76	USMARC-53873	2	12,273	373	1,1694	206	KP941582
OK 92	USMARC-55476	2	12,267	371	1,1694	202	KP941585
TX 83	USMARC-60768	2	12,236	373	1,1694	169	KT832821
MO 61	USMARC-60767	2	12,267	371	1,1694	202	KT832820
MO 79	USMARC-60766	2	12,271	373	1,1694	204	KT832819
UNK 106	USMARC-60765	2	12,245	378	1,1694	173	KT832818
UNK 37	USMARC-60764	2	12,235	372	1,1694	169	KT832817
UNK 15	USMARC-60779	2	12,239	372	1,1694	173	KT832822
OK 65	USMARC-60780	2	12,272	375	1,1694	203	KT832823

CDS = complete coding RNA sequence; nt = nucleotides; UNK = unknown; UTR = untranslated region.

The BVDV-1b strain USMARC-53847 from animal “125 MO” had a 3-nucleotide in-frame deletion in the E2 surface glycoprotein gene when compared to strains from this study as well as reference strains. No other deletions, insertions, gene duplications, or recombination events were detected in the coding sequences of the 19 genomes; however, 0–50 nucleotides were missing from the genome termini compared to reference strains.

Phylogenetic and recombination analyses of BVDV genomes

The CDSs of these 19 genomes were also compared to those of 75 full-length BVDV genomes isolated from around the world and available in GenBank. Phylogenetic analyses of the CDSs resulted in well-supported consensus trees, although there were some differences in the phylogenetic relationships inferred by the ML (Fig. 2) and NJ (Supplemental Fig. 2, available online at http://vdi.sagepub.com/content/by/supplemental-data) tree–based inferences. For BVDV-1, 100% concordance was observed between strain subtyping based on phylogenetic analyses of sequences from the 5′-UTR (Fig. 1) or the CDS (Fig. 2) for the 10 isolates from this study. Furthermore, the NJ and ML tree topology for BVDV-1 isolates was similar, although differences in the inferred relationships between the non–BVDV-1a or -1b isolates were observed.

Figure 2.

Mid-point rooted maximum likelihood (ML) tree of Bovine viral diarrhea virus (BVDV) complete coding sequences. A mid-point rooted ML tree was constructed on the basis of the complete coding RNA sequence of 19 strains from this study along with 75 full-length viral genomes available in GenBank (listed by their strain ID followed by country of origin). ML phylogeny construction was carried out using the general time reversible model of nucleotide substitution with a gamma shape parameter of 1.70 and a proportion of invariant sites equal to 0.41. The log likelihood of the inferred phylogeny is −190,930. Numbers indicate the percentage of 1,000 bootstrap replicates that support each group. The scale bar represents substitutions per site. Colors indicate BVDV subtypes.

For BVDV-2 stains, phylogenetic analysis of the complete CDS provided more conclusive genetic classification and revealed several distinct genetic clusters supported by high bootstrap values (Fig. 2). The first group contained the prototypic U.S. BVDV-2a strain 890⁵ and 4 strains from our study (USMARC-60766, -60779, -60767, and -55476). The second group contained prototypic German BVDV-2c strains¹⁹ and 2 strains from our study (USMARC-53873 and -60768). The third group contained BVDV-2b strains and no strains from our study. There were also several U.S. and Chinese strains as well as 3 strains from our study (USMARC-60764, -60765, and -60780) that were distinct from these previously defined clades (Fig. 2). JX-04/China, USMARC-60764, PI99/USA, and PI103/USA formed a group sister to the established BVDV-2c clade. This relationship was strongly supported in both NJ and ML trees. Similarly, USMARC-60780 and the recombinant strain JZ05-1/China⁴² formed a group sister to the established BVDV-2a clade. Positioning of this group was identical in both NJ and ML trees. In addition, there were 3 single lineages (SH-28/China, USMARC-60765, and HLJ-10/China) that were distinct from the well-supported BVDV-2 clades. These strains had low bootstrap support and variable placement between NJ and ML trees. Therefore, the relationships among these strains could not be confidently resolved.

The low bootstrap values and variable placement of certain BVDV-2 strains between NJ and ML trees could be indicative of viral recombination.^20,42 Therefore, sequences were tested for recombination using several algorithms within a recombination detection software package.^m Strains 3156/China, JZ05-1/China, and ILLNC/USA were identified as recombinants as previously reported^20,42; however some differences in parental strains were identified.⁴² In particular, JZ05-1/China was previously found to have BVDV-2a 11F011/South Korea as the major parent and an unknown second strain as the minor parent. In contrast, this dataset found the major parent to be USMARC-60780. Neighbor-net phylogenetic networks were constructed to more accurately depict the phylogenetic relationships of viral recombinants (Fig. 3). No strains from our study tested positive for recombination.

Figure 3.

Neighbor-joining (NJ) phylogenetic network of Bovine viral diarrhea virus (BVDV) complete coding RNA sequences (CDSs). A. NJ phylogenetic network was constructed in Splits Tree on the basis of the CDSs of 19 individuals from this study along with 75 full-length viral genomes available in GenBank (listed by their strain ID followed by country of origin). B. NJ phylogenetic network was constructed on the basis of the CDSs of 9 BVDV-2 strains from this study and 40 full-length BVDV-2 genomes available in GenBank. Numbers indicate the percentage of 1,000 bootstrap replicates that support each node. Colors indicate BVDV subtypes.

The percentage of similarity of pairwise distances within and between BVDV groups was calculated. The nucleotide sequence divergence among virus isolates of the same group was 6% for BVDV-2a, 8% for BVDV-2c, and 11% for BVDV-2b. The ranges of sequence divergence between subgroups were 10% between BVDV-2a and -2c, and 16% between BVDV-2b and the other BVDV-2 clades. This is in contrast to BVDV-1 subtypes, which have a nucleotide sequence divergence of ≤16% among virus isolates of the same group and 17–23% divergence between subtypes (data not shown).

Individual gene analyses of BVDV-2 stains

To find a genomic region that provides a more stable, well-supported phylogeny than the 5′-UTR for the subtyping of BVDV-2 strains, phylogenetic relationships of individual genes were reconstructed by NJ and ML approaches. Similar tree topology and support were obtained using the 2 methods, and only the ML trees are shown (Fig. 4). Parameters used in ML phylogeny reconstruction are listed in Table 3. Several of the coding genes had strongly supported separation of BVDV-2a, -2b, and -2c clades with bootstrap values >95% (nodes marked by an asterisk [*] in Fig. 4). The strains that fall outside these clades were more weakly supported and had variable positions between individual gene trees, although they always fell outside the previously defined clades. Coding regions for NPro, E2, NS5A, and NS5B had the highest resolution and statistical support, with tree topologies most similar to the CDS (Fig. 4). NS2/3 also had high resolution and statistical support; however, recombinant strain JZ05-1/China was placed in the BVDV-2b clade for this genomic region.

Figure 4.

Individual gene phylogeny inferred by maximum likelihood (ML) analyses. Bovine viral diarrhea virus 2 (BVDV-2) phylogeny inferred by ML analyses of the complete coding RNA sequence (CDS; A) or individual genes (B). Support for tree nodes was determined by 1,000 pseudo-alignments. A. The tree based on the CDS is shown to display the names and relationships of BVDV-2 strains. Branches are color coded for each subgroup where 2a is green, 2b is blue, and 2c is red. The strains that fall outside these clades are black. B. ML consensus trees collapsed to nodes with at least 70% support. The asterisk (*) at internal nodes represents strong statistical support for the clade with bootstrap values >95%. Listed in each panel are the number of nucleotides analyzed and the percent of identical sites in each alignment.

Table 3.

Maximum likelihood best-fit models used in Figure 4.*

Genomic region	Substitution model	Gamma shape parameter	Proportion invariant sites	BIC	Log likelihood
5′-UTR	K2 +G + I	0.53	0.69	3,564	−1,763
Npro	K2 + G	0.36	NA	6,810	−2,909
Capsid	TN93 + I	NA	0.54	4,219	−1,615
E^RNS	TN93 + G	0.38	NA	8,522	−3,759
E1	K2 + G + I	1.57	0.44	7,736	−3,419
E2	TN93 + G	0.48	NA	15,388	−7,204
NS2/3	TN93 + G + I	1.91	0.5	38,498	−19,328
NS4A	TN93 + G	0.25	NA	21,439	−1,027
NS4B	TN93 + G	0.33	NA	12,257	−5,567
NS5A	TN93 + G	0.47	NA	18,387	−8,648
NS5B	GTR + G	0.38	NA	24,442	−11,670
CDS	GTR + G + I	1.91	0.46	135,437	−68,847

BIC = Bayesian information criterion; CDS = complete coding RNA sequence; GTR = general time reversible; G = gamma distribution; I = invariant sites; K2 = Kimura 2-parameter; NA = not applicable; TN93 = Tamura–Nei; UTR = untranslated region. Models with the lowest BIC scores were considered to describe the substitution pattern the best, and were thus used for maximum likelihood phylogeny reconstruction.

Discussion

The continuous and systematic characterization of BVDV strains circulating in the United States is needed by scientists, veterinarians, and ranchers to evaluate both prevention and control strategies, and to make informed management decisions. Therefore, the purpose of our study was to estimate the current prevalence and diversity of BVDV subtypes by sampling PI beef calves in the central United States Phylogenetic analyses of 119 viral strains from at least 5 central states revealed that BVDV-1b is the most common subtype infecting PI cattle in the region. These results are consistent with the high prevalence of BVDV-1b reported in previous studies.^13,14,33,36 Interestingly, the genotype profile reported herein is comparable to that reported 10 years earlier at the same U.S. feedyard.¹¹ This suggests that the prevalence of subtypes in PI beef calves in the region has remained similar in the last decade despite potential changes in the source of animals, cattle management practices, vaccination protocols, or implementation of BVDV control programs during this time. Information regarding the use of vaccines in the herd of origin (breeding herd) for each animal was unavailable; thus, the impact of vaccination on subtype prevalence in this study is unknown.

Accurate genetic typing of viruses is important for the precise classification of new viral isolates, to understand what genetic diversity exists in a population, and to track virus evolution over time; however, uncertainty exists in BVDV classification caused by a lack of specific criteria for subtyping new virus isolates.²⁹ Previous phylogenetic studies on BVDV have been based on partial genomic sequences of the 5′-UTR or coding regions such as NPro, E2, or NS2/3. Of these, 5′-UTR is most commonly used and is ~380 nucleotides in length and contains both highly conserved regions as well as short variable regions.^16,31 The conserved regions are important to allow detection of all pestivirus species; however, these regions contribute to the low bootstrap values observed in 5′-UTR phylogenetic trees.^3,44 Despite this, it has been demonstrated that classification of BVDV-1 subtypes is consistent regardless of the genomic region used for analysis,^3,37,39–41 with few exceptions resulting from homologous recombination.^20,26,42 In contrast to what has been observed for BVDV-1 subtyping, our study revealed that BVDV-2 strains could not be confidently subtyped based on phylogenetic analyses of the 5′-UTR alone.

To better define the variability of U.S. BVDV isolates and further investigate the division of BVDV-2 isolates into subtypes, complete genome sequences were obtained by sequencing viral RNA from the plasma of infected calves. These full-length viral sequences were then compared with viral genomes from around the world. BVDV-2 subtypes were found to have much lower variability between subtypes (10–16%) than between BVDV-1 subtypes (20–24%). Nevertheless, phylogenetic analyses revealed 100% confidence in the subdivision of subtypes 2a, 2b, and 2c by NJ and ML methods. In addition, 3 strains from our study as well as strains previously classified as 2a were found to be distinct from these previously defined clades (Figs. 2 –4). These strains had unstable placement in NJ and ML trees and could represent genetic intermediates; however, they tested negative for recombination, and thus are not genetic mosaics of the other subgroups. Additional full-length BVDV-2 sequences may help clarify the relationships between these strains.

Although analysis using the complete genome provides the highest resolution of the relationships among subtypes, it is not practical for large surveys. Therefore, single-gene analyses were carried out to determine which genomic segment most clearly resolves the relationship between subtypes and most accurately reflects the tree-based inferences using complete genome sequencing (Fig. 4). Although the relationships between the genetic intermediate strains were not fully resolved in the individual gene trees, many coding regions strongly supported the division of BVDV-2a, -2b, and -2c (Fig. 4). In addition, several regions further supported the division of the intermediate strains. However, recombinant strain JZ05-1/China placed in different clades depending on the genomic region analyzed, highlighting the problems with single-gene analysis for typing viruses capable of recombination.

In summary, the present study shows that the 5′-UTR can be used to segregate BVDV into major genotypes, but may be less suited for the segregation of BVDV-2 strains into subgroups. Coding regions provided more accurate subtyping of BVDV-2 strains; however, recombination events can lead to misclassification of strains typed by individual gene analyses. The clinical significance of BVDV-2 subtype division is currently unknown; further work is needed to determine the serologic relatedness of these subtypes. Regardless, the enhanced resolution of BVDV-2 subtypes provided by these analyses is critical for studying virus evolution and provides information for more accurate monitoring. This information is important to control programs and may influence the choice of BVDV strains to be used in vaccines.

Footnotes

Acknowledgements

We thank Brad Sharp, Sue Hauver, Bob Lee, Stacy Bierman, and Jacky Carnahan for technical support and Jan Watts for secretarial support. We also thank Steve Johnson for his assistance in the acquisition of samples for this study.

Authors’ contributions

AM Workman drafted the manuscript. All authors contributed to conception and design of the study; contributed to acquisition, analysis, and interpretation of data; critically revised the manuscript; gave final approval; and agreed to be accountable for all aspects of the work in ensuring that questions relating to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Authors’ note

The use of product and company names is necessary to accurately report the methods and results; however, the United States Department of Agriculture (USDA) neither guarantees nor warrants the standard of the products. The use of names by the USDA implies no approval of the product to the exclusion of others that may also be suitable. The USDA is an equal opportunity provider and employer.

a.

TRIzol, Life Technologies, Grand Island, NY.

b.

OneStep RT-PCR kit, Qiagen Inc., Valencia, CA.

c.

Exo I, Life Technologies, Grand Island, NY.

d.

BigDye, PE Applied Biosystems, Foster City, CA.

e.

ABI 3730, PE Applied Biosystems, Foster City, CA.

f.

Geneious (version 8.0.5²¹), Biomatters Inc., Newark, NJ.

g.

Molecular Evolutionary Genetics Analysis version 6 (MEGA6 version 6.06-mac³⁸), Arizona State University, Tempe, AZ.

h.

TruSeq RNA sample preparation kit, Illumina Inc., San Diego, CA.

i.

MiSeq, Illumina Inc., San Diego, CA.

j.

Cutadapt (version 1.9.1²²), Python Software Foundation, Beaverton, OR.

k.

UniVec_Core database, NCBI, Bethesda, MD.

l.

Splits Tree (version 4.13.1), University of Tübingen, Tübingen, Germany.

m.

Recombination detection program version 4, University of Cape Town, Cape Town, South Africa.

n.

BLAST, NCBI, Bethesda, MD.

Declaration of conflicting interests

The author(s) declared the following potential conflicts of interest with respect to the research, authorship, and/or publication of this article: David Sjeklocha is employed by Cattle Empire LLC. Cattle Empire LLC did not have a role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Funding for this research was provided by the United States Department of Agriculture, Agricultural Research Service (CRIS 3040-32000-031-00D).

References

Alpay

Yeşilbağ

Serological relationships among subgroups in bovine viral diarrhea virus genotype 1 (BVDV-1). Vet Microbiol 2015;175:1–6.

Baker

JC.

The clinical manifestations of bovine viral diarrhea infection. Vet Clin North Am Food Anim Pract 1995;11:425–445.

Becher

. Phylogenetic analysis of pestiviruses from domestic and wild ruminants. J Gen Virol 1997;78(Pt 6):1357–1366.

Beer

. Phylogenetic analysis of the 5’-untranslated region of German BVDV type II isolates. J Vet Med B Infect Dis Vet Public Health 2002;49:43–47.

Bolin

Ridpath

JF.

Differences in virulence between two noncytopathic bovine viral diarrhea viruses in calves. Am J Vet Res 1992;53:2157–2163.

Bruen

. A simple and robust statistical test for detecting the presence of recombination. Genetics 2006;172:2665–2681.

Collett

. Comparisons of the pestivirus bovine viral diarrhoea virus with members of the Flaviviridae. J Gen Virol 1988;69(Pt 10):2637–2643.

Flores

. Phylogenetic analysis of Brazilian bovine viral diarrhea virus type 2 (BVDV-2) isolates: evidence for a subgenotype within BVDV-2. Virus Res 2002;87:51–60.

Fulton

RW.

Host response to bovine viral diarrhea virus and interactions with infectious agents in the feedlot and breeding herd. Biologicals 2013;41:31–38.

10.

Fulton

RW.

Impact of species and subgenotypes of bovine viral diarrhea virus on control by vaccination. Anim Health Res Rev 2015;16:40–54.

11.

Fulton

. Evaluation of diagnostic tests used for detection of bovine viral diarrhea virus and prevalence of subtypes 1a, 1b, and 2a in persistently infected cattle entering a feedlot. J Am Vet Med Assoc 2006;228:578–584.

12.

Fulton

. Multiple diagnostic tests to identify cattle with Bovine viral diarrhea virus and duration of positive test results in persistently infected cattle. Can J Vet Res 2009;73:117–124.

13.

Fulton

. Bovine viral diarrhoea virus (BVDV) subgenotypes in diagnostic laboratory accessions: distribution of BVDV1a, 1b, and 2a subgenotypes. Vet Microbiol 2005;111:35–40.

14.

Fulton

. Prevalence of bovine viral diarrhea virus (BVDV) in persistently infected cattle and BVDV subtypes in affected cattle in beef herds in south central United States. Can J Vet Res 2009;73:283–291.

15.

Guindon

. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol 2010;59:307–321.

16.

Harasawa

Giangaspero

A novel method for pestivirus genotyping based on palindromic nucleotide substitutions in the 5′-untranslated region. J Virol Methods 1998;70:225–230.

17.

Hessman

. Evaluation of economic effects and the health and performance of the general cattle population after exposure to cattle persistently infected with bovine viral diarrhea virus in a starter feedlot. Am J Vet Res 2009;70:73–85.

18.

Huson

Bryant

Application of phylogenetic networks in evolutionary studies. Mol Biol Evol 2006;23:254–267.

19.

Jenckel

. Mixed triple: allied viruses in unique recent isolates of highly virulent type 2 bovine viral diarrhea virus detected by deep sequencing. J Virol 2014;88:6983–6992.

20.

Jones

Weber

EL.

Homologous recombination in bovine pestiviruses. Phylogenetic and statistic evidence. Infect Genet Evol 2004;4:335–343.

21.

Kearse

. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 2012;28:1647–1649.

22.

Martin

Rybicki

RDP: detection of recombination amongst aligned sequences. Bioinformatics 2000;16:562–563.

23.

Martin

. A modified bootscan algorithm for automated identification of recombinant sequences and recombination breakpoints. AIDS Res Hum Retroviruses 2005;21:98–102.

24.

Martin

Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnetjournal 2011;17:10–12.

25.

McClurkin

. Production of cattle immunotolerant to bovine viral diarrhea virus. Can J Comp Med 1984;48:156–161.

26.

Nagai

. Phylogenetic analysis of bovine viral diarrhea viruses using five different genetic regions. Virus Res 2004;99:103–113.

27.

Padidam

. Possible emergence of new geminiviruses by frequent recombination. Virology 1999;265:218–225.

28.

Pellerin

. Identification of a new group of bovine viral diarrhea virus strains associated with severe outbreaks and high mortalities. Virology 1994;203:260–268.

29.

Peterhans

. Cytopathic bovine viral diarrhea viruses (BVDV): emerging pestiviruses doomed to extinction. Vet Res 2010;41:44.

30.

Posada

Crandall

KA.

Evaluation of methods for detecting recombination from DNA sequences: computer simulations. Proc Natl Acad Sci U S A 2001;98:13757–13762.

31.

Ridpath

Bolin

SR.

Differentiation of types 1a, 1b and 2 bovine viral diarrhoea virus (BVDV) by PCR. Mol Cell Probes 1998;12:101–106.

32.

Ridpath

. Segregation of bovine viral diarrhea virus into genotypes. Virology 1994;205:66–74.

33.

Ridpath

. Prevalence and antigenic differences observed between Bovine viral diarrhea virus subgenotypes isolated from cattle in Australia and feedlots in the southwestern United States. J Vet Diagn Invest 2010;22:184–191.

34.

Ridpath

. Change in predominance of Bovine viral diarrhea virus subgenotypes among samples submitted to a diagnostic laboratory over a 20-year time span. J Vet Diagn Invest 2011;23:185–193.

35.

Smith

JM.

Analyzing the mosaic structure of genes. J Mol Evol 1992;34:126–129.

36.

Tajima

Dubovi

EJ.

Genetic and clinical analyses of bovine viral diarrhea virus isolates from dairy operations in the United States of America. J Vet Diagn Invest 2005;17:10–15.

37.

Tajima

. Prevalence of genotypes 1 and 2 of bovine viral diarrhea virus in Lower Saxony, Germany. Virus Res 2001;76:31–42.

38.

Tamura

. MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Mol Biol Evol 2013;30:2725–2729.

39.

Vilcek

Discrepancies in genetic typing of bovine viral diarrhoea virus isolates. Vet Microbiol 2005;106:153–154.

40.

Vilcek

. Genetic diversity of BVDV: consequences for classification and molecular epidemiology. Prev Vet Med 2005;72:31–35; discussion 215–219.

41.

Vilcek

. Bovine viral diarrhoea virus genotype 1 can be separated into at least eleven genetic groups. Arch Virol 2001;146:99–115.

42.

Weber

. Homologous recombination in pestiviruses: identification of three putative novel events between different subtypes/genogroups. Infect Genet Evol 2015;30:219–224.

43.

Workman

. Full-length coding sequences for 12 bovine viral diarrhea virus isolates from persistently infected cattle in a feedyard in Kansas. Genome Announc 2015;3:pii:e00487–15.

44.

Xia

. Molecular phylogenetic analysis of bovine viral diarrhoea virus: a Bayesian approach. Virus Res 2007;130:53–62.