Genetic and Antigenic Characterization and Retrospective Surveillance of Bovine Influenza D Viruses Identified in Hokkaido, Japan from 2018 to 2020

Influenza D virus (IDV), which is a new member of the Orthomyxoviridae family, is potentially involved in bovine respiratory diseases (BRDs). Bovine IDVs (BIDVs) from Japan have been distributed nationwide since 2010 and are genetically distinct from foreign IDVs. We isolated BIDVs from three BRD outbreaks, in Hokkaido during 2018–2020, to understand their genetic and antigenic characteristics. Retrospective surveillance was performed using sera collected throughout the last decade in Hokkaido to investigate BIDV existence. Three BIDVs were isolated using cell culture. Comparative and phylogenetic analyses using sequence data of the three BIDVs and IDVs from Japan and other countries available in GenBank demonstrated that Japanese BIDVs, including the three BIDV isolates, were genetically distinct from other IDVs. Genotype classifications based on the rotavirus genotype classification revealed multiple genotypes of RNA segments 1–7. Two BIDVs were of a new genotype, different from those of other Japanese BIDVs. Neutralization assays against two BIDVs with different genotypes using sera collected in acute and recovery phases of BRD revealed differences in cross-reactivity to heterogenous BIDVs. Retrospective surveillance suggested that BIDV existed in Hokkaido, in 2009. Our findings suggest that BIDVs of different genotypes and antigenicity are distributed and maintained in Hokkaido and provide new insights into molecular characteristics and the evolution of IDVs.


Introduction
Influenza viruses are enveloped, segmented, single-stranded, negative-sense RNA viruses, which belong to the family Orthomyxoviridae, and are currently classified into the following four species: influenza A, B, C, and D (IAV-IDV). The genomes of IAV and IBV consist of eight RNA segments, whereas ICV and IDV have seven segments. Both IAV and IBV contain two major surface glycoproteins, i.e., hemagglutinin (HA) which binds to host cell receptors and mediates membrane fusion, and neuraminidase (NA) which cleaves receptor sialic acids resulting in the release of newly assembled virus particles [1]. In contrast, ICV and IDV have only one major glycoprotein, the hemagglutinin-esterase-fusion (HEF) protein, which possesses "all-in-one" activities of receptor binding, receptor cleavage, and membrane fusion [2][3][4][5]. The IDV HEF glycoprotein is structurally and functionally similar to the ICV HEF glycoprotein and is closely associated with the antigenicity and pathogenicity of the virus [5,6].
IDV was first isolated from swine exhibiting influenza-like symptoms in Oklahoma, United States of America (USA), in 2011 [7]. Subsequent observations have demonstrated that cattle, rather than pigs, were the natural reservoirs of the virus [8,9]. Experimental infection indicated that the bovine IDV (BIDV) caused mild respiratory symptoms in cattle [10,11]. However, metagenomic analyses of feedlot cattle suffering from bovine respiratory disease (BRD) suggested that this virus was a causative pathogen of BRD [12,13]. Moreover, a serological survey revealed that, apart from cattle and pigs, sheep, goats, horses, and camels were also susceptible to infection with IDV, suggestive of its broad cell tropism [14][15][16]. Furthermore, other serological studies have reported the detection of virus-specific antibody titers in human sera [7,17,18]. Thus, these findings suggest that this virus has zoonotic potential.
In this study, we found that BIDVs were associated with BRD outbreaks the occurred in Hokkaido, in recent years. We also demonstrated their genetic and antigenic characteristics via phylogenetic and antigenic analyses of the BIDVs obtained in this study and other IDVs reported in previous studies. Furthermore, we performed a retrospective analysis of the existence of BIDV in Hokkaido using sera collected and stored in the last 10 years and observed that BIDV existed as early as 2009.

