Genetic Evolution of Bovine Coronavirus Detected in the Republic of Korea

Bovine coronavirus (BCoV) is associated with severe diarrhea in calves, winter dysentery in adult cattle, and respiratory diseases in cattle. However, there is limited information regarding its molecular characterization in the Republic of Korea (KOR). This study was conducted to investigate the prevalence of BCoV in diarrheic calves, perform global comparison of BCoV genome sequences, and suggest a new nomenclature. A total of 521 fecal samples were collected from diarrheic pre-weaned native Korean calves (age £60 days) from 100 herds in the KOR. The samples were divided into three groups based on age as follows: 1-21 days (n = 433), 21-40 days (n = 64), and 41-60 days (n = 24). The full-length Spike (S) gene was amplied and sequenced. The phylogenetic tree was constructed using the maximum-likelihood method and the evolutionary rates were estimated. BCoV infection was detected in 20 (3.8%) calves by real-time RT-PCR analysis, and nine full-length genome sequences were obtained from the 20 BCoV-positive samples. Genomic comparison analysis showed that these 2019-2020 variants exhibited the highest nucleotide sequence identity (98.6%-99.2%) with that of water deer (Hydropotes inermis) isolates. Phylogenetic analysis of the full genomes of the S gene revealed the following four BCoV groups: G1, classical BCoV strains; G2, 2002 Korean, Japan, Vietnam, Cuba, and USA strains/isolates; G3, European strains/isolates; and G4, Korean isolates (2004 and 2006 Korean isolates, 2010 BCoV-like, 2017 water deer, and 2019-2020 variants). The evolutionary rates accelerated from G1 to G4. This grouping was also closely related to the nucleotide substitution rate. Using molecular clock analysis of the S gene, the most recent common ancestor of each group was estimated to have originated in 1953, 1979, 1986, and 1993, respectively. Recently identied BCoV variants in the KOR are undergoing slow evolution. These ndings provide useful information for understanding the molecular characterization of BCoVs. Further research is necessary to conduct recombination analysis to support BCoV evolution.

