Molecular Characteristics of Bovine Viral Diarrhea Virus Strains Isolated from Persistently Infected Cattle

Simple Summary Bovine viral diarrhea virus (BVDV) is a pathogen that is widespread throughout the world, causing serious economic losses. Most of the current BVDV infection research focuses on genotyping based on the 5′-UTR sequence, neglecting the molecular characteristics of the whole genome of BVDV strains. Cows infected with non-cytopathic (NCP) BVDV in early pregnancy may give birth to persistently infected (PI) cattle, which are the main sources of virus transmission in the herd. The present study isolated and identified nine BVDV strains from PI cattle in China. The complete genome of the new isolates was sequenced for phylogenetic analysis, recombination analysis and sequence analysis in 5′-UTR, E0 and E2. Abstract In this study, we reported the isolation, identification, and molecular characteristics of nine BVDV strains that were isolated from the serum of persistently infected cattle. The new strains were designated as BVDV TJ2101, TJ2102, TJ2103, TJ2104, TJ2105, TJ2106, TJ2107, TJ2108 and TJ2109. The TJ2102 and TJ2104 strains were found to be cytopathic BVDV, and the other strains were non-cytopathic BVDV. An alignment and phylogenetic analysis showed that the new isolates share 92.2–96.3% homology with the CP7 strain and, thus, were classified as the BVDV-1b subgenotype. A recombination analysis of the genome sequences showed that the new strains could be recombined by the major parent BVDV-1a NADL strain and the minor parent BVDV-1m SD-15 strain. Some genome variations or unique amino acid mutations were found in 5′-UTR, E0 and E2 of these new isolates. In addition, a potential linear B cell epitopes prediction showed that the potential linear B cell epitope at positions 56–61 is highly variable in BVDV-1b. In conclusion, the present study has identified nine strains of BVDV from persistently infected cattle in China. Further studies on the virulence and pathogenesis of these new strains are recommended.


Introduction
Bovine viral diarrhea virus (BVDV), the pathogen that causes bovine viral diarrhea mucosal disease (BVD-MD), is a member of the Pestivirus genus in the Flavivirus family [1]. viral shedding from excrement and secretions [33]. In this study, we successfully isolated nine BVDV strains from the serum of calves with diarrhea. To understand the genetic characteristics and possible evolutionary origins of circulating BVDV strains, we amplified and analyzed their whole genome sequences, and reported the molecular characteristics of the new isolates. These findings could contribute to the prevention and control of BVDV.

Cell Culture and Virus Isolation
Madin-Darby Bovine Kidney (MDBK) cells were grown in DMEM (Gibco, New York, NY, USA) containing 10% BVDV-free fetal bovine serum (Royacel, Lanzhou, China). The serum samples were obtained from PI cattle with mucosal disease at a field in Tianjin in 2021. The PI cattle were tested to be positive for the BVDV antigen and negative for antibodies. For virus isolation, the serum samples were diluted one-fold with 10 mM phosphate buffered saline (PBS), then were filtered with a 0.22 µm filter before virus infection. Monolayer MDBK cells were infected with the treated samples in 6-well culture plates. The cytopathic conditions were observed and recorded every day. After 96 h, the infected cells were lysed using three cycles of freeze-thaw, and centrifuged at 5000 rpm for 3 min. Then, the supernatant was taken for the next infection. The 5th passages of viruses were used for further research.

Indirect Immunofluorescence Assay (IFA)
To identify the isolated BVDV strains, monolayer MDBK cells were infected with the 5th passages of the isolated viruses. Mock-infected cells were used as the negative controls. After culturing for 48 h, the cells were fixed with 4% paraformaldehyde, incubated with 0.2% TritonX-100, and blocked with 1% bovine serum albumin. Then, the cells were incubated with BVDV monoclonal antibodies (VMRD, Washington, WA, USA) and Goat anti-Rabbit IgG/FITC (Solarbio ® , Beijing, China) before examination using the fluorescence microscope.