Samples, RNA Extraction, Diagnostic Test, and Ethics Statement
To investigate whether BIDV was associated with BRD, we diagnosed cattle from three BRD outbreaks in herds occurring in Hokkaido, Japan from 2018 to 2020, as follows: In Outbreak 1, seven of 26 calves (<3 months of age) developed fever and respiratory distress at cattle farm A on November 2018. Thereafter, BRD slowly spread across most cattle in herds, between November 2018 and January 2019. In Outbreak 2, 80 of 90 calves (1-4 months of age) developed fever, severe cough, and nasal discharge at cattle farm B, in January 2019. Furthermore, some cattle developed pneumonia and became debilitated. In Outbreak 3, 300 weaned calves (2-3 months of age) developed nasal discharge at cattle farm C, in October 2019. Subsequently, 58 of the cattle died due to pneumonia, by early January 2020.
Viral RNA was extracted from the nasal swab samples of five or six individuals from each outbreak using the High Pure Viral RNA kit (Roche, Basel, Switzerland), according to the manufacturer's instructions. These samples were subjected to conventional reverse transcription-PCR (RT-PCR) analysis specific for bovine viral diarrhea virus (BVDV), bovine respiratory syncytial virus (BRSV), and bovine coronavirus (BCoV) based on the recent prevalence of viruses associating with BRD, in Japan, reported in previous studies using the PrimesScript One Step RT-PCR Kit Ver. 2 (Takara, Shiga, Japan), according to previously reported methodologies [28][29][30][31][32]. In addition, all samples were subjected to real-time (RT)-PCR (qRT-PCR) specific for BIDVs, according to a previous report [31]. Obtained samples were also tested for the presence of three bacterial species (Mycoplasma bovis, Pasteurella multocida, and Mannheimia haemolytica) based on the recent prevalence of bacteria associating with BRD in Japan in previous studies, according to a routine methodology [31][32][33][34].
All samples were collected as a part of routine diagnostic procedures, hence, permission with regard to animal ethics was not required.

Sequence and Phylogenetic Analyses
Viral RNA was extracted from the supernatant of each BIDV isolate originating from the three BRD outbreaks, according to the methodology described above. Genomic sequences of individual RNA segments were amplified via RT-PCR, using a set of primers originally designed by reference to the sequences of other IDVs available in GenBank (Table S2). The primer sets were confirmed to not amplify PCR products from other species of IVs. RT-PCR was carried out using the PrimeScript II High Fidelity One Step RT-PCR Kit (Takara, Shiga, Japan) with the following cycling conditions: 45 • C for 10 min and 94 • C for 2 min; 35 cycles of 98 • C for 10 s, 55 • C for 15 s, and 68 • C for 30 s; final extension step at 68 • C for 7 min. PCR products were sequenced using the BigDye Terminator v3.1 Cycles Sequencing Kit on an automated ABI prism 3130 Genetic Analyzer (Thermo Fisher Scientifics, Carlsbad, CA, USA). Each genomic sequence from the three BIDVs determined herein has been submitted to the DNA Data Bank of Japan and is retrievable via GenBank (Table S1, GenBank accession number LC565467-LC565487).
The sequence data were aligned using the Clustal W method in the MEGA X software [36]. Genetic distances for the seven RNA segments were calculated using the Kimura two-parameter correction at the nucleotide level. Phylogenetic analyses for all RNA genomes, including the three BIDVs and other previously reported IDVs, were performed using the maximum-likelihood method with the general time reversible nucleotide substitution model and 1000 bootstrap replicates implemented in the MEGA X program [36]. Genotype classification of individual IDV genes was conducted using cut-off values calculated based on the definition that was used in the full genome-based genotype classification of rotavirus, a segmented RNA virus like IDV [37,38]. Briefly, the cut-off values for all genes were estimated as the percentage separating intra-genotype identity (nucleotide identity among strains belonging to the same genotypes), and inter-genotype identity (nucleotide identity among strains belonging to different genotypes). However, in cases when inter-and intra-genotype identity partially overlapped, the most appropriate cut-off value was chosen as the percentage at which the ratio of inter-genotype identity and intra-genotype identity dropped below 1.