a reservoir host for other domestic and wild ruminants [21]. Considering the current situation regarding the coronavirus pandemic, the zoonotic transmission potential of BCoV due to the rapid evolution of CoVs cannot be excluded.
In the Republic of Korea (KOR), research on BCoVs was primarily conducted in the 2000s [28][29][30][31][32], and there is limited information regarding their more recent genetic variation and evolution. Therefore, the classi cation criteria for BCoV have not yet been clearly established. To date, several studies have focused on the interspecies transmission of BCoV [27,[33][34][35]. Hence, the present study was conducted to investigate the prevalence of BCoVs in diarrheic calves and to compare the genomic sequences of BCoVs distributed worldwide, as well as to suggest a new nomenclature.

Sample collection
Between January 2019 and June 2020, a total of 521 fecal samples were collected by a veterinarian directly from the rectum of pre-weaned native Korean calves (age £60 days) suffering from diarrhea from 100 different herds in the KOR. The collected samples were divided into the following three groups according to age: 1-21 days (n = 433), 21-40 days (n = 64), and 41-60 days (n = 24). These samples were placed in plastic containers and immediately transported to the Animal Immunology Laboratory at Kyungpook National University. Most calves were not vaccinated against BCoV.

RNA extraction and real-time PCR
Total RNA was extracted from fecal suspensions using RNAiso Plus Reagent (Takara, Shiga, Japan) according to the manufacturer's instructions. Real-time polymerase chain reaction (RT-PCR) was conducted to detect BCoV infection using the Path-IDTM Multiplex One-Step RT-PCR Kit (Life Technologies, Carlsbad, CA, USA) as described previously [36]. First, partial ampli cation (824 bp) of the S gene was performed as described previously [29,37]. Next, the ampli cation values (Ct values) were determined for each sample. Ct values £35 indicated a positive test result.
RT-PCR, sequencing, and phylogenetic tree analysis To obtain sequences, positive samples detected by real-time PCR assay were further ampli ed using 2×AccuPower® RocketScript™ RT-PCR Master Mix (Bioneer, Daejeon, KOR) with the full-length S gene (4092 bp) of six primer pairs as described previously [29]. Distilled water was used as a negative control. The S1 to S6 genes of BCoV were attached with overlapping sequences. The PCR products were separated by gel electrophoresis on a 1.5% agarose gel stained with ethidium bromide and visualized. All PCR products were puri ed using the AccuPrep PCR Puri cation Kit (Bioneer) and then used for direct sequencing (Macrogen, Daejeon, KOR). The nucleotide sequences of the 2019-2020 variants obtained in this study were aligned using ClustalX (version 1.8) and compared with 81 BCoV isolates/strains obtained from the GenBank database. The phylogenetic tree was constructed based on the nucleotide alignments using the maximum-likelihood method implemented in the MEGA 7 software [38], and bootstrap analysis was used to evaluate the robustness with 1000 replicates. The cut-off values for bootstrap replication were ³70%. The nucleotide sequences of the 2019-2020 variants obtained in this study were assigned the accession numbers MW881213-MW881221.

Evolutionary rate
The evolutionary rates of nucleotide per site and the year of virus sampling based on the S gene were estimated using the Bayesian Markov chain Monte Carlo approach implemented in the Bayesian Evolutionary Analysis Sampling Trees v1.10.4 [39]. The Hasegawa-Kishino-Yano model of nucleotide substitution was speci ed for a proportion of gammadistributed rate heterogeneity with partitions into codon positions of 1 + 2 and 3. The MCMC runs consisted of 8 ´ 10 7 generations, and the posterior probability distribution of the chains was sampled at every 1000 steps. The mean evolutionary rate and mean time of the most recent common ancestor (TMRCA) were calculated. The posterior distribution was summarized in terms of 95% highest posterior density (HPD) after the exclusion of a burn-in equal to 20% of the run length. The effective sample size (ESS) was evaluated using Tracer v1.6.0, and ESS values of ³200 for all parameters were accepted [40,41]. The mutation rates were measured using DNA Sequence Polymorphism (DnaSP) v. 6.12. The mutation rate per RNA sequence per generation was obtained by DnaSP, which was then analyzed using DNA Polymorphism [42,43].

Assessment of genetic relatedness
Principal coordinate analysis (PCoA) was used to determine the relationships among BCoVs. PCoA is a distance-based model, jointly using a dissimilarity matrix calculated with a simple-matching index and factorial analysis. PCoA was conducted using the Dissimilarity Analysis and Representation for windows v. 6.0.021 to reveal the genetic divergence between BCoV strains/isolates through the Kimura method.

Statistical analysis
The PCR results for each fecal sample were recorded as being either positive or negative for BCoV and categorized according to calf age (i.e., 1-20, 21-40, or 41-60 days). Data were analyzed using the SPSS Statistics 25 software package (SPSS Inc., Chicago, IL, USA). The association between BCoV infection and each age range was evaluated using the chi-square test. P < 0.05 was considered signi cant.

Phylogenetic analysis
Based on the phylogenetic tree constructed using the S gene, a total of 90 BCoV isolates/strains, including BCoV-like and the nine above-described sequences, were divided into four groups (Fig. 2). Currently, there are no accurate classi cation criteria for BCoVs, so they were randomly named (G1-G4) in this study. Our sequences were most closely related to those of the water deer isolate (MG518518) and clustered with the BCoV-like isolates found in zoo ruminants ( Fig. 2). Furthermore, these sequences formed a different subgroup from the 2004 KCD isolates and 2006 KWD isolates within the same group (Fig. 2). Although it was detected in the KOR, the isolates/strains from 2002 belonged to a different group from those detected after 2004 (Fig. 2). Including KCD isolates from 2004, KWD isolates from 2006, BCoV-like isolates from 2010, and the contemporary 2019-2020 variants, these Korean isolates/strains were designated as belonging to G4 (Fig. 2). The KWD isolates found in 2002 formed the same clade with that of the whitetailed deer (FJ425187) and were separated from BCoV and BCoV-like isolates identi ed from Japan, Vietnam, Cuba, and USA, all of which were designated as G2 (Fig. 2). G3 was primarily composed of BCoV strains/isolates detected in Europe, which were slightly different from those found in Asia and the USA (Fig. 2). G1 represented the classical BCoV strains that share little genetic relationship with other BCoVs (Fig. 2). The results of PCoA revealed four well-separated clusters among 90 BCoV isolates/strains, including BCoV-like isolates (Fig. 3), which corresponded to the results from the phylogenetic tree. Speci cally, the PCoA analysis demonstrated clear divergence within the G4 group, i.e., 2004-2006 BCoVs vs. 2019-2020 BCoVs.

Evolutionary rate
The mean evolutionary rate and TMRCA based on the established phylogenetic tree were analyzed ( Table 3). The evolutionary rate of 68 BCoV global strains/isolates from 1965 to 2017, except for BCoV-like and the 9 variants described in this study, was estimated to be 0.3897 ´ 10 -3 substitutions/site/year (95% HPD; 0.2-0.5 ´ 10 -3 ), and the TMRCA of 68 BCoVs was dated to be 1947 (95% HPD; 1939-1954). As shown in Table 3, the nucleotide substitution rates of G1 and G2 are relatively low, whereas within G3 and G4, there is a signi cant variation in the nucleotide substitution rates. G4 was composed of Korean BCoV strains/isolates and showed the presence of signi cant variation in the nucleotide substitution rates compared with the evolutionary rate of European BCoVs ( Table 3). The TMRCA of classical BCoV strains (G1) was accordingly estimated to have originated in the year 1953, and that of G2 consisting of the USA and Asian strains/isolates was in 1979. The European BCoV strains (G3) appeared in 1986. The TMRCA of G4, including the sequences described in this study as well as Korean BCoV-like and BCoV isolates/strains, was estimated to be from 1993 (95% HPD; 1990-1996) ( Table 3

Discussion
This study demonstrated a relatively low prevalence of BCoV in diarrheic calves. This nding was similar to that (5.4%) of our previous study [44] and those from several countries [17][18][19][20]45]. Compared with the past, the incidence of BCoV appears to have constantly decreased. Although there was no statistical signi cance between diarrhea caused by BCoV and calf age, the incidence of BCoV infection was highest in calves aged < 20 days. This result was consistent with a previous study performed by our group [46]. Although the incidence of BCoV was low in the present study, the importance of BCoV as the primary causative agent of diarrhea in calves should not be ignored because it was still detected in calves with diarrhea. Especially, recent studies report that BCoV is also considered as a key pathogen associated with respiratory diseases [47][48][49][50][51]. Except for the KOR, other countries have focused on BCoV as a respiratory pathogen rather than as a diarrheal pathogen [52][53][54][55]. Further study is required to explore the prevalence of BCoV by its associations with various clinical symptoms and to identify the genetic characteristics between enteric and respiratory BCoVs circulating in the KOR.
We compared the full genome sequences of 42 Korean BCoV strains/isolates, including recent BCoVs (2019−2020) and detected genetic variations within nine sequences. Our results clearly demonstrated that the genetic variation appeared signi cantly over time (Fig. 1). Due to a lack of information on BCoVs since 2006, it is unclear when these genetic variants of BCoV have occurred in the KOR; nevertheless, our results indicate that the genetic variation of BCoV has been slowly occurring since a long period of time. Furthermore, we cannot conclude whether these variations are involved in the pathogenicity of BCoV. Further investigations are required to evaluate the association between genetic variations and pathogenicity.
The 2019−2020 BCoV variants were most similar to the isolates obtained from the water deer [24]. Water deer are widely distributed in the KOR and are known to transmit various diseases [24]. As the scarcity of food and natural habitat has resulted in frequent invasion to farmlands, the possibility that BCoV was transmitted from cattle to water deer or vice versa cannot be ruled out. Interestingly, the sequences found in the present study exhibited 98.7%−99.2% nucleotide identity with BCoV-like isolates identi ed in zoo animals suffering from diarrhea, and these BCoV-like isolates were known to be derived from cattle in 2005 [56], but the transmission route of the virus to zoo animals remains unknown. Moreover, in the United States, BCoVs have been detected in various wild ruminants, such as giraffe, antelope, alpaca, and deer [33, [57][58][59]. However, when the results of genetic analysis were compared between captive wild ruminants from the USA and the zoo ruminants from the KOR, signi cant differences were observed. It is speculated that the viral host reservoir differs between these viruses. Amer reported that the difference between BCoVlike isolates and BCoV is due to the host range, and not distinct virus species [4]; however, our results demonstrated that there were differences between BCoV-like and recent BCoVs variants in genome alignments (data not shown), indicating that the virus species might be different. Our ndings suggest that BCoV can infect non-captive animals as well as captive animals. Further epidemiological studies on wild animals are required to clarify this issue.
The classi cation criteria for BCoV have not yet been clearly established. A recent study reported that BCoV is classi ed into two major types−European and USA−which are subdivided into 11 and 3 genotypes, respectively [60]. However, our results were different from those reported by that study. According to the phylogenetic analysis, BCoVs were largely divided into four groups, G1−G4. There are at least two subgroups within each group. In our study, we could not accurately subdivide further within each group because the basis for subdivision is currently not well-de ned su ciently. G1 includes the classical BCoV strains, including Mebus, which are considered as the ancestors of all BCoVs. G2 consists of BCoV isolates/strains from Asia (Korea, Japan, and Vietnam) and North America (Cuba and USA). Within the G2 group, the 2002 Korean BCoVs were divergent from those of Japan and Vietnam, even though they belong to the same continent. However, the most common feature of G2 is that these BCoV isolates were derived from the US wild ruminants [60]. G3 included strains/isolates identi ed from European countries. This could be due to differences in the evolution of BCoVs, and consequently, it is speculated that this evolution plays a vital role in the genetic diversity of BCoVs from country to country. The contemporary BCoVs (2019−2020) identi ed in this study were classi ed as belonging to G4, and most Korean isolates, except for 2002 BCoVs, were included in this group. These results provide evidence that the 2019−2020 BCoVs circulating in the KOR have undergone recent genetic evolution.
The evolutionary analysis of 90 BCoVs from 1965 to 2020 revealed the genetic variations according to groups and, in particular, over time. When compared with other RNA viruses, the evolutionary rate of BCoVs was within the standard range [61, 62]. Interestingly, a relatively higher rate of nucleotide substitutions was detected in Korean BCoVs, with slightly higher evolutionary rates than those observed in European countries, indicating that the BCoVs detected in the KOR are undergoing more rapid evolution than that in other countries. This is probably because the sequencing information of BCoVs available according to the year of collection in the KOR is rather large compared with that in other countries. The phylogeny shown in this study was found to be strongly associated with this evolutionary rate.
However, due to the lack of whole genome sequences of Korean BCoVs, we could not con rm recombination events in this study. Further studies are required to conduct recombination analysis to support the evolution of BCoVs.
The TMRCA of BCoVs revealed that the USA and Asia share a common origin for BCoVs. The ancestor of BCoVs was historically assumed to be the same and was con rmed via molecular clock analysis to have diverged into the USA and Europe around 1971. Since then, it has been estimated that BCoVs have evolved separately on these continents. Especially, the USA appears to have more impact on the introduction and spread of BCoVs to Asian countries. This rationale may also be explained by the phylogenetic analysis. According to the TMRCA results, Korean BCoVs diverged from the USA in 1994. Hence, the 2002 Korean BCoVs are believed to have been in uenced by those originating from the USA and have been properly adapted to the domestic environment. Unfortunately, we could not perform positive selection analysis within the scope of this study. Further research is necessary to clarify the relationship between positive selection and virus evolution in each country and to provide important information for the development of an effective vaccine against BCoVs.
In summary, our study showed a low prevalence of BCoV in diarrheic calves and the lack of association between BCoV infection and calf age in the KOR. Recently identi ed BCoVs in the KOR are undergoing genetic variation and continuous evolution. These phenomena should be clearly supported through analysis of recombination events and positive selection. We have proposed a new nomenclature (G1−G4) for BCoVs based on the rate of nucleotide substitution. Our ndings provide useful information for understanding the molecular characterization of BCoVs. Further investigation is required to establish an exact classi cation criterion for subdivisions within each group.   Figure 1 Schematic diagram of the genomic sequences of Korean BCoV isolates generated using the Geneious software version 10.2.4. The genomic regions are shown above, with the black bars representing the spike glycoproteins. Lightly shaded regions are identical to the consensus sequence, and the vertical black bars indicate differences from the consensus sequence. The sequences identi ed in this study are indicated by a red line box.

Figure 2
Phylogenetic tree based on the S gene from 90 BCoV isolates/strains reported worldwide, including the BCoV sequences obtained in this study. The tree was constructed using the neighbor-joining method implemented in MEGA version 7. Bootstrap analysis was performed with 1000 replicates. The sequences identi ed in this study are marked with bold circular symbols. Bootstrap values of ≥70 are presented.