RNA Extraction and RT-PCR
Total RNA was extracted from the lysate of infected cells using a TransZol Up Plus RNA Kit (TransGen, Beijing, China) according to the manufacturer's instructions. Briefly, the samples were lysed using a triple volume of TransZol reagent, mixed with 0.2 volume of an RNA Extraction Agent, and shaken vigorously for 5 min. After centrifugation at 10,000× g, 4 • C for 15 min, the upper aqueous phase was washed with anhydrous ethanol. Then, the RNA was resuspended in 30 µL Diethyl pyrocarbonate (DEPC) treated water. Reverse transcriptase reactions were carried out using EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen, Beijing, China) in a 20 µL reaction volume containing 6 µL of extracted RNA, 1 µL Anchored Oliogo(dt)18 Primer, 10 µL of 2 × ES Reaction Mix, 1 µL of EasyScript RT/R1 Enzyme Mix, 1 µL of gDNA Remover, and 1 µL of RNase-free Water. The following reaction was 42 • C for 30 min and 85 • C for 5 s. PCR amplification was conducted using the 2 × Pfu PCR Mix (Tiangen, Beijing, China) in a 20 µL reaction volume containing 10 µL of 2 × Pfu PCR Mix, 1 µL of upstream and downstream primers (5 -TCAGCGAAGGCCGAAAAGAGG-3 , 5 -TCCATGTGCCATGTACAGCAGAG-3 ), 2 µL of cDNA template, and 6 µL ddH 2 O. The following reaction was 94 • C for 5 min, followed by 35 cycles at 94 • C for 30 s, 46 • C for 30 s, and 72 • C for 1 min, and an extension at 72 • C for 10 min.

Complete Genomic Amplification and Sequencing
Twelve pairs of primers (Table A1) were designed using the Primer Premier 5 according to the BVDV JL-1 strain (GenBank No. KF501393.1) and BVDV CC13B strain (GenBank No. KF772785.1) and were used to amplify the BVDV complete genomic sequence. The amplified products were cloned into a pMD-18T vector using pMD™ 18-T Vector Cloning Kit (TaKaRa, Beijing, China). Briefly, 4 µL of the purified amplification product was mixed with 1 µL of pMD18-T Vector and 5 µL of Solution I and incubated at 16 • C for 12 h. After the recombinant vector was transformed into DH5α competent cells (TaKaRa, Beijing, China), the single colony with AMP resistance was picked for amplification, and the plasmid was extracted using a TIANprep Mini Plasmid Kit (Tiangen, China) and submitted to Sangon Biotech Company (Shanghai, China) for DNA sequencing. To obtain the complete genomic sequence, the fragment sequences were trimmed and assembled using the SeqMan program in DNASTAR software.