Comparison of Antigenicity between Two BIDV Isolates Using A Neutralization Assay
We collected acute (pre) and recovery (post) phase serum samples from 8 and 5 cattle from farms A and B, respectively. To compare antigenicity of BIDVs (HKD1, D/bovine/Hokkaido/HKD1/2018 and HKD2, D/bovine/Hokkaido/HKD2/2019) isolated from the two farms, we performed cross-reactive neutralization tests using these serum samples. Heat-inactivated sera (50 µL) were serially two-fold diluted with DMEM and mixed with an equal volume of 200 TCID 50 of HKD1 or HKD2, at 37 • C, for 1 h. Then, the mixture was added to HRT-18G cells (2.5 × 10 4 cells/100 µL per well in 96-well plates), and cells were incubated, at 37 • C, for 10 days. On the basis of microscopic observation, the highest dilution of sera completely protecting the cells from CPEs was recorded as the viral neutralizing (VN) antibody titer.
Antigenic cross-reactivity using pre-and post-BRD outbreak serum samples from both farms against two different BIDV isolates were compared and analyzed with the Wilcoxon rank sum test. A p-value < 0.05 indicated a significant difference.

Retrospective Surveillance
We carried out neutralization assays against HKD1 and HKD2 using a total number of 960 serum samples collected from cattle (over 24 months of age) selected randomly at 96 different farms in Hokkaido, each year during 2009 to 2018, in order to investigate the existence of BIDV in Hokkaido in the past (Table S3).

Diagnosis of Cattle from Three BRD Outbreaks
We diagnosed cattle from each BRD outbreak occurring at three cattle farms in Hokkaido between 2018 and 2020. We tested five or six nasal samples from each outbreak through pathogen-specific RT-PCR and BIDV-specific qRT-PCR, and virus and bacteria isolation (Table 1). In Outbreak 1, three viruses (BCoV, BRSV, and BIDV) were detected by RT-PCR and qRT-PCR. Furthermore, three viruses (BPIV3, BCoV, and BIDV) and two bacterial species (P. multocida and Myc. bovis) were isolated in cell culture and agar, respectively. In Outbreak 2, BRSV and BIDV were identified by RT-PCR and qRT-PCR, respectively. In addition, BIDV and M. haemolytica were isolated using HRT-18G cells and agar, respectively. In Outbreak 3, three viruses (BRSV, BCoV, and BIDV) and three bacteria (P. multocida, Myc. Bovis, and M. haemolytica) were detected in the six used nasal samples tested. Moreover, BCoV and BIDV were also detected in HRT-18G cell culture. In summary, we isolated five BIDVs from all BRD outbreaks, which were caused by multiple viruses and bacteria.
The two BIDV isolates (HKD1 and HKD2) were confirmed by a HA assay using turkey RBCs and TEM observation (Table S1 and Figure S1). In addition, viral titers of HKD1, HKD2, and HKD3 determined according to the Reed and Muench's method were 8.0, 8.1, and 7.8 TCID 50 / mL , respectively (Table S1).

