Virus isolation and identification
Following inoculation of MDBK cells with three blind serial passages of the virus, two of the samples were detected as BVDV positive by RT-PCR. The IF assay revealed that MDBK cells inoculated with the third passaged cell culture were BVDV positive (Fig. 1). Virus isolated from serum sample collected in 2013 was named as BVDV BJ–2013 and virus isolated from commercial fetal bovine serum was named as BVDV BJ–2016. The significant cytopathic effects (CPE) were observed in MDBK cells infected with isolate BVDV BJ–2013 (Fig. S1). There were no CPEs visible in MDBK cells infected with isolate BVDV BJ–2016. The negatively stained virus particles collected from the third passaged cell cultures were approximately 40 nm–60 nm in length, and displayed a typical BVDV morphology (Fig. S2).
Genome analysis
The two complete genomes of BVDV isolates were sequenced. The two sequences, herein as BVDV BJ–2013 (MH490942) and BVDV BJ–2016 (MH490943), were identified as having 99% homology with the nucleotide identity of various Pestivirus A strains in GenBank. The BVDV BJ–2013 isolate had a genome of 12,305 nt in length and contained a single open reading frame (ORF) that encoded a polyprotein of 3,910 amino acids starting at position 385 and ending at position 12,078. The genome was flanked by 5’ and 3’ UTRs located at positions 1–384 and 12,079–12,305 nt, respectively. Isolate BVDV BJ–2016 was a full length genome of 12,273 nt in length that encoded a polyprotein of 3898 amino acids starting at position 384 and finishing at position 12077. The genome was flanked by 5’ and 3’ UTRs, located at positions 1–383 and 11972–12,273, respectively.
Phylogenetic analysis
The maximum likelihood Phylogenetic tree inferred from the complete genomes showed that isolate BVDV BJ–2016 (identified by a white arrow in Fig. 2) was located within the clade composed of Pestvirus A strains of genotype 1B (identified by a grey rectangle in Fig. 2). This clade was the largest in the genome phylogenetic tree, which was consistent with existing data showing that, in many countries, BVDV of genotype 1B was the most prevalent subtype[4]. Further, isolate BVDV BJ–2016 was closely related to isolates from South Korea (JQ418634) and the US (JN380088), thus corroborating predictions about the ability of subtype 1B to spread worldwide. Isolate BVDV BJ–2013 (identified by a black arrow in Fig. 2) was clustered within a group composed of Pestivirus A strains of genotype 1A. Isolate BVDV BJ–2013 was related to isolates from the USA (AF091605) and Germany (AF041040). This subtype was reportedly the second most prevalent in many countries, including in Asian countries [4]. So, the two BVDV isolates in this study, BVDV BJ–2013 and BVDV BJ–2016, were the members of Pestivirus A.
Recombination analysis of the polyprotein product of Pestivirus A
Initially, we used default parameters to screen for a recombination signal in the alignment of the 48 Pestivirus Astrains. This screening process correctly identified strain KJ541471 (isolated in China in 2013) that had previously been detected by the phylogenetic analysis (indicated by a black triangle in Fig. 2). During screening, certain parameters were modified, such as the probability threshold (P-value cut-offs), windows and step sizes, in order to avoid false positive signals of recombination. A strong (based on P-values) signal of recombination was only detected when strain KJ541471 was included in the alignment. On the contrary, alignments containing all strains or exclusively genotype 1A or genotype 1B sequences did not indicate recombination, even using the upper threshold of P < 0.05. By doing this we were able to determine that KJ541471 was a mosaic of sequences derived from two isolates of the genotype 1A. To determine the breakpoints of the recombinant strain, an alignment containing only genotype 1A strains was used. For this step we selected the method that provided the highest probability of recombination (the maxchi method) to locate the breakpoints (Fig. 3b). The results indicated that strain KJ541471 had two breakpoints (6980/6981 and 9565/9566) in the polyprotein region. Subsequently, the polyprotein alignment was partitioned into two fragments: the first was a concatenated alignment comprising nucleotides 1 to 6980 plus nucleotides 9567 to 11682. The second fragment corresponded to nucleotides 6981 to 9566 (Fig. 3b). Maximum likelihood trees were constructed for each fragment and the results were in agreement with the predicted locations of the breakpoints because they showed that the strain KJ541471 was located in conflicting positions on the trees. In particular, the tree constructed from the concatenated fragment (Fig. 3a) indicated that strain KJ541471 was related to strain KP941584 that was isolated in 2014 from the US. On the other hand, the tree constructed with fragment 6981 to 9566 (Fig. 3c) indicated that strain KJ541471 and isolate BJ–2013 are closely related. Apart from strain KJ541471, mosaic genomes were also detected in other strains (indicated by black rectangles in Fig. 2), but these were not regarded as recombinants due to their weak signals and unclear locations of the breakpoints. These problematic strains were excluded from the analyses of evolutionary rates and the coalescent analysis.
Evolutionary rates of Pestivirus A genotype1A and 1B
In order to explore the differences between genotype 1A and 1B we estimated their genetic distances based on pairwise distances. In particular, isolate BVDV BJ–2016 was compared to South Korean isolate JQ418634 (indicated by a white dot in Fig. 2). The nucleotide distance between these two isolates was, on average, 0.052 ± 0.003. The nucleotide distances among the isolates of genotype 1B were, on average, 0.048 ± 0.001. The same measurements were performed with strains of clade 1A. The nucleotide distance between isolate BVDV BJ–2013 and isolate AF091605 from Oregon in the US (indicated by a black circle in Fig. 2) was, on average 0.072 ± 0.002 while the overall distance among all 1A strains was 0.10 ± 0.021. The higher mean distances observed in genotype 1A might be related to the reduced sample size. For this reason, additional analyses were performed. Two maximum likelihood methods were used to estimate the evolutionary rates in each codon of the polyprotein region of BVDV genomes. The results indicated that the relative rates (Fig. 4a) and the dN/dS ratios (Fig. 4b) were both similar between subtypes (based on pairwise comparisons using the Kruskall Wallis test).
Comparison analyses and phylogenetic analyses of genes
To explore the differences between BVDV BJ–2013/BJ–2016 and other BVDV isolates, comparison analyses were performed based on nucleotide and amino acid sequences. The nucleotide identity of sequenced genomes was 80% when compared to each other. BVDV BJ–2013 showed highest nucleotide identity with BVDV 08GB44–1 (JQ418633.1) (99.9%)which isolated from South Korea. BVDV BJ–2016 shared highest nucleotide identity (96.4%) with BVDV HJ–1 from China. The amino acid similarity analysis of the two isolates showed the BVDV BJ–2013 shared highest amino acid similarity BVDV 08GB44–1 and BVDV BJ–2016 showed highest amino acid identity with BVDV HJ–01. Interestingly, there were 12 amino acid (EGNWLVNADRLM) insertion in non-structural protein 2 (NS2) of BVDV BJ–2013 compared with BVDV BJ–2016 (Fig. 5). It was similar with the insertion observed in other BVDV strains of cytopathogenic biotype BVDV, such as VE138cp05(LT907991), CP1741-A(AF220247), NADL(AJ133738),Vedevac(AJ585412).
The maximum likelihood phylogenetic trees (Fig. S3-Fig. S8) inferred from nucleotide sequences of genes Npro, Erns, E2, NS2, NS3, NS5B showed that isolate BVDV BJ–2013 was closely related with BVDV 08GB44–1(South Korea) within the clade composed of Pestivirus A strains of genotype 1a. Isolate BVDV BJ–2016 was closely related to BVDV 08GB45–2 (South Korea), PJ(USA), Powder(USA), HJ–1 (China). These results were consistent with that of phylogenetic analysis based on polyprotein nucleotide sequence.