Complete Genomic Sequence Analysis
He sequence alignment for the isolated BVDV, with the known BVDV strains listed in Table A2, was performed using Clustal W software. Phylogenetic analysis was carried out by the neighbor-joining method using MEGA7 software. The E2 glycosylation sites were analyzed using the NetNGlyc 1.0 Server (http://www.cbs.dtu.dk/services/NetNGlyc, accessed on 7 November 2022). The potential liner B cell epitopes on the E2 protein were predicted using the Bepipred Linear Epitope Prediction 2.0 method in IEDB Analysis Resource (http://tools.iedb.org/bcell/, accessed on 7 November 2022). Tertiary structures on the E2 protein of the new isolates were predicted using SWISS-MODEL (https:// swissmodel.expasy.org/, accessed on 28 April 2023) and the images were generated using Pymol software.

Recombination Analysis
The recombination analysis between the new isolates and the reference strains were carried out using RDP4 software and SimPlot software. Seven methods, including RDP, GENECONV, Chimaera, MaxChi, Bootscan, SiScan, and 3Seq were used in RDP4. Fragments with recombination signals detected by at least four methods are considered to take place as recombination events.

Virus Isolation and Identification
Nine treated serum samples were cultured and passaged in MDBK cells. The presence of BVDV in the MDBK cells was confirmed with RT-PCR and an immunofluorescence assay. The amplification with 255 bp ( Figure 1A) and the specific green signals ( Figure 1B,C) were detected in all infected cells, suggesting the existence of BVDV. After five passages of the virus, nine isolates were obtained and named BVDV TJ2101, TJ2102, TJ2103, TJ2104, TJ2105, TJ2106, TJ2107, TJ2108, and TJ2109, respectively. Cytopathic effects were only observed in the TJ2102 and TJ2104 infected plates ( Figure 1E), suggesting the BVDV TJ2102 and TJ2104 strains were CP BVDV strains, and the other isolates were NCP BVDV strains.

Genome Sequencing and Phylogenetic Analysis
The complete genomic sequence of nine isolates were assembled and submitted to GenBank with accession numbers OR004805 to OR004813. The complete genomic sequence of the 9 new isolates were compared with 31 reference strains. Results showed that the genomic nucleotide identity between the new isolates and other BVDV strains was 78.5-92.6%, and the new isolates were most similar to the CP7 strain, sharing the 92.2-92.6% homology with it ( Figure A1). As shown in Figures 2 and 3, the new isolates and the CP7 and CC13B strains were in the same branch, suggesting that the new isolates could be classified as a BVDV-1b subgenotype.

Sequence Analysis of 5 -UTR
As shown in Figure 4, compared with the isolates of other subgenotypes, several molecular characteristics were found in the new isolates, including amino acid T at position 86 (T 86 ) and A 137 , which were consistent with other 1b strains. In addition, several unique mutations, including nucleotide T82C in the 5 -UTR of the new isolates, BVDV-1d and BVDV-1k, and nucleotide G302A in the BVDV-1b and BVDV-1c strains, were found.
Compared with the BVDV-1d strain (10JJ-SKR), TJ2102, TJ2104 and TJ2108 showed more continuous nucleotide deletions than other strains at positions 50-57. In addition, compared with the BVDV-1c and BVDV-1h strains, two nucleotide deletions were found in all of the BVDV-1b strains at positions 223-224, which had one more deletion than the strains of other subgenotypes. In addition, a nucleotide deletion at position 309 was observed in most BVDV-1b strains, which was not found in other strains.

Amino Acid Analysis of E2
As shown in Figure 5A, some molecular characteristics, including T 24 , T 173 , V 198 , A 199 , S 252 and M 371 , were detected in the E2 protein of the new isolates, which was consistent with other BVDV-1b strains. In addition, I 27 , L 84 , D 167 , A 194 , Q 236 , and L 254 were also found in the new isolates, which was similar to several BVDV-1b strains. Furthermore, the new isolates had several unique amino acid sequence characteristics, including A 43 , EA 49,50 , V 182 and S 251 , which were different from other strains ( Figure 5C). Compared with other BVDV-1b strains, CP7, CC13B and HJ-1, the amino acid replacements E45V and E52K were observed in the new strains. The potential linear B cell epitopes on the E2 protein of the new isolates and three BVDV-1b strains were predicted using BepiPred-2.0. As shown in Figure 5B, amino acid sequence variations affected 6 out of 10 liner B cell potential epitopes of the new isolate and the E2-2 epitope was found to be the most variant region, followed by E2-6, E2-8, and E2-9. The conserved regions were observed in the E2-1, E2-3, E2-4, and E2-10 epitopes. In addition, the glycosylation sites were conserved in the new isolates except 230 NET, which was consistent with other BVDV-1b strains ( Figure 5D).

Genome Sequencing and Phylogenetic Analysis
The complete genomic sequence of nine isolates were assembled and submitted to GenBank with accession numbers OR004805 to OR004813. The complete genomic sequence of the 9 new isolates were compared with 31 reference strains. Results showed that the genomic nucleotide identity between the new isolates and other BVDV strains was 78.5-92.6%, and the new isolates were most similar to the CP7 strain, sharing the 92.2-92.6% homology with it ( Figure A1). As shown in Figures 2 and 3, the new isolates and the CP7 and CC13B strains were in the same branch, suggesting that the new isolates could be classified as a BVDV-1b subgenotype.

Sequence Analysis of 5′-UTR
As shown in Figure 4, compared with the isolates of other subgenotype lecular characteristics were found in the new isolates, including amino acid 86 (T 86 ) and A 137 , which were consistent with other 1b strains. In addition, s

Amino Acid Analysis of E2
As shown in Figure 5A, some molecular characteristics, including T 24 , T 173 , V 198 , A 199 , S 252 and M 371 , were detected in the E2 protein of the new isolates, which was consistent with other BVDV-1b strains. In addition, I 27 , L 84 , D 167 , A 194 , Q 236 , and L 254 were also found in the new isolates, which was similar to several BVDV-1b strains. Furthermore, the new isolates had several unique amino acid sequence characteristics, including A 43 , EA 49,50 , V 182 and S 251 , which were different from other strains ( Figure 5C). Compared with other BVDV-1b strains, CP7, CC13B and HJ-1, the amino acid replacements E45V and E52K were observed in the new strains. The potential linear B cell epitopes on the E2 protein of the new isolates and three BVDV-1b strains were predicted using BepiPred-2.0. As shown in Figure 5B, amino acid sequence variations affected 6 out of 10 liner B cell potential epitopes of the new isolate and the E2-2 epitope was found to be the most variant region, followed by E2-6, E2-8, and E2-9. The conserved regions were observed in the E2-1, E2-3, E2-4 ,and E2-10 epitopes. In addition, the glycosylation sites were conserved in the new isolates except 230 NET, which was consistent with other BVDV-1b strains ( Figure 5D).

Amino Acid Analysis of E0
E0 amino acid sequences of the new isolates were compared with strains of other BVDV subgenotypes. Seven linear epitopes exist in the BVDV E0 protein, including the epitopes 31GIWPEKIC38, 65NYTCCKLQ72, 127QARNRPTT134, 145SFAGTVIE152, 161VEDILY166, 114CRYDKNTDVNV124 and 116YDKNTDVNV124 [21]. As shown in Figure 6, all of the linear epitopes were found in the new isolates, while a unique amino acid replacement, a hydrophilic amino acid T mutating into a hydrophobic amino acid A, was only found at position 39 closed to the E0 epitope 31GIWPEKIC38 of the new isolates. The residues W 33 , L 71 , Q 127 , N 130 , S 145 , G 148 and T 102 were conserved in all BVDV-1 strains including the new isolates. However, an amino acid replacement at position 107 was only observed in the new isolates. In addition, several unique replacements, including the residues Q22L and P107L in the E0 of the new isolates, and R194M in the TJ2102 strain, were found. In addition, several amino acid sequence characteristics were observed, including S173 in the new isolates and BVDV-1a, I175 in BVDV-1b, and V161 in the BVDV-1b and BVDV-1o strains.
Vet. Sci. 2023, 10, x FOR PEER REVIEW 11 of 18 Figure 6. Amino acid sequence alignment of E0 of the 9 new BVDV isolates and 13 reference strains. Seven linear epitopes of E0 were found and are less lightly shadowed than the amino acid mutations region. * represent a number. For example, the * between 60 and 80 is 70.

Recombination Analysis
A recombination analysis of the full-length genome sequences of the new isolates and reference strains was conducted using RDP4 and SimPlot software. As shown in Table 1, recombination signals were found in most new isolates except TJ2103 around the 691- Figure 6. Amino acid sequence alignment of E0 of the 9 new BVDV isolates and 13 reference strains. Seven linear epitopes of E0 were found and are less lightly shadowed than the amino acid mutations region. * represent a number. For example, the * between 60 and 80 is 70.

Recombination Analysis
A recombination analysis of the full-length genome sequences of the new isolates and reference strains was conducted using RDP4 and SimPlot software. As shown in Table 1, recombination signals were found in most new isolates except TJ2103 around the 691-1155 nucleotides region, suggesting most new isolates may be recombined by the major parent BVDV-1a NADL strain and the minor parent BVDV-1m SD-15 strain ( Figure 7A,B).

Discussion
BVD is a significant bovine disease causing severe economic losses to the cattle industry. With the improvement of people's living standards, the demand for beef and milk is increasing gradually. Multinational cattle transportation and commodity trading have also promoted the spread and mutation of the virus [5]. In recent years, the research on BVDV has been increasing, but most research is focused on the genotyping and epidemiological investigation of BVDV, while studies on the molecular characteristics and genome sequence of BVDV have been neglected. PI cattle will emit a large amount of virus during their lifetime, which seriously threatens the health of herds and is the main reason for the widely spread BVDV [34]. In this study, nine BVDV strains were isolated from the serum of PI cattle, and their full-length genome sequences and molecular characteristics were analyzed in order to provide a potential genome variation mode of the virus isolated from PI cattle.
The TJ2102 and TJ2104 strains were CP BVDV strains, and the other isolates were NCP BVDV strains. Surprisingly, no nucleotide insertions were found in the NS2-3 protein of the two CP BVDV isolates. Studies have shown that the insertion between NS2 and NS3 is necessary for its cleavage in CP BVDV strains [35], whether this characteristic of the TJ2102 and TJ2104 strains will affect the expression of NS3 remains to be further studied. In this study, the genome sequences of the nine isolates were compared with other reference strains. The new isolates had the highest similarity with the CP7 strain, sharing 92.2-96.3%, higher than other subgenotypes of BVDV such as BVDV-1o IS26/01ncp (77.9-78.2%) and BVDV-1n Shitara/02/06 (78.5-78.8%). The phylogenetic tree based on the genome sequence comparison showed that the new isolates were in the same branch with the BVDV-1b strains (CP7, JL-1 and CC13B), which was consistent with the analysis based on the 5 -UTR and N pro sequences.
The 5 -UTR is a conserved region in pestiviruses that is often used to identify virus genotypes. In this study, some molecular characteristics (T 86 and A 137 ) were found in the 5 -UTR of every BVDV-1b strains. This finding could probably contribute to the identification of BVDV-1b strains. There exists an internal ribosomal entry site (IRES) in the secondary hairpin structure formed by the BVDV 5 -UTR, which is crucial for virus replication; changes in the secondary hairpin structure could directly affect the infectivity of the virus [36,37]. It has been found that deleting the 32-65 nucleotides or inserting 4 nucleotides at 63-66 in the stem-loop Ib did not affect the efficiency of the IRES, indicating that Ib was dispensable for BVDV translation. Deletions or insertions in any of the other four stem-loop regions, including II, IIIa, IIIc, and IIId, caused sharp decreases in IRES activity [38]. Therefore, deletions of the new isolates at the nucleotide positions 50-57 might not influence the IRES's activity; although, whether the base mutations in other stem-loops of the new isolates affects BVDV translation remains to be confirmed experimentally.
Genome recombination plays an important role in viral genome diversification and genetic variation and affects pathogenicity and environmental adaptability. In this study, recombination events were identified using RDP4 and SimPlot software and were verified by at least four methods. Recombination signals were found in most of the new isolates except the TJ2103 strain. The major and minor parents of the recombinant strains were found to be BVDV-1a (the NADL strain isolated in USA, in 1988) and BVDV-1m (the SD-15 strain isolated from cattle in northern China, in 2015). As the center of the cattle industry in China, most of the cattle and dairy products were transported from Inner Mongolia to the eastern and central regions where the products were in high demand. Studies have shown that BVDV-1m might be the dominant subgenotype spreading in Inner Mongolia [39]. Therefore, the BVDV-1m strains might have been spread from Inner Mongolia to other regions. Moreover, in the past few decades, China has imported a large amount of Holstein semen and embryos from North America and Europe [40]. Therefore, we speculate that frequent commodity trade and transportation might directly or indirectly lead to the gene recombination and mutation of BVDV.
E0 is a relatively conserved structural protein in BVDV, containing seven linear epitopes, which can induce the production of neutralizing antibodies. Prediction based on hydrophilicity is a classic method used for predicting potential B cell epitopes; so changes in the hydrophilicity of amino acids may affect epitope predictions [41]. The amino acid at position 39 of the new isolates' E0 protein had an opposite change in hydrophilicity; whether this change could affect the potential B cells' epitopes remains to be further verified. For CSFV, the residues W 33 , L 71 , Q 127 , N 130 , S 145 , G 148 , T 102 and D 107 may be crucial for interactions with antibodies [42] In the present study, one amino acid replacement, P107L, has been found in the E0 protein of the new isolates, which could influence the effectiveness of vaccination. In addition, binding to the cell surface glycosaminoglycans (GAGs) and attachment to the cell surface of E0 plays an important role in the viral infection of host cells [22,23]. The cluster of residues 480KKLENKSK487 localized at the C-terminal end of E0 was proved to be a binding site of GAGs. Replacements at the amino acids 481 and 485 in E0 resulted in the loss of the ability to bind to cells and even the failure to infect cells [24]. The new isolates were found to be conserved in this region.
The E2 protein is an envelope glycoprotein in BVDV that is necessary for virus recognition, adsorption, and the entry into host cells. After the C-terminus of E0 binds to cell surface glycosaminoglycans (GAGs), E2 binds to CD46 to allow the virus to enter the cell [43,44]. E2 is the most antigenic protein in BVDV, and it is the main site that induces the production of neutralizing antibodies [45]. In the present study, some unique molecular features were found in the new isolates. Whether the epitopes are arranged in the correct conformation determines the efficiency of antibody binding. A change of a certain amino acid residue may affect the sequence or three-dimensional structure of the epitope [46]. These amino acid variations altered eight B cell linear epitopes of the isolate's E2 protein, especially around position 60. Linear epitopes of viral proteins affect the antigenic structure and virus-antibody interactions, and the changes may affect the effectiveness of vaccination. Some special epitopes can be used in the development of new vaccines but the antigenic epitopes predicted using software cannot fully reflect the actual site of virus-induced neutralizing antibody production [47]. N-glycosylation is the most common form of protein glycosylation and has two forms, Asn-X-Ser and Asn-X-Thr. In the present study, a NET-to-NES substitution at position 230 was found in E2 of the new isolates and in other BVDV-1b strains. Studies have found that the replacement of serine with threonine greatly reduces the efficiency of N-glycosylation [48]. Whether this change affects the structure and function of the E2 protein remains to be further verified.
Consistent with previous findings, the BVDV genome isolated from PI cattle had more substitutions, and these amino acid substitutions mainly occurred in the structural proteins E0 and E2 [49]. In the present study, some molecular characteristics of PI bovine isolates were found; whether these features are conserved in PI bovine isolates of other subgenotypes needs further study.

Conclusions
In this study, nine isolates were isolated from PI calves and identified as the BVDV-1b subgenotype. Some unique variations that may be related to viral replication and immune response were found in the new isolates. A novel B cell potential linear epitope was found in the E2 protein of the new isolates, suggesting that the new isolate may be a potential candidate for vaccine optimization. Recombination analysis suggested that commodity trade across regions may have contributed to the emergence of viral diversity. Whether these unique molecular characteristics are conserved in the PI isolates of other subgenotypes remains to be further verified. This study provides a potential genetic variation model of PI isolates in order to provide a reference for the genetic evolution of BVDV.  Institutional Review Board Statement: Serum was sampled from the jugular vein and no further procedures were carried out on the study cows. The sample collection was approved by the relevant farms and the current study was approved by the Animal Care and Use Committee of Hebei Agricultural University (Protocol number 2020044).
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The complete genome sequences of the nine BVDV isolates (Accession number: OR004805-OR004813) have been submitted to GenBank.

Conflicts of Interest:
The authors declare no conflict of interest.  Figure A1. The genomic nucleotide identity between the new isolates and reference BVDV strains. The numbers from 1 to 27 represent the strains from TJ2109 to SD0803. The black cells represent the dividing line between percent identity and distance.