Sequence and Phylogenetic Analyses
Amplification by RT-PCR, using a set of primers originally designed by reference to the complete genomes of other IDVs available in GenBank, successfully determined the nearly full-length nucleotide sequences of all RNA segments, excluding several nucleotides at the 5' and 3' termini, of the five BIDVs isolated in this study (Table S1). We defined one of them as a representative strain, HKD2, because the nucleotide sequences of seven RNA segments of three BIDVs isolated from Outbreak 2 were identical (data not shown). The lengths of the open reading frames (ORFs) of all genes of the three BIDVs (HKD1, HKD2, and HKD3) were almost identical to those of reference BIDVs without insertions and deletions.
Comparative sequence analyses among the three BIDVs identified in this study, as well as among these and other BIDVs from Japan detected in previous studies, demonstrated Japanese BIDVs had high genetic diversity, especially in HEF gene (Table 2). A comparison of the nucleotide sequences of the seven RNA segments of Japanese BIDVs with those of IDVs from other countries revealed that Japanese BIDVs are genetically distinct from IDVs from other countries.  Phylogenetic analyses using ORFs of individual genes were performed by adding data from the three BIDVs to other IDV data available in GenBank (Figure 1). In addition, we also carried out the phylogenetic analyses using ORF nucleotide sequences of NS1 and NS2, because NS gene encodes two proteins (NS1 and NS2). However, these dendrograms revealed similarity with regard to that of NS gene (data not shown). According to the definition for genotype classification of rotavirus, cut-off values for the genotype classification of the PB2, PB1, P3, HEF, NP, M, and NS genes were calculated from the frequency distribution of pairwise sequence identities and were set to 97.5%, 97.2%, 97.6%, 97.4%, 98.1%, 97.8%, and 98.1%, at the nucleotide level, respectively (Table 3 and Table S4). On the basis of these cut-off values, we revealed the existence of 5, 4, 5, 6, 7, 4, and 5 genotypes for the PB2, PB1, P3, HEF, NP, M, and NS genes, respectively ( Figure 1 and Table 3). In the analysis of all genes, BIDVs from Japan were clearly distinct from IDVs from other countries and were classified into one or two genotypes. In addition, the HKD2 and HKD3 were grouped into a new genotype, different from the genotype of other Japanese BIDVs, including HKD1, based on analyses of the PB2, P3, HEF, NP, and NS genes. Moreover, HKD3 was classified into a PB1 genotype different from other IDVs.

Comparison of Antigenicity between Two BIDV Isolates Using Neutralization Assay
There were significant differences (serum samples from farm A, p = 0.003 and serum samples from farm B, p = 0.026) in cross-reactivity to heterogenous BIDV isolates in the neutralization assay using serum samples from farms A and B (Table 4). Briefly, the neutralization assay using serum samples collected from pre-and post-BRD outbreak at farm A showed a clear increase of VN antibody titers against HKD1 isolated from farm A, but not against HKD2 isolated from farm B. In contrast, the assay using five serum samples from pre-and post-BRD outbreak at farm B exhibited high increases of VN antibody titers against HKD2 but smaller increases of titers against HKD1.

Retrospective Surveillance
Retrospective surveillance was performed through a neutralization assay for two BIDV isolates, HKD1 and HKD2, using 960 serum samples collected in Hokkaido over the past decade. When considering VN antibody titers of more than 10 as positive, the detection rates of antibodies against BIDV HKD1 and HKD2 were 55% and 57%, respectively, almost identical between both assays (coincidence rate = 96%, 515 and 406 samples were positive and negative against both HKD1 and HKD2, respectively; 16 samples were positive against HKD1, not HKD2; and 23 samples were positive against HKD2, not HKD1). Our analysis revealed that antibodies against BIDV in cattle sera have been continuously (in a range from 45% to 71%) detected in Hokkaido from 2009 to 2018 (Table 5). In addition, our analysis also suggested that the type of dominant virus tended to fluctuate between 2009 and 2018. However, we could not confirm the trend, because the VN antibody titers were affected by timing of sample collection and individual differences. Table 4. Viral neutralizing antibody titers of serum samples collected in the acute (pre) and recovery (post) phases of BRD outbreaks that occurred at farms A and B against bovine influenza D viruses (HKD1 and HKD2) isolated from the two farms, as measured using a neutralization assay.  2018 and 2020. Interestingly, we also isolated BIDVs from samples by using an HRT-18G cell culture. These findings strongly support the notion that BRD is caused by interactions between viruses and bacteria and that BIDV is one of the causative pathogens of BRD, as reported in previous studies [12,13]. Sets of primers originally designed by reference to sequences at the 5' and 3' termini of the segments of other IDVs available in GenBank successfully amplified the PCR products of seven nearly full-length RNA segments from the three BIDVs obtained in this study. This suggests that the nucleotide sequences of the 5' and 3' termini are conserved between different IDVs, as commonly observed for segmented RNA viruses [39].

