- Email: [email protected]
Evolution of Bovine viral diarrhea virus in Canada from 1997 to 2013
Virology 509 (2017) 232–238
Contents lists available at ScienceDirect
Virology journal homepage: www.elsevier.com/locate/yviro
Evolution of Bovine viral diarrhea virus in Canada from 1997 to 2013 ⁎
Adam Chernick , Frank van der Meer
MARK
⁎
Ecosystem and Public Health, Faculty of Veterinary Medicine, University of Calgary, 3330 Hospital Dr NW, Calgary, AB T2N 1N4, Canada
A R T I C L E I N F O
A BS T RAC T
Keywords: Bovine viral diarrhea virus Evolution Next-generation sequencing ssRNA viruses Phylogenetics Canada
Bovine viral diarrhea virus (BVDV) is a rapidly evolving, single-stranded RNA virus and a production limiting pathogen of cattle worldwide. 79 viral isolates collected between 1997 and 2013 in Canada were subjected to next-generation sequencing. Bayesian phylogenetics was used to assess the evolution of this virus. A mean substitution rate of 1.4×10−3 substitutions/site/year was found across both BVDV1 and BVDV2. Evolutionary rates in the E2 gene were slightly faster than other regions. We also identified population structures below the sub-genotype level that likely have phenotypic implications. Two distinct clusters within BVDV2a are present and can be differentiated, in part, by a tyrosine to isoleucine mutation at position 963 in the E2 protein, a position implicated in the antigenicity of BVDV1 isolates. Distinct clustering within all sub-genotypes, particularly BVDV2a, is apparent and could lead to new levels of genotypic classification. Continuous monitoring of emerging variants is therefore necessary.
1. Introduction Bovine viral diarrhea virus (BVDV) is a single-stranded RNA (ssRNA) virus that is a production limiting pathogen of cattle with a global distribution (Ridpath, 2010). It can be categorized in two genotypes or species, BVDV1 and BVDV2. Numerous sub-genotypes of BVDV1 have been reported as well as two proposed BVDV2 subgenotypes (Vilček et al., 2001; Flores et al., 2002). In Canada and the United States BVDV1a, BVDV1b and BVDV2a are currently in circulation. Previous studies have described notable genetic diversity within and between each of these sub-genotypes (Vilcek et al., 2004) as well as within BVDV1a in western Canada (Chernick et al., 2014). In the latter study, the substitution rate of BVDV1a isolates was estimated to be 1.26×10−3 substitutions/site/year in the E1-E2 gene across all codon positions. This rapid evolution drives many important phenotypic changes in BVDV including extensive antigenic variability (Ridpath et al., 2000). High substitution rates are not unique to BVDV and are often noted in other ssRNA viruses including type 1a Hepatitis C virus (Yuan et al., 2013), HIV (Suzuki et al., 2000) and others (Duffy et al., 2008). This rapid evolution stems largely from the inherent properties of these viruses and their error-prone replication. In the case of BVDV specifically, the unique immune environment associated with persistent infections also appears to play an important role in generating novel viral variants (Dow et al., 2015).
BVDV is a Pestivirus within the family Flaviviridae. It shares a genome structure similar to these related viruses with 5′ and 3′ untranslated regions (UTRs), structural genes towards the 5′ end and most nonstructural genes to the 3′ end. The UTRs are both highly conserved with the 5′UTR acting as a ribosomal entry site (Fletcher and Jackson, 2002) and 3′UTR also essential to replication (Yu et al., 1999). Following the 5′UTR is the NPro, a viral protease responsible for cleaving some viral proteins from one another (Wiskerchen et al., 1991). The structural proteins follow and include the envelope proteins (Erns, E1 and E2) and the capsid. Erns is secreted from infected cells (Rumenapf et al., 1993) and binds to host cell surfaces (Iqbal et al., 2000) while E1 and E2 are both surface glycoproteins important to host cell binding and entry (Ronecker et al., 2008) as well as being highly antigenic (Paton et al., 1992; Deregt et al., 1998a, 1998b). The p7 gene sits between the structural and non-structural genes and is required for infectivity (Harada et al., 2000). NS2 and NS3 are non-structural genes encoding for a protein that performs multiple functions including a serine protease (Tautz et al., 1997) and helicase (Warrener and Collett, 1995). NS4A is a co-factor for the NS2/3 serine protease activity (Xu et al., 1997) and NS4B appears to play a role in viral genome replication and host membrane remodelling (Weiskircher et al., 2009). NS5A (Isken et al., 2014) and NS5B (Zhong et al., 1998) are both components of the RNA-dependent RNA polymerase complex. Several phylogenetic and genetic diversity studies of pestiviruses
Abbreviations: BVDV, Bovine viral diarrhea virus; E, envelope; ESS, effective sample size; Gbp, gigabase pairs; GTR, general time-reversible; HIV, human immunodeficiency virus; HPD95, 95% highest probability density; MRCA, most recent common ancestor; NS, non-structural; RB, reversible-jump based; ssRNA, single-stranded RNA; SYM, symmetrical; TIM, transition; UK, United Kingdom; UTR, untranslated region; VSACC, Veterinary Sciences Animal Care Committee ⁎ Corresponding authors. E-mail addresses: [email protected] (A. Chernick), [email protected] (F. van der Meer). http://dx.doi.org/10.1016/j.virol.2017.06.024 Received 8 April 2017; Received in revised form 20 June 2017; Accepted 21 June 2017 0042-6822/ © 2017 Elsevier Inc. All rights reserved.
Virology 509 (2017) 232–238
A. Chernick, F. van der Meer
and BVDV have also been conducted, although none exhaustively examine Canadian field isolates. With respect to the age of the entire BVDV lineage, all BVDV1 and BVDV2 isolates are believed the share a most recent common ancestor about the year 1743 (95% highest probability density (HPD95) 1373–1926) (Liu et al., 2009). Other studies have examined the spatial movement of BVDV while inferring evolutionary rates about 9.3×10−3 substitutions/site/year in the 5′ untranslated region (5′UTR) among Italian isolates (Luzzago et al., 2012). Phylogenetic and genetic diversity studies of European isolates are available in large part due to the need to inform BVDV-eradication and control programs in those countries. For example, isolates in Denmark span BVDV1 sub-genotypes a, b, d, e, f, g and h (Uttenthal et al., 2005) while a collection from France included sub-genotypes b, d, e and l (Jackova et al., 2008). Prior to these phylogenetic studies numerous experiments were conducted using anti-BVDV antibodies to assess strain variation. Although many antibodies were found to be broadly neutralizing against genetically distinct isolates (Deregt et al., 1990), significant differences that correlated with genetic distance were also noted (Deregt et al., 1992). Panels of antibodies capable of distinguishing BVDV1 from BVDV2 as well as some sub-genotypic variation have also been generated (Deregt et al., 1998a, 1998b). Although there is merit in classifying isolates based on antibody binding assays, it has been noted that such a practice is “suggestive at best”, implying that sequence data is a more reliable tool even if it does not produce such an explicit connection to the phenotypic/antigenic properties of a given isolate (Dubovi, 1992). Sequence data also allows for more comprehensive evolutionary studies that describe the origins and evolutionary rates of BVDV. In this study we describe the evolution of BVDV in Canada using a collection of viral isolates archived from 1997 to 2013 and obtained from academic and diagnostic laboratories across the country as well as from active field sampling efforts. We employed a next generation sequencing (NGS) approach in conjunction with Bayesian phylogenetics to generate an integrated model of BVDV1 and BVDV2 evolution within the Canadian cattle population. This both describes the evolutionary history of this pathogen in Canada as well as important variability among currently circulating BVDV isolates. Apart from the practical use of this information in the development and maintenance of virus control programs including vaccination, these results also contribute to the understanding of the evolutionary capacity of ssRNA viruses which include a number of important human and animal pathogens.
Table 1 Primer sequences and references for the four amplicon method used to generate DNA for next-generation sequencing. Name
Sequence
Reference
Position
Frag1_fwd
CATGCCCATAGTAGGAC
107
Frag1_rev
GGGCAWACCATYTGGAAGGCYGG
Frag2_fwd
TGATAACAGGGGTACAAGGGC
Frag2_rev Frag3_fwd
RTTGYCCCTAGCGGTATAT GCAGATTTTGAAGAAAGACACTA
Frag3_rev
AAGAAGCCATCATCMCCACA
Frag4_fwd
AAGATCCACCCTTATGARGC
Frag4_rev
CTGTGTGCATTRARTGTAGTGTT
(Harada et al., 2000) (Tautz et al., 1997) (Warrener and Collett, 1995) N/A (Warrener and Collett, 1995) (Xu et al., 1997) (Xu et al., 1997) N/A
2853 2442 5314 4934 11547 10385 12511
Primer sequences and references for the four amplicon method used to generate DNA for next-generation sequencing. All positions are relative to the NADL reference genome (GenBank Accession M31182).
subjected to library preparation using the Nextera XT DNA Library Preparation kit (Illumina, San Diego, USA). All libraries were prepared with dual indices and sequenced with paired-end, 250 bp reads on a MiSeq using a 500-cycle v2 cartridge (Illumina, San Diego, USA). 2.2. Data Analysis Raw reads were imported into Geneious R9 (Kearse et al., 2012) for pre-processing and assembly. Reads were paired and assembled using the Geneious assembler set to high sensitivity and iterated up to 5 times. Consensus sequences were extracted from assemblies for all individual genes except p7 (leaving a total of 10 genes). These are available in GenBank via accession numbers KX169935 through KX170703 along with the relevant protein translations (Table S2). Gene-specific consensus sequence alignments were built using MUSCLE (Edgar, 2004). Alignments were visually inspected and corrected for errors. The HyPhy software package (Kosakovsky Pond et al., 2005) running on the Datamonkey server (Delport et al., 2010) was used to test for recombination. Specifically, the single breakpoint recombination (SBP) analysis (Kosakovsky Pond et al., 2006) was run using the HKY85 substitution model on each gene-specific alignment. Sequence alignments were imported into BEAUti to prepare the run parameter files for BEAST 2 (Bouckaert et al., 2014). Viral genes were analyzed independently to explicitly accommodate any differences in their phylogenies. All sequences were dated with their year of collection and separated into three partitions for the codon positions. Year was used as it was the most precise date data available for all isolates. A reversible-jump based (RB) substitution model accounting for ambiguities in the consensus sequences was employed to simultaneously explore different substitution models and infer the phylogeny using the mixture that best fit the data (Bouckaert et al., 2014, 2013). Six substitution models were evaluated including the Felsenstein (F81) (Felsenstein, 1981), HKY (Hasegawa et al., 1985), Tamura and Nei (TN93) (Tamura and Nei, 1993), transition (TIM), symmetrical (SYM) (Zharkikh, 1994) and generalized time-reversible (GTR) (Lanave et al., 1984) models. Separate relaxed exponential molecular clocks were used for each codon partition. Runs were allowed to continue until adequate convergence (ESS > 200 for most parameters) was reached. Phylogenetic trees were built using FigTree v1.4.2 (http://tree.bio.ed. ac.uk/software/figtree/) and Tracer was used to extract parameters of interest while verifying run quality. Two additional trees were built. The first using only the BVDV2a isolates in addition to BVDV2a and BVDV2b reference strains to further refine the clusters observed in the larger trees using partial 5′UTR sequences. The second includes all isolates in this study as well
2. Methods 2.1. Sample collection and processing Viral isolates were obtained from diagnostic and research laboratories across Canada (Table S1). This was the result of an extensive survey of laboratories working on BVDV and likely represents most historical isolates available for evaluation. All samples were collected in accordance with the animal care protocol VSACC AC11-0081/AC140144. All isolates had been passaged at least once at the laboratory of origin in Madin Darby Bovine Kidney or Bovine Turbinate cell culture prior to further characterization in our laboratory. RNA was extracted from an aliquot of each isolate using the E.Z.N.A. Viral RNA kit (Omega Bio-tek, Norcross, USA). A series of four reverse transcriptase PCRs were used to amplify overlapping fragments of the protein-coding region of the genome. See Table 1 for primer sequences and amplicon regions (Ridpath and Bolin, 1998; Couvreur et al., 2002; Greiser-Wilke et al., 1993; Gilbert et al., 1999). The amplicons were gel purified using the E.Z.N.A. Gel Extraction kit (Omega Bio-tek, Norcross, USA) and quantified using the Qubit 3.0 and the Qubit dsDNA HS Assay kit (Thermo Fisher Scientific, Waltham, USA). Amplicons were pooled in equimolar amounts, adjusted for amplicon length. Each pool was 233
Virology 509 (2017) 232–238
A. Chernick, F. van der Meer
as a selection of reference strains to place these isolates within the global context based on E2 sequences. Both were built as described above but without tip dates (5′UTR and E2) or codon partitioning (5′UTR only). To compare different gene trees two approaches were employed. First, to assess the overall similarity of the trees both visual inspections as well as a phylogenetic network analysis were employed. Although phylogenetic networks are not as easily interpreted as standard phylogenetic trees they do provide a means to visualize phylogenetic uncertainty or the incongruity between multiple trees. Second, to assess differences in branch lengths, Pearson correlations between gene trees were calculated using the branch length distances between all pairs of taxa in SPSS v22 (IBM, Armonk, USA). In addition to the correlation coefficients, scatterplots were drawn to assess the relationship between branch lengths for each pair of genes. To compliment the tree-based assessment, sequence alignments were also generated to highlight differences between certain clusters. Consensus sequences for these clusters of interest were generated from MUSCLE sequence alignments and compared to consensus sequences derived from other clusters and to older reference strain sequences from GenBank. These consensus sequences were allowed to contain ambiguities to account for any diversity that exists within subgenotypes.
Fig. 1. Mean clock rates with HPD95 intervals for each codon position in each gene phylogeny. Clock rates are measured in nucleotide substitutions per site per year. NS2 and NS4A have been excluded due to poor convergence during the Bayesian phylogenetic analysis.
3. Results The sequencing run generated a total of 4.49 Gbp across all samples and controls. The cluster density was 488 k/mm2 with 92.97% of all clusters passed quality filters. Based on re-sequencing of the PhiX control library a machine error rate of 2.04% was calculated. Across all clusters that passed the quality filter, an average of 74.94% of the resulting bases had a Q-score ≥ 30. Following assembly, it was found that full protein-coding region sequences were not obtained for all isolates due to inadequate coverage in certain regions, especially in the region amplified by the third RTPCR covering positions 4934–11547. This RT-PCR most commonly failed for BVDV2 isolates, indicating a possible genotypic bias in the primers. Excluding the p7 gene, consensus sequences for all 10 genes in the coding region were generated for 42 of 79 isolates. Nine out of 10 genes were sequenced successfully for another 10 isolates with the remaining 27 isolates yielding eight or fewer genes. See Table S1 for a break down of the genes sequenced for each isolate. Phylogenetic reconstruction was performed for each gene using all available isolate sequences. The trees for NPro, Erns, E1 and E2 contain the full set of 79 isolates. Regarding recombination, no signs of recombination were detected using the small sample AIC indicator in the SBP analysis except for in NS2 and NS5B. There was significant support for recombination in NS2 at nucleotide positions 374–375 in that alignment and NS5B at nucleotide positions 420–421 and 615 in that alignment. As the following phylogenetic analysis was performed on individual genes (with the exception of the phylogenetic network) no isolates were removed. The NS2 and NS5B phylogenies are for comparison only and cannot be used to infer the relatedness of these sequences. Bayesian evolutionary analysis was performed using BEAST 2 where each sequence alignment was assessed separately until satisfactory ESS ( > 200) was obtained for key parameters such as the posterior, tree likelihoods, tree height and clock rates. Notable exceptions were the NS2 and NS4A genes which demonstrated unconverged or bimodal distributions in several key parameters, respectively. The following statistics exclude NS2 and NS4A for this reason. The mean most recent common ancestor (MRCA) for all phylogenies was dated to 1498 with a mean HPD95 range of 1182–1750. All individual HPD95 ranges overlapped as well. For individual sub genotypes the mean MRCA (mean HPD95 range) across all genes were 1930 (1890–1966), 1955 (1915–1980) and 1950 (1902–1979) for BVDV1a, BVDV1b and
BVDV2a, respectively. The mean uncorrelated exponential clock rates across all sub-genotypes for each codon position were 7.57×10−4, 7.04×10−4 and 2.74×10−3 substitutions per site per year, respectively. See Fig. 1 for a break down of these rates by each gene. It should be noted that the NS2 rates are inflated, likely due to the unique placement of a type 1b isolate (V060) in the type 1a clade and ambiguous placement of V083 as seen using Densitree when compared to the same figure based on E2. Recombination between viral isolates may also explain these unusual findings. The remainder of the NS2 gene tree is consistent with other gene trees. Excluding NS2, the evolutionary rates of the first and second codon positions of E2 appear to be the only ones notably higher than their counterparts in other genes. The RB model implemented in BEAST 2 was employed to both select substitution models and to infer phylogenies. Fig. 2 indicates the proportion of the BEAST 2 run each codon was evaluated using each substitution model for each gene, excluding those that did not properly converge as described above. Across all genes analyzed, codon position three tended to most often be best modelled using the SYM model while positions one and two were more variable depending on the gene being analyzed. Phylogenetic trees were built using the full post-burnin set of trees even though they were derived from multiple substitution models explored during the analysis. Fig. 3 shows the phylogenetic tree based on the E2 gene. The clusters based on BVDV1a, BVDV1b and BVDV2a are clearly visible. There are also several more contemporary clusters that are not defined at the sub-genotype level. A notable example are the two clusters within the BVDV2a isolates that are consistently present regardless of the phylogeny being examined. For reference, an additional phylogeny using E2 sequences but also including several common reference strains was prepared (Fig. S1). The remaining gene-specific phylogenies are in supplemental (Figs. S2–S10). The phylogenetic network in Fig. 4 visualizes the differences and similarities between the phylogenies of each individual gene. Lines and narrow reticulations indicate agreement among the trees. The general topology of all gene trees is quite similar with some disagreement among BVDV1a isolates. Likewise, the pairwise distances between taxa calculated using different gene trees are always significantly correlated 234
Virology 509 (2017) 232–238
A. Chernick, F. van der Meer
sequences identified many differences. The consensus sequence from all BVDV1a E2 sequences differed at 47 amino acid positions from the Singer reference sequence (GenBank accession DQ088995) while the BVDV1b E2 consensus sequence differed from the CC13B reference sequence (GenBank accession KF772785) at 23 amino acid positions (Figs. S11 and S12). This variability appears evenly distributed across the gene with the exception of a six amino acid cluster starting at amino acid 30 in the gene where the Singer reference sequence encodes LEGLLA and the BVDV1a consensus sequence encodes WKDYSH. On mature viral particles this region is most distal from the viral surface on the E2 protein, accessible to the host immune system. Additionally, the consensus sequence alignment of the two BVDV2a clusters noted above identified 21 amino acid differences between those clusters in the E2 gene. A phylogenetic analysis of these BVDV2a isolates with BVDV2a and BVDV2b 5′UTR reference sequences also indicated that all of the strains detected in Canada cluster with BVDV2a; the clustering described in this data set occurs below the sub-genotype level.
4. Discussion This study expands upon and supports previous work on the evolution of BVDV in Canada (Chernick et al., 2014). The evolutionary rate of BVDV is 1.4×10−3 substitutions per site per year in this data set which is consistent with previous findings and with ssRNA viruses in general (Yuan et al., 2013; Suzuki et al., 2000; Duffy et al., 2008). Our analysis also indicates that the collection of BVDV isolates that can be found circulating in Canadian cattle is dynamic and has changed with respect to common reference strains based on sequence alignments. Although there is likely some bias in these samples based on the convenience sampling method, this finding still strongly supports the need for ongoing herd sampling and evaluation/revision of BVDV vaccines. Although some strains currently in circulation share similarities with older reference strains (as shown by the consensus sequence alignments) and there is some cross protection between vaccines against genetically distinct strains (Xue et al., 2010), there can also
Fig. 2. For each gene the proportion of time that the Bayesian phylogenetic analysis spent using each possible substitution model at each codon position. This is indicative of how well the model fit that particular data set. Note that the Felsenstein (F81) model was never employed so is not noted in the legend even though it is part of the RB model.
with one another. This is true even in cases where the pairwise scatter plots indicate some disagreement between gene trees, particularly for NS2 and NS5B (Fig. 5), both of which may contain recombinant sequences. Sequence alignments to compare the cluster-specific consensus sequences from these isolates to one another and to reference
Fig. 3. Phylogenetic tree of all isolates using E2 gene sequences. Branches are shaded based on their evolutionary rate (min=blue, max=red). The symbols beside each taxa are associated with the province of origin for that viral isolate (line=BC, triangle=AB, square=SK, hexagon=MB, circle=ON, star=QB). The vertical bars are spaced apart by 50 years each.
235
Virology 509 (2017) 232–238
A. Chernick, F. van der Meer
Fig. 4. A splits network consisting of all available gene trees. The BVDV1a, BVDV1b and BVDV2a sub-genotypes are indicated.
due to a lack of phylogenetic signal as indicated by bimodal and unconverged posterior distributions. It appears that in order to properly assess a BVDV phylogeny of this nature using NS2 and NS4A gene sequences sampling over a much larger time period would be required in order to observe a sufficient number of genetic changes. The selection method employed by the RB substitution model found that different models better fit the data at each codon position and at each gene (Fig. 2). It is well documented that the tolerance for nucleotide changes at specific codon positions varies due to how most of the degeneracy in the codon table is at the third codon position (Crick, 1966). Also, lower substitution rates are better suited to infer more ancestral events over larger sampling periods as indicated by the lack of convergence for genes with lower substitution rates in this analysis (NS2 and NS4A, specifically). Higher substitution rates, such as those generally observed at the third codon position and other genes in this analysis, allow for a finer grained phylogenetic inference of contemporary evolutionary events over smaller sampling periods. This is especially true within closely related isolates within sub-genotypes. This does not mean that the codon positions with lower substitution rates are experiencing a different evolutionary history, rather that the history cannot be fully observed at those positions. Furthermore, the gene-specific variation in preferred substitution models results from genes being subjected to unique selective pressures. Although they are all part of the same genome, gene-specific selective pressures cause them to each evolve in a different manner and therefore be best described using gene-specific substitution models. Also of note is the lack of geographic clustering within the phylogenetic trees based on visual inspection of the trees. The small clusters of isolates from the same province are generally also from the same farm and represent farm-specific strains. This phenomenon has been previously described using 5′UTR sequences from BVDV isolates in Sweden (Vilcek et al., 1999). As in this previous research, when these farm-specific patterns are excluded from the analysis there is no evidence that BVDV is geographically restricted. We propose this to be a result of the extensive transportation of the Canadian cattle population within Canada and to the United States and the ensuing mixing of viral isolates (Dube et al., 2010; Gonzalez et al., 2012). Epidemiological studies in Scotland have emphasized the importance of cattle movement in BVDV transmission (Gates et al., 2013). Additionally, it has also been found that cattle movement in the UK is an important contributor to BVDV transmission (Booth et al., 2013). Although more intensive sampling may elucidate some phylogeographical patterns in Canada, this study cannot support any such interpretation. We therefore conclude that any interventions to limit the circulation of BVDV in livestock populations should be implemented at a national level. As detailed in the results, there are phylogenetic population
Fig. 5. A matrix of scatter plots based on the phylogenetic distances between each pair of viral isolates calculated using each gene tree. NS2 and NS4A have been included in this plot to indicate the consequences of poor convergence and the incorrect placement of taxa during phylogenetic inference.
be significant antigenic differences (Ridpath et al., 2010). For example, the six amino acid stretch of changes (LEGLLA to WKDYSH) identified in the BVDV1a isolates represents a notable deviation from the Singer strain which is commonly used in vaccine production. The difference at this location between the NADL reference strain and these samples is more subtle (WKEYSPG to WKDYSHE) and includes a seventh position, but is still important. Since these amino acid changes occurred in an immunologically accessible region of this protein it can be expected that the antigenic properties of these isolates will differ from those of the reference strain. This and other changes like it may lead to reduced vaccine efficacy over time and eventual vaccine escape. The lack of robust phylogenetic inference in the NS2 and NS4A genes does not negate the findings based on other viral genes. Genes that are highly conserved are not suitable for use in a phylogenetic analysis over a relatively short sampling period. A previous study was not able to use NS2 in a Bayesian phylogenetic approach either due to a lack of variability (Xia et al., 2007). Similar to their observations, applying our approach did not identify a specific evolutionary pattern
236
Virology 509 (2017) 232–238
A. Chernick, F. van der Meer
Appendix A. Supporting information
structures below the level of sub-genotype in this data that are observed in all gene-specific trees. Of particular interest are the two clusters of BVDV2a isolates. Other work has proposed the segregation of BVDV2 into two sub-genotypes (Flores et al., 2002). Based on this classification all of the BVDV2 isolates in this study would be classified as BVDV2a using partial 5′UTR sequences, indicating that all variation in Canadian BVDV2 isolates are below the sub-genotypic, BVDV2a level. Although this study did not assess the antigenic differences between these clusters in vitro, it did identify at least one amino acid difference between BVDV2a clusters in a site that has a known antigenic impact in BVDV1 viruses. Specifically, a change from tyrosine to isoleucine at amino acid position 963 (relative to the Oregon reference strain, GenBank AF091605) reduced or eliminated binding of three different monoclonal antibodies (Paton et al., 1992). In this data set the groups are differentiated at the same amino acid position by a lysine in one group and an asparagine in the other. It should also be noted that all of the variation in these isolates is below the subgenotype level. It is also likely that some of the other mutations that exist across the genome result in different antigenic properties between these clusters. Other characteristics including viral enzyme stability and fidelity, transmissibility, host-cell binding affinity and so on could be impacted in a similar manner. Introducing more fine-grained levels of classification may prove to be useful but should be supported by additional research that phenotypically characterizes and differentiates these groups. In conclusion, this study provides one of the most comprehensive phylogenetic studies of BVDV in Canada in terms of the number of samples included. We observed diversity which is the result of ongoing temporal changes in the BVDV isolates in circulation and were unable to demonstrate a spatial element to this variability. These analyses provide insight into how rapidly evolving ssRNA viruses can evolve. Additional research that more explicitly links the small-scale evolutionary processes in persistent infections to the long-term evolution of BVDV would provide further insight while more rigorous, ongoing sampling of the isolates in circulation could help to confirm or reject the lack of spatial structure in the phylogenies presented here.
Supplementary data associated with this article can be found in the online version at doi:10.1016/j.virol.2017.06.024. References Booth, R.E., Thomas, C.J., El-Attar, L.M., Gunn, G., Brownlie, J., 2013. A phylogenetic analysis of Bovine Viral Diarrhoea Virus (BVDV) isolates from six different regions of the UK and links to animal movement data. Vet. Res. 44 (1), 43. Bouckaert, R., Alvarado-Mora, M.V., Pinho, J.R., 2013. Evolutionary rates and HBV: issues of rate estimation with Bayesian molecular methods. Antivir. Ther. 18 (3 Pt B), 497–503. Bouckaert, R., Heled, J., Kuhnert, D., Vaughan, T., Wu, C.H., Xie, D., et al., 2014. BEAST 2: a software platform for Bayesian Evolutionary Analysis. PLoS Comput. Biol. 10 (4), e1003537. Chernick, A., Godson, D.L., van der Meer, F., 2014. Metadata beyond the sequence enables the phylodynamic inference of bovine viral diarrhea virus type 1a isolates from Western Canada. Infect. Genet. Evol.: J. Mol. Epidemiol. Evolut. Genet. Infect. Dis. 28, 367–374. Couvreur, B., Letellier, C., Collard, A., Quenon, P., Dehan, P., Hamers, C., et al., 2002. Genetic and antigenic variability in bovine viral diarrhea virus (BVDV) isolates from Belgium. Virus Res. 85 (1), 17–28. Crick, F.H., 1966. Codon–anticodon pairing: the wobble hypothesis. J. Mol. Biol. 19 (2), 548–555. Delport, W., Poon, A.F., Frost, S.D., Kosakovsky Pond, S.L., 2010. Datamonkey 2010: a suite of phylogenetic analysis tools for evolutionary biology. Bioinformatics 26 (19), 2455–2457. Deregt, D., Masri, S.A., Cho, H.J., Ohmann, H.B., 1990. Monoclonal antibodies to the P80/125 and Gp53 proteins of bovine viral diarrhea virus: their potential use as diagnostic reagents. Can. J. Vet. Res. 54 (3), 343–348. Deregt, D., Smithson, S., Kozub, G.C., 1992. A short incubation serum neutralization test for bovine viral diarrhea virus. Can. J. Vet. Res. 56 (2), 161–164. Deregt, D., Bolin, S.R., van den Hurk, J., Ridpath, J.F., Gilbert, S.A., 1998a. Mapping of a type 1-specific and a type-common epitope on the E2 (gp53) protein of bovine viral diarrhea virus with neutralization escape mutants. Virus Res. 53 (1), 81–90. Deregt, D., van Rijn, P.A., Wiens, T.Y., van den Hurk, J., 1998b. Monoclonal antibodies to the E2 protein of a new genotype (type 2) of bovine viral diarrhea virus define three antigenic domains involved in neutralization. Virus Res. 57 (2), 171–181. Dow, N., Chernick, A., Orsel, K., van Marle, G., van der Meer, F., 2015. Genetic variability of Bovine Viral Diarrhea Virus and evidence for a possible genetic bottleneck during vertical transmission in persistently infected cattle. PloS One 10 (7), e0131972. Dube, C., Ribble, C., Kelton, D., 2010. An analysis of the movement of dairy cattle through 2 large livestock markets in the province of Ontario, Canada. Can. Vet. J. 55 (11), 1254–1260. Dubovi, E.J., 1992. Genetic diversity and BVD virus. Comp. Immunol., Microbiol. Infect. Dis. 15 (3), 155–162. Duffy, S., Shackelton, L.A., Holmes, E.C., 2008. Rates of evolutionary change in viruses: patterns and determinants. Nat. Rev. Genet. 9 (4), 267–276. Edgar, R.C., 2004. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32 (5), 1792–1797. Felsenstein, J., 1981. Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17 (6), 368–376. Fletcher, S.P., Jackson, R.J., 2002. Pestivirus internal ribosome entry site (IRES) structure and function: elements in the 5′ untranslated region important for IRES function. J. Virol. 76 (10), 5024–5033. Flores, E.F., Ridpath, J.F., Weiblen, R., Vogel, F.S.F., Gil, L.H.V.G., 2002. Phylogenetic analysis of Brazilian bovine viral diarrhea virus type 2 (BVDV-2) isolates: evidence for a subgenotype within BVDV-2. Virus Res. 87 (1), 51–60. Gates, M.C., Woolhouse, M.E., Gunn, G.J., Humphry, R.W., 2013. Relative associations of cattle movements, local spread, and biosecurity with bovine viral diarrhoea virus (BVDV) seropositivity in beef and dairy herds. Prev. Vet. Med. 112 (3–4), 285–295. Gilbert, S.A., Burton, K.M., Prins, S.E., Deregt, D., 1999. Typing of bovine viral diarrhea viruses directly from blood of persistently infected cattle by multiplex PCR. J. Clin. Microbiol. 37 (6), 2020–2023. Gonzalez, L.A., Schwartzkopf-Genswein, K.S., Bryan, M., Silasi, R., Brown, F., 2012. Benchmarking study of industry practices during commercial long haul transport of cattle in Alberta, Canada. J. Anim. Sci. 90 (10), 3606–3617. Greiser-Wilke, I., Haas, L., Dittmar, K., Liess, B., Moennig, V., 1993. RNA insertions and gene duplications in the nonstructural protein p125 region of pestivirus strains and isolates in vitro and in vivo. Virology 193 (2), 977–980. Harada, T., Tautz, N., Thiel, H.J., 2000. E2-p7 region of the bovine viral diarrhea virus polyprotein: processing and functional studies. J. Virol. 74 (20), 9498–9506. Hasegawa, M., Kishino, H., Yano, T.-a., 1985. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22 (2), 160–174. Iqbal, M., Flick-Smith, H., McCauley, J.W., 2000. Interactions of bovine viral diarrhoea virus glycoprotein E(rns) with cell surface glycosaminoglycans. J. Gen. Virol. 81 (Pt 2), 451–459. Isken, O., Langerwisch, U., Schonherr, R., Lamp, B., Schroder, K., Duden, R., et al., 2014. Functional characterization of bovine viral diarrhea virus nonstructural protein 5A by reverse genetic analysis and live cell imaging. J. Virol. 88 (1), 82–98. Jackova, A., Novackova, M., Pelletier, C., Audeval, C., Gueneau, E., Haffar, A., et al., 2008. The extended genetic diversity of BVDV-1: typing of BVDV isolates from France. Vet. Res. Commun. 32 (1), 7–11.
Author statements Funding Information University of Calgary, Faculty of Veterinary Medicine (10004757). University of Calgary, Margaret Gunn Endowment for Animal Research (10005041). Alberta Livestock and Meat Agency (2011F096R). Conflicts of interest The authors declare that they hold no conflicts of interest regarding this research. Ethical statement All viral isolates collected from animals in this study were collected in accordance with the approved animal care protocol VSACC AC110081/AC14-0144 at the University of Calgary. Acknowledgments We would like to acknowledge the invaluable contribution of samples from the following groups: John Robinson and Melissa Trapp at the Plant and Animal Health Centre, BC Ministry of Agriculture, Dale Godson at Prairie Diagnostic Services, Carl Gagnon at Université de Montreal, and Davor Ojkic at the University of Guelph. Additionally, we thank Sandeep Atwal and Gurleen Brar for their assistance with sample processing. 237
Virology 509 (2017) 232–238
A. Chernick, F. van der Meer
pestiviruses: determination of cleavage sites. J. Virol. 71 (7), 5415–5422. Uttenthal, A., Stadejek, T., Nylin, B., 2005. Genetic diversity of bovine viral diarrhoea viruses (BVDV) in Denmark during a 10-year eradication period. APMIS 113 (7–8), 536–541. Vilcek, S., Alenius, S., Paton, D.J., Mittelholzer, C., Belak, S., 1999. Genetic clustering of bovine viral diarrhoea viruses in cattle farms: genetic identification and analysis of viruses directly from cattle sera. Vet. J. 158 (1), 33–38. Vilcek, S., Durkovic, B., Kolesarova, M., Greiser-Wilke, I., Paton, D., 2004. Genetic diversity of international bovine viral diarrhoea virus (BVDV) isolates: identification of a new BVDV-1 genetic group. Vet. Res. 35 (5), 609–615. Vilček, Š., Paton, D.J., Durkovic, B., Strojny, L., Ibata, G., Moussa, A., et al., 2001. Bovine viral diarrhoea virus genotype 1 can be separated into at least eleven genetic groups. Arch. Virol. 146 (1), 99–115. Warrener, P., Collett, M.S., 1995. Pestivirus NS3 (p80) protein possesses RNA helicase activity. J. Virol. 69 (3), 1720–1726. Weiskircher, E., Aligo, J., Ning, G., Konan, K.V., 2009. Bovine viral diarrhea virus NS4B protein is an integral membrane protein associated with Golgi markers and rearranged host membranes. Virol. J. 6, 185. Wiskerchen, M., Belzer, S.K., Collett, M.S., 1991. Pestivirus gene expression: the first protein product of the bovine viral diarrhea virus large open reading frame, p20, possesses proteolytic activity. J. Virol. 65 (8), 4508–4514. Xia, H., Liu, L., Wahlberg, N., Baule, C., Belak, S., 2007. Molecular phylogenetic analysis of bovine viral diarrhoea virus: a Bayesian approach. Virus Res. 130 (1–2), 53–62. Xu, J., Mendez, E., Caron, P.R., Lin, C., Murcko, M.A., Collett, M.S., et al., 1997. Bovine viral diarrhea virus NS3 serine proteinase: polyprotein cleavage sites, cofactor requirements, and molecular model of an enzyme essential for pestivirus replication. J. Virol. 71 (7), 5312–5322. Xue, W., Mattick, D., Smith, L., Umbaugh, J., Trigo, E., 2010. Vaccination with a modified-live bovine viral diarrhea virus (BVDV) type 1a vaccine completely protected calves against challenge with BVDV type 1b strains. Vaccine 29 (1), 70–76. Yu, H., Grassmann, C.W., Behrens, S.E., 1999. Sequence and structural elements at the 3′ terminus of bovine viral diarrhea virus genomic RNA: functional role during RNA replication. J. Virol. 73 (5), 3638–3648. Yuan, M., Lu, T., Li, C., Lu, L., 2013. The evolutionary rates of HCV estimated with subtype 1a and 1b sequences over the ORF length and in different genomic regions. PloS One 8 (6), e64698. Zharkikh, A., 1994. Estimation of evolutionary distances between nucleotide sequences. J. Mol. Evol. 39 (3), 315–329. Zhong, W., Gutshall, L.L., Del Vecchio, A.M., 1998. Identification and characterization of an RNA-dependent RNA polymerase activity within the nonstructural protein 5B region of bovine viral diarrhea virus. J. Virol. 72 (11), 9365–9369.
Kearse, M., Moir, R., Wilson, A., Stones-Havas, S., Cheung, M., Sturrock, S., et al., 2012. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28 (12), 1647–1649. Kosakovsky Pond, S.L., Frost, S.D., Muse, S.V., 2005. HyPhy: hypothesis testing using phylogenies. Bioinformatics 21 (5), 676–679. Kosakovsky Pond, S.L., Posada, D., Gravenor, M.B., Woelk, C.H., Frost, S.D., 2006. Automated phylogenetic detection of recombination using a genetic algorithm. Mol. Biol. Evol. 23 (10), 1891–1901. Lanave, C., Preparata, G., Saccone, C., Serio, G., 1984. A new method for calculating evolutionary substitution rates. J. Mol. Evol. 20 (1), 86–93. Liu, L., Xia, H., Wahlberg, N., Belak, S., Baule, C., 2009. Phylogeny, classification and evolutionary insights into pestiviruses. Virology 385 (2), 351–357. Luzzago, C., Ebranati, E., Sassera, D., Lo Presti, A., Lauzi, S., Gabanelli, E., et al., 2012. Spatial and temporal reconstruction of bovine viral diarrhea virus genotype 1 dispersion in Italy. Infect., Genet. Evol.: J. Mol. Epidemiol. Evolut. Genet. Infect. Dis. 12 (2), 324–331. Paton, D.J., Lowings, J.P., Barrett, A.D., 1992. Epitope mapping of the gp53 envelope protein of bovine viral diarrhea virus. Virology 190 (2), 763–772. Ridpath, J.F., 2010. Bovine viral diarrhea virus: global status. Vet. Clin. North Am. Food Anim. Pract. 26 (1), 105–121. Ridpath, J.F., Bolin, S.R., 1998. Differentiation of types 1a, 1b and 2 bovine viral diarrhoea virus (BVDV) by PCR. Mol. Cell Probes. 12 (2), 101–106. Ridpath, J.F., Neill, J.D., Frey, M., Landgraf, J.G., 2000. Phylogenetic, antigenic and clinical characterization of type 2 BVDV from North America. Vet. Microbiol. 77 (1– 2), 145–155. Ridpath, J.F., Fulton, R.W., Kirkland, P.D., Neill, J.D., 2010. 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. 22 (2), 184–191. Ronecker, S., Zimmer, G., Herrler, G., Greiser-Wilke, I., Grummer, B., 2008. Formation of bovine viral diarrhea virus E1-E2 heterodimers is essential for virus entry and depends on charged residues in the transmembrane domains. J. Gen. Virol. 89 (Pt 9), 2114–2121. Rumenapf, T., Unger, G., Strauss, J.H., Thiel, H.J., 1993. Processing of the envelope glycoproteins of pestiviruses. J. Virol. 67 (6), 3288–3294. Suzuki, Y., Yamaguchi-Kabata, Y., Gojobori, T., 2000. Nucleotide substitution rates of HIV-1. AIDS Rev. 2, 39–47. Tamura, K., Nei, M., 1993. Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Mol. Biol. Evol. 10 (3), 512–526. Tautz, N., Elbers, K., Stoll, D., Meyers, G., Thiel, H.J., 1997. Serine protease of
238
Contents lists available at ScienceDirect
Virology journal homepage: www.elsevier.com/locate/yviro
Evolution of Bovine viral diarrhea virus in Canada from 1997 to 2013 ⁎
Adam Chernick , Frank van der Meer
MARK
⁎
Ecosystem and Public Health, Faculty of Veterinary Medicine, University of Calgary, 3330 Hospital Dr NW, Calgary, AB T2N 1N4, Canada
A R T I C L E I N F O
A BS T RAC T
Keywords: Bovine viral diarrhea virus Evolution Next-generation sequencing ssRNA viruses Phylogenetics Canada
Bovine viral diarrhea virus (BVDV) is a rapidly evolving, single-stranded RNA virus and a production limiting pathogen of cattle worldwide. 79 viral isolates collected between 1997 and 2013 in Canada were subjected to next-generation sequencing. Bayesian phylogenetics was used to assess the evolution of this virus. A mean substitution rate of 1.4×10−3 substitutions/site/year was found across both BVDV1 and BVDV2. Evolutionary rates in the E2 gene were slightly faster than other regions. We also identified population structures below the sub-genotype level that likely have phenotypic implications. Two distinct clusters within BVDV2a are present and can be differentiated, in part, by a tyrosine to isoleucine mutation at position 963 in the E2 protein, a position implicated in the antigenicity of BVDV1 isolates. Distinct clustering within all sub-genotypes, particularly BVDV2a, is apparent and could lead to new levels of genotypic classification. Continuous monitoring of emerging variants is therefore necessary.
1. Introduction Bovine viral diarrhea virus (BVDV) is a single-stranded RNA (ssRNA) virus that is a production limiting pathogen of cattle with a global distribution (Ridpath, 2010). It can be categorized in two genotypes or species, BVDV1 and BVDV2. Numerous sub-genotypes of BVDV1 have been reported as well as two proposed BVDV2 subgenotypes (Vilček et al., 2001; Flores et al., 2002). In Canada and the United States BVDV1a, BVDV1b and BVDV2a are currently in circulation. Previous studies have described notable genetic diversity within and between each of these sub-genotypes (Vilcek et al., 2004) as well as within BVDV1a in western Canada (Chernick et al., 2014). In the latter study, the substitution rate of BVDV1a isolates was estimated to be 1.26×10−3 substitutions/site/year in the E1-E2 gene across all codon positions. This rapid evolution drives many important phenotypic changes in BVDV including extensive antigenic variability (Ridpath et al., 2000). High substitution rates are not unique to BVDV and are often noted in other ssRNA viruses including type 1a Hepatitis C virus (Yuan et al., 2013), HIV (Suzuki et al., 2000) and others (Duffy et al., 2008). This rapid evolution stems largely from the inherent properties of these viruses and their error-prone replication. In the case of BVDV specifically, the unique immune environment associated with persistent infections also appears to play an important role in generating novel viral variants (Dow et al., 2015).
BVDV is a Pestivirus within the family Flaviviridae. It shares a genome structure similar to these related viruses with 5′ and 3′ untranslated regions (UTRs), structural genes towards the 5′ end and most nonstructural genes to the 3′ end. The UTRs are both highly conserved with the 5′UTR acting as a ribosomal entry site (Fletcher and Jackson, 2002) and 3′UTR also essential to replication (Yu et al., 1999). Following the 5′UTR is the NPro, a viral protease responsible for cleaving some viral proteins from one another (Wiskerchen et al., 1991). The structural proteins follow and include the envelope proteins (Erns, E1 and E2) and the capsid. Erns is secreted from infected cells (Rumenapf et al., 1993) and binds to host cell surfaces (Iqbal et al., 2000) while E1 and E2 are both surface glycoproteins important to host cell binding and entry (Ronecker et al., 2008) as well as being highly antigenic (Paton et al., 1992; Deregt et al., 1998a, 1998b). The p7 gene sits between the structural and non-structural genes and is required for infectivity (Harada et al., 2000). NS2 and NS3 are non-structural genes encoding for a protein that performs multiple functions including a serine protease (Tautz et al., 1997) and helicase (Warrener and Collett, 1995). NS4A is a co-factor for the NS2/3 serine protease activity (Xu et al., 1997) and NS4B appears to play a role in viral genome replication and host membrane remodelling (Weiskircher et al., 2009). NS5A (Isken et al., 2014) and NS5B (Zhong et al., 1998) are both components of the RNA-dependent RNA polymerase complex. Several phylogenetic and genetic diversity studies of pestiviruses
Abbreviations: BVDV, Bovine viral diarrhea virus; E, envelope; ESS, effective sample size; Gbp, gigabase pairs; GTR, general time-reversible; HIV, human immunodeficiency virus; HPD95, 95% highest probability density; MRCA, most recent common ancestor; NS, non-structural; RB, reversible-jump based; ssRNA, single-stranded RNA; SYM, symmetrical; TIM, transition; UK, United Kingdom; UTR, untranslated region; VSACC, Veterinary Sciences Animal Care Committee ⁎ Corresponding authors. E-mail addresses: [email protected] (A. Chernick), [email protected] (F. van der Meer). http://dx.doi.org/10.1016/j.virol.2017.06.024 Received 8 April 2017; Received in revised form 20 June 2017; Accepted 21 June 2017 0042-6822/ © 2017 Elsevier Inc. All rights reserved.
Virology 509 (2017) 232–238
A. Chernick, F. van der Meer
and BVDV have also been conducted, although none exhaustively examine Canadian field isolates. With respect to the age of the entire BVDV lineage, all BVDV1 and BVDV2 isolates are believed the share a most recent common ancestor about the year 1743 (95% highest probability density (HPD95) 1373–1926) (Liu et al., 2009). Other studies have examined the spatial movement of BVDV while inferring evolutionary rates about 9.3×10−3 substitutions/site/year in the 5′ untranslated region (5′UTR) among Italian isolates (Luzzago et al., 2012). Phylogenetic and genetic diversity studies of European isolates are available in large part due to the need to inform BVDV-eradication and control programs in those countries. For example, isolates in Denmark span BVDV1 sub-genotypes a, b, d, e, f, g and h (Uttenthal et al., 2005) while a collection from France included sub-genotypes b, d, e and l (Jackova et al., 2008). Prior to these phylogenetic studies numerous experiments were conducted using anti-BVDV antibodies to assess strain variation. Although many antibodies were found to be broadly neutralizing against genetically distinct isolates (Deregt et al., 1990), significant differences that correlated with genetic distance were also noted (Deregt et al., 1992). Panels of antibodies capable of distinguishing BVDV1 from BVDV2 as well as some sub-genotypic variation have also been generated (Deregt et al., 1998a, 1998b). Although there is merit in classifying isolates based on antibody binding assays, it has been noted that such a practice is “suggestive at best”, implying that sequence data is a more reliable tool even if it does not produce such an explicit connection to the phenotypic/antigenic properties of a given isolate (Dubovi, 1992). Sequence data also allows for more comprehensive evolutionary studies that describe the origins and evolutionary rates of BVDV. In this study we describe the evolution of BVDV in Canada using a collection of viral isolates archived from 1997 to 2013 and obtained from academic and diagnostic laboratories across the country as well as from active field sampling efforts. We employed a next generation sequencing (NGS) approach in conjunction with Bayesian phylogenetics to generate an integrated model of BVDV1 and BVDV2 evolution within the Canadian cattle population. This both describes the evolutionary history of this pathogen in Canada as well as important variability among currently circulating BVDV isolates. Apart from the practical use of this information in the development and maintenance of virus control programs including vaccination, these results also contribute to the understanding of the evolutionary capacity of ssRNA viruses which include a number of important human and animal pathogens.
Table 1 Primer sequences and references for the four amplicon method used to generate DNA for next-generation sequencing. Name
Sequence
Reference
Position
Frag1_fwd
CATGCCCATAGTAGGAC
107
Frag1_rev
GGGCAWACCATYTGGAAGGCYGG
Frag2_fwd
TGATAACAGGGGTACAAGGGC
Frag2_rev Frag3_fwd
RTTGYCCCTAGCGGTATAT GCAGATTTTGAAGAAAGACACTA
Frag3_rev
AAGAAGCCATCATCMCCACA
Frag4_fwd
AAGATCCACCCTTATGARGC
Frag4_rev
CTGTGTGCATTRARTGTAGTGTT
(Harada et al., 2000) (Tautz et al., 1997) (Warrener and Collett, 1995) N/A (Warrener and Collett, 1995) (Xu et al., 1997) (Xu et al., 1997) N/A
2853 2442 5314 4934 11547 10385 12511
Primer sequences and references for the four amplicon method used to generate DNA for next-generation sequencing. All positions are relative to the NADL reference genome (GenBank Accession M31182).
subjected to library preparation using the Nextera XT DNA Library Preparation kit (Illumina, San Diego, USA). All libraries were prepared with dual indices and sequenced with paired-end, 250 bp reads on a MiSeq using a 500-cycle v2 cartridge (Illumina, San Diego, USA). 2.2. Data Analysis Raw reads were imported into Geneious R9 (Kearse et al., 2012) for pre-processing and assembly. Reads were paired and assembled using the Geneious assembler set to high sensitivity and iterated up to 5 times. Consensus sequences were extracted from assemblies for all individual genes except p7 (leaving a total of 10 genes). These are available in GenBank via accession numbers KX169935 through KX170703 along with the relevant protein translations (Table S2). Gene-specific consensus sequence alignments were built using MUSCLE (Edgar, 2004). Alignments were visually inspected and corrected for errors. The HyPhy software package (Kosakovsky Pond et al., 2005) running on the Datamonkey server (Delport et al., 2010) was used to test for recombination. Specifically, the single breakpoint recombination (SBP) analysis (Kosakovsky Pond et al., 2006) was run using the HKY85 substitution model on each gene-specific alignment. Sequence alignments were imported into BEAUti to prepare the run parameter files for BEAST 2 (Bouckaert et al., 2014). Viral genes were analyzed independently to explicitly accommodate any differences in their phylogenies. All sequences were dated with their year of collection and separated into three partitions for the codon positions. Year was used as it was the most precise date data available for all isolates. A reversible-jump based (RB) substitution model accounting for ambiguities in the consensus sequences was employed to simultaneously explore different substitution models and infer the phylogeny using the mixture that best fit the data (Bouckaert et al., 2014, 2013). Six substitution models were evaluated including the Felsenstein (F81) (Felsenstein, 1981), HKY (Hasegawa et al., 1985), Tamura and Nei (TN93) (Tamura and Nei, 1993), transition (TIM), symmetrical (SYM) (Zharkikh, 1994) and generalized time-reversible (GTR) (Lanave et al., 1984) models. Separate relaxed exponential molecular clocks were used for each codon partition. Runs were allowed to continue until adequate convergence (ESS > 200 for most parameters) was reached. Phylogenetic trees were built using FigTree v1.4.2 (http://tree.bio.ed. ac.uk/software/figtree/) and Tracer was used to extract parameters of interest while verifying run quality. Two additional trees were built. The first using only the BVDV2a isolates in addition to BVDV2a and BVDV2b reference strains to further refine the clusters observed in the larger trees using partial 5′UTR sequences. The second includes all isolates in this study as well
2. Methods 2.1. Sample collection and processing Viral isolates were obtained from diagnostic and research laboratories across Canada (Table S1). This was the result of an extensive survey of laboratories working on BVDV and likely represents most historical isolates available for evaluation. All samples were collected in accordance with the animal care protocol VSACC AC11-0081/AC140144. All isolates had been passaged at least once at the laboratory of origin in Madin Darby Bovine Kidney or Bovine Turbinate cell culture prior to further characterization in our laboratory. RNA was extracted from an aliquot of each isolate using the E.Z.N.A. Viral RNA kit (Omega Bio-tek, Norcross, USA). A series of four reverse transcriptase PCRs were used to amplify overlapping fragments of the protein-coding region of the genome. See Table 1 for primer sequences and amplicon regions (Ridpath and Bolin, 1998; Couvreur et al., 2002; Greiser-Wilke et al., 1993; Gilbert et al., 1999). The amplicons were gel purified using the E.Z.N.A. Gel Extraction kit (Omega Bio-tek, Norcross, USA) and quantified using the Qubit 3.0 and the Qubit dsDNA HS Assay kit (Thermo Fisher Scientific, Waltham, USA). Amplicons were pooled in equimolar amounts, adjusted for amplicon length. Each pool was 233
Virology 509 (2017) 232–238
A. Chernick, F. van der Meer
as a selection of reference strains to place these isolates within the global context based on E2 sequences. Both were built as described above but without tip dates (5′UTR and E2) or codon partitioning (5′UTR only). To compare different gene trees two approaches were employed. First, to assess the overall similarity of the trees both visual inspections as well as a phylogenetic network analysis were employed. Although phylogenetic networks are not as easily interpreted as standard phylogenetic trees they do provide a means to visualize phylogenetic uncertainty or the incongruity between multiple trees. Second, to assess differences in branch lengths, Pearson correlations between gene trees were calculated using the branch length distances between all pairs of taxa in SPSS v22 (IBM, Armonk, USA). In addition to the correlation coefficients, scatterplots were drawn to assess the relationship between branch lengths for each pair of genes. To compliment the tree-based assessment, sequence alignments were also generated to highlight differences between certain clusters. Consensus sequences for these clusters of interest were generated from MUSCLE sequence alignments and compared to consensus sequences derived from other clusters and to older reference strain sequences from GenBank. These consensus sequences were allowed to contain ambiguities to account for any diversity that exists within subgenotypes.
Fig. 1. Mean clock rates with HPD95 intervals for each codon position in each gene phylogeny. Clock rates are measured in nucleotide substitutions per site per year. NS2 and NS4A have been excluded due to poor convergence during the Bayesian phylogenetic analysis.
3. Results The sequencing run generated a total of 4.49 Gbp across all samples and controls. The cluster density was 488 k/mm2 with 92.97% of all clusters passed quality filters. Based on re-sequencing of the PhiX control library a machine error rate of 2.04% was calculated. Across all clusters that passed the quality filter, an average of 74.94% of the resulting bases had a Q-score ≥ 30. Following assembly, it was found that full protein-coding region sequences were not obtained for all isolates due to inadequate coverage in certain regions, especially in the region amplified by the third RTPCR covering positions 4934–11547. This RT-PCR most commonly failed for BVDV2 isolates, indicating a possible genotypic bias in the primers. Excluding the p7 gene, consensus sequences for all 10 genes in the coding region were generated for 42 of 79 isolates. Nine out of 10 genes were sequenced successfully for another 10 isolates with the remaining 27 isolates yielding eight or fewer genes. See Table S1 for a break down of the genes sequenced for each isolate. Phylogenetic reconstruction was performed for each gene using all available isolate sequences. The trees for NPro, Erns, E1 and E2 contain the full set of 79 isolates. Regarding recombination, no signs of recombination were detected using the small sample AIC indicator in the SBP analysis except for in NS2 and NS5B. There was significant support for recombination in NS2 at nucleotide positions 374–375 in that alignment and NS5B at nucleotide positions 420–421 and 615 in that alignment. As the following phylogenetic analysis was performed on individual genes (with the exception of the phylogenetic network) no isolates were removed. The NS2 and NS5B phylogenies are for comparison only and cannot be used to infer the relatedness of these sequences. Bayesian evolutionary analysis was performed using BEAST 2 where each sequence alignment was assessed separately until satisfactory ESS ( > 200) was obtained for key parameters such as the posterior, tree likelihoods, tree height and clock rates. Notable exceptions were the NS2 and NS4A genes which demonstrated unconverged or bimodal distributions in several key parameters, respectively. The following statistics exclude NS2 and NS4A for this reason. The mean most recent common ancestor (MRCA) for all phylogenies was dated to 1498 with a mean HPD95 range of 1182–1750. All individual HPD95 ranges overlapped as well. For individual sub genotypes the mean MRCA (mean HPD95 range) across all genes were 1930 (1890–1966), 1955 (1915–1980) and 1950 (1902–1979) for BVDV1a, BVDV1b and
BVDV2a, respectively. The mean uncorrelated exponential clock rates across all sub-genotypes for each codon position were 7.57×10−4, 7.04×10−4 and 2.74×10−3 substitutions per site per year, respectively. See Fig. 1 for a break down of these rates by each gene. It should be noted that the NS2 rates are inflated, likely due to the unique placement of a type 1b isolate (V060) in the type 1a clade and ambiguous placement of V083 as seen using Densitree when compared to the same figure based on E2. Recombination between viral isolates may also explain these unusual findings. The remainder of the NS2 gene tree is consistent with other gene trees. Excluding NS2, the evolutionary rates of the first and second codon positions of E2 appear to be the only ones notably higher than their counterparts in other genes. The RB model implemented in BEAST 2 was employed to both select substitution models and to infer phylogenies. Fig. 2 indicates the proportion of the BEAST 2 run each codon was evaluated using each substitution model for each gene, excluding those that did not properly converge as described above. Across all genes analyzed, codon position three tended to most often be best modelled using the SYM model while positions one and two were more variable depending on the gene being analyzed. Phylogenetic trees were built using the full post-burnin set of trees even though they were derived from multiple substitution models explored during the analysis. Fig. 3 shows the phylogenetic tree based on the E2 gene. The clusters based on BVDV1a, BVDV1b and BVDV2a are clearly visible. There are also several more contemporary clusters that are not defined at the sub-genotype level. A notable example are the two clusters within the BVDV2a isolates that are consistently present regardless of the phylogeny being examined. For reference, an additional phylogeny using E2 sequences but also including several common reference strains was prepared (Fig. S1). The remaining gene-specific phylogenies are in supplemental (Figs. S2–S10). The phylogenetic network in Fig. 4 visualizes the differences and similarities between the phylogenies of each individual gene. Lines and narrow reticulations indicate agreement among the trees. The general topology of all gene trees is quite similar with some disagreement among BVDV1a isolates. Likewise, the pairwise distances between taxa calculated using different gene trees are always significantly correlated 234
Virology 509 (2017) 232–238
A. Chernick, F. van der Meer
sequences identified many differences. The consensus sequence from all BVDV1a E2 sequences differed at 47 amino acid positions from the Singer reference sequence (GenBank accession DQ088995) while the BVDV1b E2 consensus sequence differed from the CC13B reference sequence (GenBank accession KF772785) at 23 amino acid positions (Figs. S11 and S12). This variability appears evenly distributed across the gene with the exception of a six amino acid cluster starting at amino acid 30 in the gene where the Singer reference sequence encodes LEGLLA and the BVDV1a consensus sequence encodes WKDYSH. On mature viral particles this region is most distal from the viral surface on the E2 protein, accessible to the host immune system. Additionally, the consensus sequence alignment of the two BVDV2a clusters noted above identified 21 amino acid differences between those clusters in the E2 gene. A phylogenetic analysis of these BVDV2a isolates with BVDV2a and BVDV2b 5′UTR reference sequences also indicated that all of the strains detected in Canada cluster with BVDV2a; the clustering described in this data set occurs below the sub-genotype level.
4. Discussion This study expands upon and supports previous work on the evolution of BVDV in Canada (Chernick et al., 2014). The evolutionary rate of BVDV is 1.4×10−3 substitutions per site per year in this data set which is consistent with previous findings and with ssRNA viruses in general (Yuan et al., 2013; Suzuki et al., 2000; Duffy et al., 2008). Our analysis also indicates that the collection of BVDV isolates that can be found circulating in Canadian cattle is dynamic and has changed with respect to common reference strains based on sequence alignments. Although there is likely some bias in these samples based on the convenience sampling method, this finding still strongly supports the need for ongoing herd sampling and evaluation/revision of BVDV vaccines. Although some strains currently in circulation share similarities with older reference strains (as shown by the consensus sequence alignments) and there is some cross protection between vaccines against genetically distinct strains (Xue et al., 2010), there can also
Fig. 2. For each gene the proportion of time that the Bayesian phylogenetic analysis spent using each possible substitution model at each codon position. This is indicative of how well the model fit that particular data set. Note that the Felsenstein (F81) model was never employed so is not noted in the legend even though it is part of the RB model.
with one another. This is true even in cases where the pairwise scatter plots indicate some disagreement between gene trees, particularly for NS2 and NS5B (Fig. 5), both of which may contain recombinant sequences. Sequence alignments to compare the cluster-specific consensus sequences from these isolates to one another and to reference
Fig. 3. Phylogenetic tree of all isolates using E2 gene sequences. Branches are shaded based on their evolutionary rate (min=blue, max=red). The symbols beside each taxa are associated with the province of origin for that viral isolate (line=BC, triangle=AB, square=SK, hexagon=MB, circle=ON, star=QB). The vertical bars are spaced apart by 50 years each.
235
Virology 509 (2017) 232–238
A. Chernick, F. van der Meer
Fig. 4. A splits network consisting of all available gene trees. The BVDV1a, BVDV1b and BVDV2a sub-genotypes are indicated.
due to a lack of phylogenetic signal as indicated by bimodal and unconverged posterior distributions. It appears that in order to properly assess a BVDV phylogeny of this nature using NS2 and NS4A gene sequences sampling over a much larger time period would be required in order to observe a sufficient number of genetic changes. The selection method employed by the RB substitution model found that different models better fit the data at each codon position and at each gene (Fig. 2). It is well documented that the tolerance for nucleotide changes at specific codon positions varies due to how most of the degeneracy in the codon table is at the third codon position (Crick, 1966). Also, lower substitution rates are better suited to infer more ancestral events over larger sampling periods as indicated by the lack of convergence for genes with lower substitution rates in this analysis (NS2 and NS4A, specifically). Higher substitution rates, such as those generally observed at the third codon position and other genes in this analysis, allow for a finer grained phylogenetic inference of contemporary evolutionary events over smaller sampling periods. This is especially true within closely related isolates within sub-genotypes. This does not mean that the codon positions with lower substitution rates are experiencing a different evolutionary history, rather that the history cannot be fully observed at those positions. Furthermore, the gene-specific variation in preferred substitution models results from genes being subjected to unique selective pressures. Although they are all part of the same genome, gene-specific selective pressures cause them to each evolve in a different manner and therefore be best described using gene-specific substitution models. Also of note is the lack of geographic clustering within the phylogenetic trees based on visual inspection of the trees. The small clusters of isolates from the same province are generally also from the same farm and represent farm-specific strains. This phenomenon has been previously described using 5′UTR sequences from BVDV isolates in Sweden (Vilcek et al., 1999). As in this previous research, when these farm-specific patterns are excluded from the analysis there is no evidence that BVDV is geographically restricted. We propose this to be a result of the extensive transportation of the Canadian cattle population within Canada and to the United States and the ensuing mixing of viral isolates (Dube et al., 2010; Gonzalez et al., 2012). Epidemiological studies in Scotland have emphasized the importance of cattle movement in BVDV transmission (Gates et al., 2013). Additionally, it has also been found that cattle movement in the UK is an important contributor to BVDV transmission (Booth et al., 2013). Although more intensive sampling may elucidate some phylogeographical patterns in Canada, this study cannot support any such interpretation. We therefore conclude that any interventions to limit the circulation of BVDV in livestock populations should be implemented at a national level. As detailed in the results, there are phylogenetic population
Fig. 5. A matrix of scatter plots based on the phylogenetic distances between each pair of viral isolates calculated using each gene tree. NS2 and NS4A have been included in this plot to indicate the consequences of poor convergence and the incorrect placement of taxa during phylogenetic inference.
be significant antigenic differences (Ridpath et al., 2010). For example, the six amino acid stretch of changes (LEGLLA to WKDYSH) identified in the BVDV1a isolates represents a notable deviation from the Singer strain which is commonly used in vaccine production. The difference at this location between the NADL reference strain and these samples is more subtle (WKEYSPG to WKDYSHE) and includes a seventh position, but is still important. Since these amino acid changes occurred in an immunologically accessible region of this protein it can be expected that the antigenic properties of these isolates will differ from those of the reference strain. This and other changes like it may lead to reduced vaccine efficacy over time and eventual vaccine escape. The lack of robust phylogenetic inference in the NS2 and NS4A genes does not negate the findings based on other viral genes. Genes that are highly conserved are not suitable for use in a phylogenetic analysis over a relatively short sampling period. A previous study was not able to use NS2 in a Bayesian phylogenetic approach either due to a lack of variability (Xia et al., 2007). Similar to their observations, applying our approach did not identify a specific evolutionary pattern
236
Virology 509 (2017) 232–238
A. Chernick, F. van der Meer
Appendix A. Supporting information
structures below the level of sub-genotype in this data that are observed in all gene-specific trees. Of particular interest are the two clusters of BVDV2a isolates. Other work has proposed the segregation of BVDV2 into two sub-genotypes (Flores et al., 2002). Based on this classification all of the BVDV2 isolates in this study would be classified as BVDV2a using partial 5′UTR sequences, indicating that all variation in Canadian BVDV2 isolates are below the sub-genotypic, BVDV2a level. Although this study did not assess the antigenic differences between these clusters in vitro, it did identify at least one amino acid difference between BVDV2a clusters in a site that has a known antigenic impact in BVDV1 viruses. Specifically, a change from tyrosine to isoleucine at amino acid position 963 (relative to the Oregon reference strain, GenBank AF091605) reduced or eliminated binding of three different monoclonal antibodies (Paton et al., 1992). In this data set the groups are differentiated at the same amino acid position by a lysine in one group and an asparagine in the other. It should also be noted that all of the variation in these isolates is below the subgenotype level. It is also likely that some of the other mutations that exist across the genome result in different antigenic properties between these clusters. Other characteristics including viral enzyme stability and fidelity, transmissibility, host-cell binding affinity and so on could be impacted in a similar manner. Introducing more fine-grained levels of classification may prove to be useful but should be supported by additional research that phenotypically characterizes and differentiates these groups. In conclusion, this study provides one of the most comprehensive phylogenetic studies of BVDV in Canada in terms of the number of samples included. We observed diversity which is the result of ongoing temporal changes in the BVDV isolates in circulation and were unable to demonstrate a spatial element to this variability. These analyses provide insight into how rapidly evolving ssRNA viruses can evolve. Additional research that more explicitly links the small-scale evolutionary processes in persistent infections to the long-term evolution of BVDV would provide further insight while more rigorous, ongoing sampling of the isolates in circulation could help to confirm or reject the lack of spatial structure in the phylogenies presented here.
Supplementary data associated with this article can be found in the online version at doi:10.1016/j.virol.2017.06.024. References Booth, R.E., Thomas, C.J., El-Attar, L.M., Gunn, G., Brownlie, J., 2013. A phylogenetic analysis of Bovine Viral Diarrhoea Virus (BVDV) isolates from six different regions of the UK and links to animal movement data. Vet. Res. 44 (1), 43. Bouckaert, R., Alvarado-Mora, M.V., Pinho, J.R., 2013. Evolutionary rates and HBV: issues of rate estimation with Bayesian molecular methods. Antivir. Ther. 18 (3 Pt B), 497–503. Bouckaert, R., Heled, J., Kuhnert, D., Vaughan, T., Wu, C.H., Xie, D., et al., 2014. BEAST 2: a software platform for Bayesian Evolutionary Analysis. PLoS Comput. Biol. 10 (4), e1003537. Chernick, A., Godson, D.L., van der Meer, F., 2014. Metadata beyond the sequence enables the phylodynamic inference of bovine viral diarrhea virus type 1a isolates from Western Canada. Infect. Genet. Evol.: J. Mol. Epidemiol. Evolut. Genet. Infect. Dis. 28, 367–374. Couvreur, B., Letellier, C., Collard, A., Quenon, P., Dehan, P., Hamers, C., et al., 2002. Genetic and antigenic variability in bovine viral diarrhea virus (BVDV) isolates from Belgium. Virus Res. 85 (1), 17–28. Crick, F.H., 1966. Codon–anticodon pairing: the wobble hypothesis. J. Mol. Biol. 19 (2), 548–555. Delport, W., Poon, A.F., Frost, S.D., Kosakovsky Pond, S.L., 2010. Datamonkey 2010: a suite of phylogenetic analysis tools for evolutionary biology. Bioinformatics 26 (19), 2455–2457. Deregt, D., Masri, S.A., Cho, H.J., Ohmann, H.B., 1990. Monoclonal antibodies to the P80/125 and Gp53 proteins of bovine viral diarrhea virus: their potential use as diagnostic reagents. Can. J. Vet. Res. 54 (3), 343–348. Deregt, D., Smithson, S., Kozub, G.C., 1992. A short incubation serum neutralization test for bovine viral diarrhea virus. Can. J. Vet. Res. 56 (2), 161–164. Deregt, D., Bolin, S.R., van den Hurk, J., Ridpath, J.F., Gilbert, S.A., 1998a. Mapping of a type 1-specific and a type-common epitope on the E2 (gp53) protein of bovine viral diarrhea virus with neutralization escape mutants. Virus Res. 53 (1), 81–90. Deregt, D., van Rijn, P.A., Wiens, T.Y., van den Hurk, J., 1998b. Monoclonal antibodies to the E2 protein of a new genotype (type 2) of bovine viral diarrhea virus define three antigenic domains involved in neutralization. Virus Res. 57 (2), 171–181. Dow, N., Chernick, A., Orsel, K., van Marle, G., van der Meer, F., 2015. Genetic variability of Bovine Viral Diarrhea Virus and evidence for a possible genetic bottleneck during vertical transmission in persistently infected cattle. PloS One 10 (7), e0131972. Dube, C., Ribble, C., Kelton, D., 2010. An analysis of the movement of dairy cattle through 2 large livestock markets in the province of Ontario, Canada. Can. Vet. J. 55 (11), 1254–1260. Dubovi, E.J., 1992. Genetic diversity and BVD virus. Comp. Immunol., Microbiol. Infect. Dis. 15 (3), 155–162. Duffy, S., Shackelton, L.A., Holmes, E.C., 2008. Rates of evolutionary change in viruses: patterns and determinants. Nat. Rev. Genet. 9 (4), 267–276. Edgar, R.C., 2004. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32 (5), 1792–1797. Felsenstein, J., 1981. Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17 (6), 368–376. Fletcher, S.P., Jackson, R.J., 2002. Pestivirus internal ribosome entry site (IRES) structure and function: elements in the 5′ untranslated region important for IRES function. J. Virol. 76 (10), 5024–5033. Flores, E.F., Ridpath, J.F., Weiblen, R., Vogel, F.S.F., Gil, L.H.V.G., 2002. Phylogenetic analysis of Brazilian bovine viral diarrhea virus type 2 (BVDV-2) isolates: evidence for a subgenotype within BVDV-2. Virus Res. 87 (1), 51–60. Gates, M.C., Woolhouse, M.E., Gunn, G.J., Humphry, R.W., 2013. Relative associations of cattle movements, local spread, and biosecurity with bovine viral diarrhoea virus (BVDV) seropositivity in beef and dairy herds. Prev. Vet. Med. 112 (3–4), 285–295. Gilbert, S.A., Burton, K.M., Prins, S.E., Deregt, D., 1999. Typing of bovine viral diarrhea viruses directly from blood of persistently infected cattle by multiplex PCR. J. Clin. Microbiol. 37 (6), 2020–2023. Gonzalez, L.A., Schwartzkopf-Genswein, K.S., Bryan, M., Silasi, R., Brown, F., 2012. Benchmarking study of industry practices during commercial long haul transport of cattle in Alberta, Canada. J. Anim. Sci. 90 (10), 3606–3617. Greiser-Wilke, I., Haas, L., Dittmar, K., Liess, B., Moennig, V., 1993. RNA insertions and gene duplications in the nonstructural protein p125 region of pestivirus strains and isolates in vitro and in vivo. Virology 193 (2), 977–980. Harada, T., Tautz, N., Thiel, H.J., 2000. E2-p7 region of the bovine viral diarrhea virus polyprotein: processing and functional studies. J. Virol. 74 (20), 9498–9506. Hasegawa, M., Kishino, H., Yano, T.-a., 1985. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22 (2), 160–174. Iqbal, M., Flick-Smith, H., McCauley, J.W., 2000. Interactions of bovine viral diarrhoea virus glycoprotein E(rns) with cell surface glycosaminoglycans. J. Gen. Virol. 81 (Pt 2), 451–459. Isken, O., Langerwisch, U., Schonherr, R., Lamp, B., Schroder, K., Duden, R., et al., 2014. Functional characterization of bovine viral diarrhea virus nonstructural protein 5A by reverse genetic analysis and live cell imaging. J. Virol. 88 (1), 82–98. Jackova, A., Novackova, M., Pelletier, C., Audeval, C., Gueneau, E., Haffar, A., et al., 2008. The extended genetic diversity of BVDV-1: typing of BVDV isolates from France. Vet. Res. Commun. 32 (1), 7–11.
Author statements Funding Information University of Calgary, Faculty of Veterinary Medicine (10004757). University of Calgary, Margaret Gunn Endowment for Animal Research (10005041). Alberta Livestock and Meat Agency (2011F096R). Conflicts of interest The authors declare that they hold no conflicts of interest regarding this research. Ethical statement All viral isolates collected from animals in this study were collected in accordance with the approved animal care protocol VSACC AC110081/AC14-0144 at the University of Calgary. Acknowledgments We would like to acknowledge the invaluable contribution of samples from the following groups: John Robinson and Melissa Trapp at the Plant and Animal Health Centre, BC Ministry of Agriculture, Dale Godson at Prairie Diagnostic Services, Carl Gagnon at Université de Montreal, and Davor Ojkic at the University of Guelph. Additionally, we thank Sandeep Atwal and Gurleen Brar for their assistance with sample processing. 237
Virology 509 (2017) 232–238
A. Chernick, F. van der Meer
pestiviruses: determination of cleavage sites. J. Virol. 71 (7), 5415–5422. Uttenthal, A., Stadejek, T., Nylin, B., 2005. Genetic diversity of bovine viral diarrhoea viruses (BVDV) in Denmark during a 10-year eradication period. APMIS 113 (7–8), 536–541. Vilcek, S., Alenius, S., Paton, D.J., Mittelholzer, C., Belak, S., 1999. Genetic clustering of bovine viral diarrhoea viruses in cattle farms: genetic identification and analysis of viruses directly from cattle sera. Vet. J. 158 (1), 33–38. Vilcek, S., Durkovic, B., Kolesarova, M., Greiser-Wilke, I., Paton, D., 2004. Genetic diversity of international bovine viral diarrhoea virus (BVDV) isolates: identification of a new BVDV-1 genetic group. Vet. Res. 35 (5), 609–615. Vilček, Š., Paton, D.J., Durkovic, B., Strojny, L., Ibata, G., Moussa, A., et al., 2001. Bovine viral diarrhoea virus genotype 1 can be separated into at least eleven genetic groups. Arch. Virol. 146 (1), 99–115. Warrener, P., Collett, M.S., 1995. Pestivirus NS3 (p80) protein possesses RNA helicase activity. J. Virol. 69 (3), 1720–1726. Weiskircher, E., Aligo, J., Ning, G., Konan, K.V., 2009. Bovine viral diarrhea virus NS4B protein is an integral membrane protein associated with Golgi markers and rearranged host membranes. Virol. J. 6, 185. Wiskerchen, M., Belzer, S.K., Collett, M.S., 1991. Pestivirus gene expression: the first protein product of the bovine viral diarrhea virus large open reading frame, p20, possesses proteolytic activity. J. Virol. 65 (8), 4508–4514. Xia, H., Liu, L., Wahlberg, N., Baule, C., Belak, S., 2007. Molecular phylogenetic analysis of bovine viral diarrhoea virus: a Bayesian approach. Virus Res. 130 (1–2), 53–62. Xu, J., Mendez, E., Caron, P.R., Lin, C., Murcko, M.A., Collett, M.S., et al., 1997. Bovine viral diarrhea virus NS3 serine proteinase: polyprotein cleavage sites, cofactor requirements, and molecular model of an enzyme essential for pestivirus replication. J. Virol. 71 (7), 5312–5322. Xue, W., Mattick, D., Smith, L., Umbaugh, J., Trigo, E., 2010. Vaccination with a modified-live bovine viral diarrhea virus (BVDV) type 1a vaccine completely protected calves against challenge with BVDV type 1b strains. Vaccine 29 (1), 70–76. Yu, H., Grassmann, C.W., Behrens, S.E., 1999. Sequence and structural elements at the 3′ terminus of bovine viral diarrhea virus genomic RNA: functional role during RNA replication. J. Virol. 73 (5), 3638–3648. Yuan, M., Lu, T., Li, C., Lu, L., 2013. The evolutionary rates of HCV estimated with subtype 1a and 1b sequences over the ORF length and in different genomic regions. PloS One 8 (6), e64698. Zharkikh, A., 1994. Estimation of evolutionary distances between nucleotide sequences. J. Mol. Evol. 39 (3), 315–329. Zhong, W., Gutshall, L.L., Del Vecchio, A.M., 1998. Identification and characterization of an RNA-dependent RNA polymerase activity within the nonstructural protein 5B region of bovine viral diarrhea virus. J. Virol. 72 (11), 9365–9369.
Kearse, M., Moir, R., Wilson, A., Stones-Havas, S., Cheung, M., Sturrock, S., et al., 2012. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28 (12), 1647–1649. Kosakovsky Pond, S.L., Frost, S.D., Muse, S.V., 2005. HyPhy: hypothesis testing using phylogenies. Bioinformatics 21 (5), 676–679. Kosakovsky Pond, S.L., Posada, D., Gravenor, M.B., Woelk, C.H., Frost, S.D., 2006. Automated phylogenetic detection of recombination using a genetic algorithm. Mol. Biol. Evol. 23 (10), 1891–1901. Lanave, C., Preparata, G., Saccone, C., Serio, G., 1984. A new method for calculating evolutionary substitution rates. J. Mol. Evol. 20 (1), 86–93. Liu, L., Xia, H., Wahlberg, N., Belak, S., Baule, C., 2009. Phylogeny, classification and evolutionary insights into pestiviruses. Virology 385 (2), 351–357. Luzzago, C., Ebranati, E., Sassera, D., Lo Presti, A., Lauzi, S., Gabanelli, E., et al., 2012. Spatial and temporal reconstruction of bovine viral diarrhea virus genotype 1 dispersion in Italy. Infect., Genet. Evol.: J. Mol. Epidemiol. Evolut. Genet. Infect. Dis. 12 (2), 324–331. Paton, D.J., Lowings, J.P., Barrett, A.D., 1992. Epitope mapping of the gp53 envelope protein of bovine viral diarrhea virus. Virology 190 (2), 763–772. Ridpath, J.F., 2010. Bovine viral diarrhea virus: global status. Vet. Clin. North Am. Food Anim. Pract. 26 (1), 105–121. Ridpath, J.F., Bolin, S.R., 1998. Differentiation of types 1a, 1b and 2 bovine viral diarrhoea virus (BVDV) by PCR. Mol. Cell Probes. 12 (2), 101–106. Ridpath, J.F., Neill, J.D., Frey, M., Landgraf, J.G., 2000. Phylogenetic, antigenic and clinical characterization of type 2 BVDV from North America. Vet. Microbiol. 77 (1– 2), 145–155. Ridpath, J.F., Fulton, R.W., Kirkland, P.D., Neill, J.D., 2010. 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. 22 (2), 184–191. Ronecker, S., Zimmer, G., Herrler, G., Greiser-Wilke, I., Grummer, B., 2008. Formation of bovine viral diarrhea virus E1-E2 heterodimers is essential for virus entry and depends on charged residues in the transmembrane domains. J. Gen. Virol. 89 (Pt 9), 2114–2121. Rumenapf, T., Unger, G., Strauss, J.H., Thiel, H.J., 1993. Processing of the envelope glycoproteins of pestiviruses. J. Virol. 67 (6), 3288–3294. Suzuki, Y., Yamaguchi-Kabata, Y., Gojobori, T., 2000. Nucleotide substitution rates of HIV-1. AIDS Rev. 2, 39–47. Tamura, K., Nei, M., 1993. Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Mol. Biol. Evol. 10 (3), 512–526. Tautz, N., Elbers, K., Stoll, D., Meyers, G., Thiel, H.J., 1997. Serine protease of
238