Genetic recombination of bovine viral diarrhea virus subgenotype ‐1a and ‐1c in persistently infected dairy cattle

The bovine viral diarrhea virus (BVDV) is a major viral pathogen in cattle worldwide. In Indonesia, diversity in subgenotypes of BVDV‐1 has been observed, with the highest proportion of subgenotype ‐1a, followed by ‐1c, ‐1b, and ‐1d. So far, phylogenetic analysis of BVDV‐1 is based on nucleotide sequences of the 5′ UTR and partial NS5B regions. Accuracy in identifying the subgenotype and antigenic type is critical for vaccine development and effective vaccination. The aim of this study was to determine genetic recombination of BVDV through phylogenetic analysis of five different regions (5′ UTR, NPro, E2, NS3, and NS5B) of BVDV in persistently infected dairy cattle. Five isolates were sequenced using next‐generation sequencing, and data were analyzed with the CLC Genomic Workbench 9.0 and MEGA‐X programs. Phylogenetic analysis based on the 5′ UTR (275 nt), NPro (504 nt), E2 (1,122 nt), NS3 (2,049 nt), and NS5B (2,157 nt) regions indicated that one BVDV isolate from Banyumas, Central Java, could be classified into different subgenotypes based on the E2 region (‐1c), but the same subgenotype based on the other four regions (‐1a), suggesting the presence of genetic recombination of the BVDV subgenotypes ‐1a and ‐1c in persistently infected dairy cattle.