Viral Neutralizing Antibody Titers for Serum Samples from Farm
Comparative sequence analyses among the three BIDVs isolated in this study and IDVs from Japan and other countries reported in previous studies suggested that these BIDVs had large genetic variation. Furthermore, HKD2 and HKD3 of the three BIDVs are genetically distinct from other IDVs. Genotype classifications following the phylogenetic analyses for all RNA genomic sequences of the three BIDVs and other IDVs revealed the existence of multiple genotypes for individual genes based on cut-off values calculated according to the definition of rotavirus genotype classification [37,38]. In addition, the analyses of all genes revealed that the genotypes to which six Japanese BIDVs belonged, including the three BIDVs from the current study, were clearly different from the genotypes of other IDVs, as reported in previous studies [24,27]. Moreover, the data presented in this study classified HKD2 and HKD3 into a new genotype, which was distinct from genotypes that other IDVs belonged to, for the five or six remaining genes, except for the PB1 or M genes. Taken together, we successfully found a novel BIDV with a unique genotype via a series of genetic analyses of BIDVs identified, in Hokkaido, in recent years.
In the three BRD outbreaks that occurred in Hokkaido between 2018 and 2020, cattle with BRD from Outbreaks 2 and 3 had more severe symptoms (a rapid spread in cattle herds, debilitation, and death due to pneumonia) than cattle from Outbreak 1. In addition, cross-reactive neutralization assays against two BIDVs (HKD1 and HKD2) using serum samples collected from infected cattle in farms A and B revealed significant differences in cross-reactivity to heterogenous BIDVs, suggestive of the antigenic heterogeneity between the two BIDVs. Moreover, the phylogenetic analyses and genetic classifications of the seven RNA segments revealed that HKD1 had a different genetic background from the two other BIDVs (HKD2 and HKD3). Especially, the three BIDVs were clearly classified into two different genotypes of the HEF gene, which have been closely associated with the antigenicity and pathogenicity of the virus [5,6]. These findings suggest that several kinds of BIDVs with different pathogenicity, antigenicity, and genotypes have been distributed and maintained in Hokkaido, Japan.
Retrospective surveillance using 960 sera collected in Hokkaido since 2009 revealed that antibodies against BIDV in sera had been detected every year for the last decade. This observation suggests that BIDV has existed in Hokkaido since 2009, which is earlier than previously reported [23].
In conclusion, we isolated three BIDVs from infected cattle with mild to severe respiratory diseases during BRD outbreaks, in Hokkaido, Japan in recent years. Two BIDVs isolated from cattle with severe symptoms were classified into a new genotype, different from the genotype of previously described Japanese BIDVs, especially with regard to the IDV HEF gene. In addition, the neutralization tests using BIDVs of two different genotypes suggested a significant difference in antigenicity between the two BIDVs. Moreover, our data demonstrated that the BIDV had been distributed in Hokkaido at least since 2009. The genetic classification of IDVs performed in this study will provide useful information for the monitoring and identification of a new IDV, which could emerge in the future.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4915/12/8/877/s1, Figure S1: Image of BIDV virus, which isolated from HRT-18G cell culture, by transmission electron microscopy observation. The nasal swab sample collected from cattle with respiratory disorder was inoculated into HRT-18G cells, and cells were kept, at 37 • C, for a week. After a week, the supernatants were harvested from the cell culture which exhibited cytopathic effects, and then observed by transmission electron microscopy. Bar indicates 100 nm. Table S1: Summary of HA and virus titers, and length of open reading frame determined by a sequence analysis and GenBank accession number of three bovine influenza D viruses isolated from this study, Table S2: A set of primers for genomic sequence determination of seven RNA segments from bovine influenza D virus, which were originally designed by reference to other influenza D viruses available in GenBank, Table S3: Summary  of 960 serum sample collected between 2009 and 2018 used in retrospective surveillance, Table S4: Open reding frame (ORF) nucleotide sequence identities among genotypes on individual RNA segments of influenza D viruses. Funding: This study was partially supported by a grant from National Agriculture and Food Research Organization (NARO).