Introduction
Bovine viral diarrhea virus (BVDV) is an important vi ral pathogen in cattle that has spread throughout the world and has an economic impact with respect to animal hus bandry (Houe et al. 1995; Deregt et al. 2005. BVDV was first reported in 1946 in North America and has now ex panded throughout the Americas and into Europe, Aus tralia, Africa, and Asia, including Indonesia (Wiyono et al. 1989; Brodersen 2014. The incidence of BVDV in In donesia was first described in Balinese cattle in South Su lawesi in 1989 (Wiyono et al. 1989).
Persistent infection (PI) in cattle occurs when the non cytopathic BVDV infects pregnant cattle during the first trimester (40-120 d) of gestation, when the fetal immune system has not completely developed (Neill 2013). PI in cattle is determined by collecting blood samples twice at intervals of at least 3-4 weeks, after which PI is de tected using antigens (Firat et al. 2002). Persistently in fected cattle are a potential source of transmitting the virus and therefore, the disease, to herd communities (Liebner Tenorio 2005).
BVDV is a member of the genus Pestivirus, and the family Flaviviridae. It is a small, approximately 12.3 kblong, enveloped, singlestranded positivesense RNA virus (Chernick and van der Meer 2017). The genome of BVDV contains one openreading frame that encodes a large polyprotein with the sequence NProCE1E2p7 NS2NS3NS4aNS4bNS5aNS5b (Neill 2013). BVDV has two genotypes: BVDV1 and BVDV2 (Nagai et al. 2004; Vilcek et al. 2005. Genotype classification of BVDV is usually based on comparison of sequences from three genetic regions: 5′ UTR, Npro, and E2. Initially, there were only two subgenotypes of BVDV1 isolates identified, namely, BVDV1a and BVDV1b (Vilcek et al. 2005). However, based on different regions (5′ UTR, Npro, E2), BVDV1 is now divided into 21 subgenotypes (BVDV1a 1u) (Gao et al. 2014; Yeşilbağ et al. 2014; Baz zucchi et al. 2017); whereas BVDV2 has four subgeno types, probably owing to the lack of analyzed virus col lections (Vilcek et al. 2005; Giangaspero et al. 2008. The distribution of BVDV1 in Indonesia is indicated by a number of samples from Java that have been evaluated based on three genes: NS5B, NPro, and 5′ UTR (Saepulloh et al. 2015; Karimy 2016; Wuryastuti et al. 2018. Mean while, BVDV2 has not been found in individual samples of BVDVPI from Central Java. It is believed that there are two main BVDV1 subgenotypes, 1a and 1c (Saepulloh et al. 2015; Wuryastuti et al. 2015, in addition to other subgenotypes, 1b and 1d (Saepulloh et al. 2015). The distribution of three subgenotypes, namely, BVDV 1a, 1b, and 1c, in Central and East Java, Indonesia, re flects the genetic variability of BVDV1 viruses based on partial NS5B genes (Irianingsih et al. 2019). Phylogenetic analysis from different regions, specifically 5′ NCR, Npro, E2, NS3, and NS5B3' NCR, of BVDV isolates has been carried out in Japan (Nagai et al. 2004) and 5'UTR nd NPro in the UK (Booth et al. 2013). The determination of BVDV genotypes and subgenotypes is very useful for the correct classification of BVDV for the purposes of di agnosis, identification, and characterization of viruses, as well as for evaluating vaccine efficacy and molecular epi demiology (Nagai et al. 2004; Vilcek et al. 2005. This study seeks to determine genetic recombination in BVDV through phylogenetic analysis based on five different re gions, namely, noncoding regions (NCR): 5′ UTR, non structural genes (NPro, NS3, and NS5B), and structural genes (E2) of BVDVs from persistently infected dairy cat tle.

Virus isolation
Five BVDV isolates were obtained from serum samples of persistently infected dairy cattle in Central and East Java. These were collected by active and passive surveillance by the Disease Investigation Centre Wates, over the course of 20132016. Serum samples were inoculated in MDBK cells that were incubated at 37°C in 5% CO 2 in a main tenance medium (Minimal Essential Medium, MEM, sup plemented with 2% bovine serum, penicillinstreptomycin 100 IU/mL-100 µg/mL, gentamycin 50 µg/mL, fungizone 2.5 µg/mL, and HEPES buffer 0.01 M) for 3-5 d. Virus isolation required four passages in cell culture and then identification using realtime RTPCR on BVDV. Virus isolates were stored in a freezer (80°C) until subsequent experiments.

Targeted sequencing
RNA from BVDV isolates was extracted using the QI Aamp Viral RNA kit (Qiagen, Germany). Using RNA as a template, singlestranded cDNAs were generated with the SuperScript III FirstStrand Synthesis System for RT PCR (Invitrogen, California, USA) according to the man ufacturer's instructions. The BVDV genome was ampli fied using a series of four primers that overlapped frag ments of the proteincoding region (Chernick and van der Meer 2017) employing Platinum Taq DNA Polymerase High Fidelity (Invitrogen, California, USA). The ampli cons were gelpurified using the DNA Clean & Con centrator (Zymo Research, California, USA) and quan tified with the Qubit 2.0 and Qubit ds DNA HS Assay kits (Thermo Fisher Scientific, Waltham, US). Amplicons were pooled into equimolar amounts and adjusted for am plicon length. Each pool was subjected to library prepa ration using a Nextera XT DNA Library Preparation kit (Illumina, San Diego, USA). All libraries were prepared with dual indices and sequenced with pairedend, 300bp reads on a MiSeq using a 600cycle V3 cartridge (Illu mina, San Diego, USA). The NGS approach has advan tages, by producing massive sequencing data, decreasing costs, and performing high throughput analysis compare to conventional sequencing. Sequencing was conducted in the biotechnology laboratory of the Disease Investiga tion Center Wates, Yogyakarta, Indonesia.

Phylogenetic analysis
Raw reads were imported into the CLC Genomics Work bench (Qiagen) for preprocessing and assembly. Reads were paired and assembled with CLC Genomics set to map reads to reference. The accession number, KF896608, was used as a reference for mapping. Consensus se quences were extracted from the assemblies for all sin gle sequences. Genespecific consensus sequence align ments of five different regions (5′ UTR, NPro, E2, NS3, and NS5B) were constructed using MUSCLE in MEGA X. Phylogenetic analysis of the five regions was performed utilizing the MaximumLikelihood method and bootstrap test (n = 1,000). The BVDV sequences available in Gen Bank were aligned with the sample sequence using blastn for phylogenetic analysis.

Result
Five BVDV isolates were sequenced using the NGS tech nique with MiSeq Illumina. The sequencing results were then assessed for five different regions, namely 5′ UTR, NPro, E2, NS3, and NS5B to determine the BVDV subgenotype. Details of the sequencing results for the five samples are shown in Table 1. Based on the sequencing results, the full length of the 5′ UTR and NPro genes could be analyzed, whereas there were two to four samples for genes E2, NS3, and NS5B. Sample sequences with only 30-45% of whole sequences were less representative of the region. More over, these five regions are categorized into untranslated regions (5′ UTR), structural genes (E2), and nonstructural genes (NPro, NS3, and NS5B). The phylogenetic tree of BVDV1 based on the 5′ UTR along 275 nucleotides (nt) against five BVDV isolates are shown in Figure  1 The nucleotide sequences of five field viruses and a number of known reference strains originating from GenBank were aligned and phylogenetic trees were constructed. Phylogenetic trees based on nonstructural regions of NPro along 504 nt against five BVDV isolates were shown in Figure 2. Two subgenotypes1a and 1c were also divided the same as phylogenetic tree of 5'UTR. The BVDV1a branch was supported by 97% bootstrap value. It showed three isolates V4 BVDV1/Indonesia CJBms/HST0783/2014, V7 BVDV1/IndonesiaCJ Bms/0415116716/2015, and V20 BVDV1/Indonesia  The sequences of 2,049 nt of nonstructural region of NS3 were determined and analyzed with 18 sequences data from GenBank ( Figure 3). Two Indonesian BVDV 1 isolates could be analyzed in full length. There were two branch in this phylogenetic tree constructed the same branch as that of 5'UTR and NPro regions, BVDV1a and BVDV1c with a bootstrap value of 97%. One isolate, V7 BVDV1/IndonesiaCJBms/04151167 16/2015 was clustered as subgnotype BVDV1a with NC001461/BVDV1/NADL and other BVDV1a refer ences. Isolate V1 BVDV1/IndonesiaEJPsn/04131368 2/2013 from different district clustered phylogenetically with BVDV1c references, KF896608/Australia/Bega like/2012 and AF052304 Pest type 1 strain Trangie. The sequence of three isolates could not be analyzed in NS3 re gion because the length of sequence had no good results, either no sequence or less than 50%.
The phylogenetic tree based on the nonstructural regions of NS5B along 2,157 nt against four BVDV isolates and 21 sequences originating from Genbank is shown in Figure 4. The phylogenetic tree of this region also showed the same topology as the 5'UTR, NPro, and NS3. There were four Indonesian BVDV1 isolates which divided into two subgenotypes1a and 1c with a bootstrap value of 97%. Subgenotype1a showed three isolates V4 BVDV1/IndonesiaCJBms/HST0783/2014, V7 BVDV1/IndonesiaCJBms/0415116716/2015, and V20 BVDV1/IndonesiaCJBms/0416164414/2016 were in one group with AF091605/BVDV1/Oregon C24V, NC001461 BVDV1/NADL, and AF145364 BVDV1/Singer. One isolate, V1 BVDV1/IndonesiaEJ Psn/041313682/2013 was included in subgenotype1c with KF896608/Australia/Begalike/2012. There was one isolate that could not be analyzed on this region because of less than 50% full length sequence.
The structural E2 region was the fifth region to be analysed of the following 1,122 bases. There were four Indonesia BVDV1 isolates could be analysed, however one isolate was less than 50% of full length sequence. The phylogenetic tree based on the structural region of E2 along 1,122 nt against Indonesian BVDV1 isolates and 20 sequences originating from Genbank is shown in Figure 5. This phylogenetic was constructed using the same references that generated for the first four phylogenetic analysis. Two subgenotypes, BVDV1a and BVDV1c were obtained from phylogenetic tree of E2 region BVDV1 with a bootstrap value of 99%. There was one isolate V20 BVDV1/IndonesiaCJBms/0416164414/2016 was featured different subgenotypes, mostly 1c in the 5'UTR, NPro, NS3, and NS5B, while their E2 region were clustered into 1a with 96% bootstrap value.

Discussion
Classification of subgenotypes of five BVDV isolates based on the five regions, 5 UTR, NPro, NS3, NS5B, and E2, resulted in two subgenotypes BVDV1a and BVDV1c, as per previously published studies in Indone sia (Wuryastuti et al. 2018). The 5'UTR phylogenetic tree using Australian BVDV reference sequences illus trated that three out of five Indonesian viruses belonged to the BVDV1c. The other subgenotypes, 1b and 1d were not found in this study even though isolated from Central Java as previous studies (Saepulloh et al. 2015).
According to Yeşilbağ et al. (2017), the comparison of global BVDV1 distribution showed a greater distribution of BVDV1b followed by BVDV1a, whereas BVDV1c was the most dominant subgenotype in Australia (Reichel et al. 2018). Comparison of BVD1a subgenotypes based on 5' UTR is threefold greater than BVD1c in Java, In donesia (Wuryastuti et al. 2018).
The conserved region could be targeted to de tect various subgenotypes of BVDV (Nagai et al. 2004). Further genetic studies using American isolates (L32875 BVDV1/Singer, NC001461 BVDV1/NADL, and AF091605 BVDV1/OregonC24V) and Australian isolates (KF896608 BVDV1/Australia/Begalike1c/2012, AF049221 BVDV1/Australia/Bega/2001 and AF049222 BVDV1/Australia/TrangieY546/2001) as references have shown that clustered BVDV1a and BVDV1c respectively. These findings are highly significant as published previously (Yeşilbağ et al. 2017) and also whole genome (Bazzucchi et al. 2017) (Ridpath 2005). The structural region, E2 as hypervariation region used to in this study for the phylogenetic analysis and could be used to predict the rate of mutations in the evolution and virulence of the virus as well as factors that influence the pathogenesis of the virus in causing disease (Goens 2002; Dow et al. 2015. A collection of five BVDV isolates from persis tently infected dairy cow samples was identified using realtime BVDV RTPCR. Persistently infected cattle were the source of the virus, and thus have the potential to transmit BVDV; therefore, they excluded from the herd group (LiebnerTenorio 2005 In this study, based on the time sequence, the 2013 virus isolate was identified as the BVDV1c subgenotype, whereas in 2015, besides being identified as 1c, the iso late was also grouped in 1a. V20 BVDV1/IndonesiaCJ Bms/0416164414/2016 isolates from 2016 cases showed the existence of both subgenotypes. According to Na gai et al. (2004), recombination between virus strains, classified into different subgroups, took place in animals through different subgenotypes. Genetic diversity can be caused by mutations, errors in RNAdependent RNA poly merase activity, and recombination of homologous RNA and heterologous viral genomes; however, it is very im portant for taxonomy, laboratory diagnosis, and vaccine design (Vilcek et al. 2004).
Recombination in Pestiviruses has been characterized in BVDVs. The recombination can occur between the host and virus, producing a host RNA sequence that in serts into the viral genome or through rearrangement, du plication, or both, of the viral genome sequences (Goens 2002). Heterolog/nonhomologous recombination in Pes tiviruses is reported to produce cytopathic biotypes that can evolve from noncytopathic viruses and cellular se quences (Becher and Tautz 2011).
Homologous RNA recombination of the BVDV genome in structural regions has been observed in trials of animals infected with mucosal diseases (MDs) (Na gai et al. 2003). According to Peterhans and Schweizer (2010), persistently infected cattle are at risk of being in fected with a mutant or different biotype BVDV, thus in creasing their potential to have MD. The discovery of two subgenotypes, 1a and 1c, in V20 BVDV1/IndonesiaCJ Bms/0416164414/2016 originating from persistently in fected cattle indicates recombination of the viral genome and the potential to cause MD events. Phylogenetic analy sis of various nonstructural and structural proteincoding regions is useful for virus characterization and further epi demiological analysis.

Conclusions
Phylogenetic analysis of five BVDV isolates collected from persistently infected cattle was carried based on five different regions: the 5' UTR, nonstructural regions (NPro, NS3, NS5B), and the E2 structural region, and showed two subgenotypes, specifically 1a and 1c. The analysis of one of these virus isolates also indicated ge netic recombination of BVD subgenotypes 1a and 1c in persistently infected dairy cattle.