# Global transmission, spatial segregation, and recombination determine the long-term evolution and epidemiology of bovine coronaviruses

Elias Salem, Vijaykrishna Dhanasekaran, Herve H. Cassard, Ben Hause, Sarah

Maman, Gilles Meyer, Mariette F Ducatez

► To cite this version:

Elias Salem, Vijaykrishna Dhanasekaran, Herve H. Cassard, Ben Hause, Sarah Maman, et al. Global transmission, spatial segregation, and recombination determine the long-term evolution and epidemiology of bovine coronaviruses. Viruses, 2020, 12 (5), pp.534. 10.3390/v12050534. hal-03160293

HAL Id: hal-03160293

https://hal.inrae.fr/hal-03160293

Submitted on 5 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Global Transmission, Spatial Segregation, and Recombination Determine the Long-Term Evolution and Epidemiology of Bovine Coronaviruses

Elias Salem<fnref n="1" />, Vijaykrishna Dhanasekaran 2,3✉, Herve Cassard<fnref n="1" />, Ben Hause<fnref n="4" />, Sarah Maman<fnref n="5" />, Gilles Meyer<fnref n="1" /> and Mariette F. Ducatez 1,*

Received: 20 April 2020; Accepted: 11 May 2020; Published: 13 May 2020

<figure><img src="image_3.png" /></figure>

check for
updates

**Abstract:** Bovine coronavirus (BCoV) is widespread in cattle and wild ruminant populations throughout the world. The virus causes neonatal calf diarrhea and winter dysentery in adult cattle, as well as upper and lower respiratory tract infection in young cattle. We isolated and deep sequenced whole genomes of BCoV from calves with respiratory distress in the south-west of France and conducted a comparative genome analysis using globally collected BCoV sequences to provide insights into the genomic characteristics, evolutionary origins, and global diversity of BCoV. Molecular clock analyses allowed us to estimate that the BCoV ancestor emerged in the 1940s, and that two geographically distinct lineages diverged from the 1960s-1970s. A recombination event in the spike gene (breakpoint at nt 1100) may be at the origin of the genetic divergence sixty years ago. Little evidence of genetic mixing between the spatially segregated lineages was found, suggesting that BCoV genetic diversity is a result of a global transmission pathway that occurred during the last century. However, we found variation in evolution rates between the European and non-European lineages indicating differences in virus ecology.

**Keywords:** bovine coronavirus; France; recombination; geographic clustering; molecular clock

## 1. Introduction

Bovine coronavirus (BCoV) is widespread in cattle and wild ruminant populations, resulting in neonatal calf diarrhea and winter dysentery in adult cattle, as well as upper and lower respiratory tract infection. BCoV is associated with the bovine respiratory disease complex (BRD), a leading cause of morbidity of young cattle worldwide with large economic costs due to mortality, treatment, and impeded growth performances, and pose a threat to public health as a result of antibiotic use. Having a multifactorial etiology, BRD is a complex disease triggered by the presence of one or several viruses and/or bacteria [1] favored by an altered state of the host immune system [2] and usually disturbed environmental factors [3] and several physical and biological stressors [2,4,5]. Viruses associated with BRD are capable of initiating a cascade of events such as immune suppression [2,6], respiratory

<footnote n="1">1 IHAP, Université de Toulouse, INRAE, ENVT, 31076 Toulouse, France; elias.salem@envt.fr (E.S.); herve.cassard@envt.fr (H.C.); gilles.meyer@envt.fr (G.M.)</footnote>

<footnote n="2">2 Biomedicine Discovery Institute and Department of Microbiology, Monash University, Melbourne, VIC 3800, Australia; vijay.dhanasekaran@monash.edu</footnote>

<footnote n="3">3 Program in Emerging Infectious Diseases, Duke-NUS Medical School, Singapore 169857, Singapore</footnote>

<footnote n="4">4 Diagnostic Medicine/Pathobiology, College of Veterinary Medicine, Kansas State University, Manhattan Kansas, KS 66506-5802, USA; bhause@cambridgetechnologies.net</footnote>

<footnote n="5">5 SIGENAE, GenPhySE, Université de Toulouse, INRAE, INPT, ENVT, 31326 Castanet Tolosan, France; sarah.maman@inrae.fr</footnote>

<footnote>* Correspondence: mariette.ducatez@envt.fr</footnote>
epithelium damage, [7] and altering commensal microbiota and local biofilm [8–10], thus leading to
complex secondary bacterial infections [8,10,11].

BCoV is a pneumoenteric virus frequently isolated from the fecal samples of diarrheic calves [12–15]
and from upper and lower respiratory tract samples during BRD episodes worldwide [16–19]. While both
forms of BCoV infection are widespread with high prevalence among cattle [12,13,15,20,21] and wild
ruminants [22], the putative differences between strains isolated from enteric and respiratory tracts
are still not clear, and it is unknown whether these differences are related to tissue tropism or simply
to distinct times and locations of isolation. Some reports indeed suggest that all the BCoV are similar
at genomic and antigenic levels [18,23], while others suggest that enteric and respiratory strains are
genetically and antigenically different [24,25]. Zhang et al. highlighted molecular differences within
a single host, with intra-host quasispecies and enteric strains more prone to genetic changes than
their respiratory counterparts [26]. The global genetic and antigenic diversity of BCoV is also poorly
understood. Using a short variable region of the BCoV genome and with sporadically collected samples
from a few countries in Europe, Asia, and North America, it has been suggested that BCoV are segregated
into European and American BCoV lineages [17,27–29], with periodic introductions of North American
BCoV observed in Asian countries (for example in Japan during the 1990s [30]) highlighting global
movement of BCoV likely through cattle trade.

Coronaviruses are enveloped particles with a large non-segmented, linear, polyadenylated, single-stranded positive sense RNA genome [31]. They belong to the family Coronaviridae in the order Nidovirales within the Coronavirinae subfamily [32]. Coronaviruses are classified into four genera [33]: *Alphacoronavirus*, *Betacoronavirus*, *Gammacoronavirus*, and *Deltacoronavirus*, and are known to infect birds and many mammals including humans [34]. BCoV are members of the *Betacoronavirus* genus within the *Coronaviridae* family, which includes several important human pathogens of zoonotic origin, such as the human coronavirus (HCoV) OC43, HKU1, MERS, SARS-CoV-1 [35], and SARS-CoV-2 [36,37], and mammalian viruses including equine coronaviruses (EqCoV) [38], canine respiratory coronaviruses (CRCoV) [39], swine hemagglutinating encephalomyelitis virus (HEV) [40], murine coronaviruses [32], and bat coronaviruses (BatCoV) HKU4, HKU5, and HKU9 [41–44]. BCoV are most closely related to members of the *Betacoronavirus1* (*Embecovirus*) species PHEV, OC43, and EqCoV and share similarity in genome organization [32].

Coronaviruses pose a constant public health threat, due to their low replication fidelity
and high genetic variability, making them prone to a constant changing pattern often leading to
emerging diseases [45]. The intensive investigation that followed the emergence of SARS-CoV1 [46]
and MERS-CoV [47] confirmed the zoonotic origins of those beta-coronaviruses [35,43,48–50],
their efficient adaptability [51], and their potency to switch species [52]. SARS-CoV2 was also
shown to have emerged from an animal reservoir and very successfully jumped into human causing
the current pandemics [36,37]. Based on antigenic similarity it is thought that the human OC43
and canine respiratory coronaviruses emerged from BCoV [53,54].

The aim of the present study was to characterize BCoV from both upper and lower respiratory tracts of symptomatic animals in the south-west of France. We aimed to isolate the viruses, sequence their full genomes by deep sequencing, and analyze the obtained genetic data using phylogeny and time scale evolution models to better understand the virus evolution. While our starting point was French field BCoV, our study aims at a better understanding of BCoV genetic evolution at a global scale, looking both at the genetic characteristics of viruses and at the collection times in order to contribute to discussions on BCoV origin (in time and space) and global spread.

**2. Materials and Methods**

All methods were carried out in accordance with relevant guidelines and regulations.
All experimental protocols were approved by the Ecole Nationale Vétérinaire de Toulouse and/or
the Genotoul bioinformatics platform Toulouse Midi-Pyrénées.
## 2.1. Sample Collection and Processing

Deep nasal swabs (NS) and bronchoalveolar lavages (BALs) were collected from herds affected with respiratory disease in the south-west of France during 2014. The calves included in this study were not treated nor vaccinated prior to sampling. The inclusion criteria were: (i) calves from breeding farms of maximum five months of age; (ii) with acute signs of pneumonia during BRD episodes and clinical signs for maximum four days; (iii) minimum 60% of the herd with clinical signs; (iv) clinical signs: hyperthermia (>40 °C) or hyperthermia (>39.5 °C) with cough, snuffles, or tachypnea; (v) respiratory distress and abnormal lung sounds (rhonchi). Table 1 summarizes the details on sampled herds. Four calves were sampled from each of the affected herds, and samples were pooled separately by tissue type (NS or BAL). Swabs were inserted deep in the calf nostril (approximatively 20 cm deep) and rotated for 30 s against the nasal mucosae. Swabs were placed in PBS, vortexed, and kept on ice for no more than 2 h until being stored at -80 °C. BALs by bronchoscopy were performed under deep anesthesia: 10 mg/kg Diazepam and 10 mg/kg ketamine were injected intravenously, and a sterile flexible bronchoscope (Olympus, Tokyo, Japan) was passed through the pharynx and visually guided to the first bronchial division corresponding to the caudal segment of the right apical lung lobe. The lung lobe was then lavaged with 100 mL of sterile isotonic saline solution and the liquid was immediately aspirated after the infusion and refrigerated until further processing.

<table><caption>**Table 1.** Presence of co-infecting pathogens in bronchoalveolar lavages (BAL) or nasal swabs (NS) pools as determined by real time using the vetmax screening kit.</caption><thead><tr><th colspan="2">Sample</th><th colspan="6">Ct Value<sup>2</sup></th></tr><tr><th>Id</th><th>Type<sup>1</sup></th><th>P. multocida</th><th>M. haemolytica</th><th>M. bovis</th><th>H. somni</th><th>BCoV</th><th>RSV</th><th>PI-3</th></tr></thead><tbody><tr><td>ICSA-4</td><td>BAL</td><td>26.75</td><td></td><td></td><td></td><td>30.17</td><td>22.74</td><td></td></tr><tr><td>ICSA-16</td><td>BAL</td><td></td><td></td><td></td><td>39.01</td><td>24.19</td><td></td><td></td></tr><tr><td>ICSA-17</td><td>BAL</td><td>27.04</td><td></td><td></td><td></td><td>38.9</td><td></td><td></td></tr><tr><td>ICSA 21</td><td>BAL</td><td>29.32</td><td>27.31</td><td>29.27</td><td></td><td>38.61</td><td></td><td></td></tr><tr><td>ICSA-4</td><td>NS</td><td>32.1</td><td></td><td></td><td></td><td>27.9</td><td>29.4</td><td>32.61</td></tr><tr><td>ICSA-16</td><td>NS</td><td></td><td></td><td></td><td></td><td>18.67</td><td></td><td></td></tr><tr><td>ICSA-17</td><td>NS</td><td>28.9</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ICSA 21</td><td>NS</td><td>24.4</td><td>26.7</td><td></td><td></td><td>34.93</td><td></td><td></td></tr></tbody></table>

## 2.2. Sample Preparation and NGS Sequencing

Pools of both sample types were filtered by a syringe driven 0.45 µm filter, and the BALs were concentrated by centrifugation (100,000× g, 2 h, 4 °C). Three hundred microliters of PBS were added to the pellet and left at 4 °C overnight before suspension. All samples were then treated with a cocktail of nucleases for 1 h at 37 °C to eliminate the host nucleic acids (Exonuclease I (20U, Thermo Fischer, Waltham, MA, USA), Benzonase (25U, Merck chemical, Darmstadt, Germany), Turbo DNAse (2U, Thermo Fisher Scientific, Waltham, MA, USA), and RNase I (10U, Thermo Fischer, Waltham, MA, USA). Viral nucleic acids were purified using the QIAamp MinElute virus spin kit (Qiagen, Hilden, Germany). Sequence-independent, single primer amplification was carried out as described by Liais et al. [55]. Briefly, reverse transcription was performed using tagged random hexamers, first and second strands were then randomly amplified by PCR. PCR products were purified and libraries processed using the Truseq nano kit (Illumina, San Diego, CA, USA) for next-generation sequencing with Illumina Miseq v2 platform (San Diego, CA, USA) generating paired 250 bp reads. Bioinformatics analysis was carried out using an in-house bioinformatics pipeline. Briefly, sequence reads from each library were filtered, demultiplexed, sorted, and pair-end reads merged. Viral composition and relative abundance were determined using GAAS [56], the metagenomes were compared to NCBI RefSeq virus database using Blastx. BCoV reads were also mapped against the Mebus BCoV genome (accession U00735.2) using the Burrows-Wheeler Alignment Tool. Reads of BCoV were

<footnote n="1">¹ BAL, bronchoalveolar lavage and NS, nasal swabs. ² Ct values >40 are not shown. P.: *Pasteurella*; M.: *Mycoplasma*; H.: *Histophilus*; BCoV: bovine coronavirus; RSV: bovine respiratory syncytial virus; PI-3: bovine parainfluenza-3 virus.</footnote>
as well assembled using CLC Genomics Workbench (Qiagen, Hilden, Germany). Specific primers were designed based on the extracted reads for RT PCR to fill sequence gaps between the contigs and the produced amplicons were sequenced on a 3130XL Applied Biosystems capillary sequencer (Applied Biosystems, Foster City, CA, USA).

## 2.3. Screening for Bovine Respiratory Pathogens

Samples were screened using the VetMAX Screening pack-Ruminant Respiratory Pathogens real-time PCR kit (Thermo Fischer, Waltham, MA, USA) according to the manufacturer instructions for the detection of *Mycoplasma bovis*, *Histophilus somni*, *Pasteurella multocida*, *Mannheimia haemolytica*, bovine coronavirus, bovine respiratory syncytial virus (BRSV), and bovine parainfluenza-3 virus (BPI-3).

## 2.4. BCoV Isolation

Virus isolation was attempted in Human Rectal Tumor-18G (HRT-18G) cells as follows: samples filtered by a syringe driven 0.22 µM filter were propagated on HRT-18G (ATCC, Manassas, VA, USA) cells in the presence of TPCK trypsin (1 µg/mL, Thermo Fisher Scientific, Waltham, MA, USA) in optiMEM Reduced Serum Media (Thermo Fisher Scientific, Waltham, MA, USA) supplemented with Penicillin (100 units/mL, Waltham, MA, USA), Streptomycin (10 mg/mL; Life technologies Carlsbad, CA, USA), Ciprofloxacin (10 µg/mL, Sigma-Aldrich, St Louis, MO, USA), Amphotericin B (2.5 µg/mL, Sigma-Aldrich, St Louis, MO, USA), and BM-cyclin (15 µg/mL, Sigma-Aldrich, St Louis, MO, USA). The cells were incubated for four days at 37 °C before pooling and filtration of the supernatant.

## 2.5. Virus sequence Datasets and Phylogenetic Analysis

Viruses sequenced in this study were combined with BCoV nucleotide sequence data available in GenBank, along with their sampling dates, host, and sample type (respiratory or fecal) (Supplementary Tables S1 and S2). We focused on the spike (S) and nucleocapsid (N) genes as the large majority of the available BCoV sequence data was on S and N genes. Multiple sequence alignments were generated with Clustal Omega [57,58] and manually optimized with BioEdit v7.1 [59]. Maximum likelihood (ML) trees were constructed in MEGA v6 [60] using the best-fit nucleotide substitution models: Tamura-Nei model with gamma-distributed rate variation (TN93+G) for N gene and the general time reversible model with gamma-distributed rate variation and a proportion of invariant sites (GTR+G+I) for the S gene. The quality of nucleotide sequences along with their temporal signal were assessed using a regression of the root-to-tip distances of the ML tree and virus sampling dates in TempEst v1.5 [61]. Samples that deviated from the regression were assumed to either be erroneous or incorrectly labelled and were removed from subsequent analysis.

## 2.6. Recombination Detection

Recombination detection was conducted for the N and S gene datasets using the Genetic Algorithm of Recombination Detection (GARD) method [62] in the Datamonkey server of HyPhy v2 [62,63]. Datasets of recombination-free sequences were constructed by dividing the original alignments at the detected breakpoints and were used in all subsequent analyses.

## 2.7. Molecular Dating

Nucleotide substitution rates and the evolutionary time scale of divergence of BCoV lineages were estimated using the year of virus sampling used as tip-calibrations in a relaxed molecular clock method under the Bayesian Markov chain Monte Carlo (MCMC) framework in BEAST v1.7.1 [64], implemented on a Galaxy workbench (http://galaxy-workbench.toulouse.inra.fr/). Molecular clock branch rates were derived from a log-normal distribution [65]. The Hasegawa-Kishino-Yano (HKY) nucleotide substitution model was specified separately for codon positions 1 + 2 and 3 and used
a constant population size coalescent model. Analysis were run for $2 \times 10^8$ generations, sampling every 10,000 generations. Defaults were used for all priors except for the clock rate for which a uniform prior distribution was specified with an initial value of $10^{-4}$ substitutions/nucleotide/year, based on previous estimates of coronavirus substitution rates [66]. Convergence of each of the parameters with a minimum effective sample size of 500, following a burning of 10%, was assessed using Tracer [64]. Maximum clade credible trees with the mean time of the most recent common ancestor (tMRCA) and their 95% highest posterior density (HPDs) were generated using the TreeAnnotator program included in the BEAST package and visualized in FigTree v1.4.3.

## 2.8. *Selective Pressure Analysis*

Selection pressure acting on the S gene were investigated by estimating the ratio between the mean number of non-synonymous (dN) and synonymous (dS) nucleotide substitution site (dN/dS) using the Single Likelihood Ancestor Counting (SLAC), Fixed Effects likelihood (FEL) [67,68] and the Mixed Effects Model of Evolution (MEME) method implemented in the datamonkey server of HyPhy v2 [63].

# 3. Results

## 3.1. *Genomic Surveillance of BCoV in France*

Upper (NS) and lower (BAL) respiratory tract samples were collected from 16 calves exhibiting acute respiratory disease during 2014 in the south-west of France (two respiratory samples per calf, four samples per pool; Table 1). Screening for other respiratory pathogens using real-time PCR resulted in the identification of co-infecting agents (Table 1), including *P. multocida*, *M. haemolytica*, *M. bovis*, and *H. somni*. All samples had at least one co-infecting pathogen, with *P. multocida* in 3/4 calves, followed by *M. haemolytica*, BRSV, *M. bovis*, BPI-3, and *H. somni* (Table 1). This suggests that the respiratory clinical signs were not due to BCoV alone but to co-infecting pathogens. A genetic focus was however made on BCoV in the present study. A greater proportion of BCoV were detected in NS samples (with lower Ct values) compared to BAL in three out of four pools studied.

Nucleic acids extracted from purified samples were subjected to NGS sequencing, focusing on BCoV genetic material arising separately from the upper and lower respiratory samples of infected cattle. In total, five complete and two partial BCoV genomes were generated (accession numbers MG757138-44), from which the nucleocapsid (N) gene and spike (S) gene sequences were extracted and utilized for subsequent evolutionary analysis (Supplementary Tables S1 and S2). Attempts of virus isolation in HRT-18G cells resulted in the successful isolation of one strain, BCoV/France/ICSA17_LBA/2014, collected from the BAL sample of infected calf in 2014.

## 3.2. *Recombination of BCoV*

Since coronavirus genomes are known to undergo recombination, we inferred the occurrence of recombinant breakpoints along the N and S gene alignments of BCoV collected since the 1960s, and found one recombination breakpoint in the S gene at nucleotide 1101 with significant differences in evolutionary histories (KH test) between the two loci at nucleotide positions 1–1100 and 1101–4089, whereas no breakpoints were detected in the N gene. These results suggest that a single phylogenetic tree may not accurately describe the evolutionary history of S gene and data were therefore partitioned for subsequent phylogenetic analysis.

## 3.3. *Phylogenetic Relationships and Spatio-Temporal Evolution of BCoV*

Phylogenetic analysis was conducted using 95 BCoV S gene sequences including seven from the present study, and 54 BCoV N gene sequences including five newly generated sequences. Our sequence dataset was predominantly restricted to sequences originating from North America, Asia and Europe highlighting the paucity of BCoV surveillance and the lack of well curated sequence datasets from around the globe. S gene alignments were partitioned at the detected recombination
breakpoint to create S gene datasets ranging from nucleotide positions 1–1100 and 1101–4089. In both data partitions, BCoV circulating since the 1990s formed two broad spatially segregated phylogenetic groups in Asia and Europe, respectively, derived from Mebus/1972-like strain collected in North America (Figures 1 and 2). There were several topological differences between the S gene data partitions that occurred both during the emergence of lineages in Europe and Asia, but we also identified topological differences within the Asian and European clade suggesting continued recombination among co-circulating lineages following spatial segregation (Figures 1 and 2).

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Dated phylogeny of BCoV S gene partition-1 (nt 1–1100). Tree generated using the relaxed clock model and a constant coalescent tree prior. Nodes correspond to mean tMRCAs. 95% highest</figcaption></figure>
posterior density (HPD) are indicated on nodes with grey boxes. French sequences (including those generated in the present study), European, Asian, and North American sequences are in red, orange, blue, and green fonts, respectively.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Dated phylogeny of BCoV S gene partition-2 (nt 1101–4089). Tree generated using the relaxed clock model and a constant coalescent tree prior. Nodes correspond to mean tMRCAs. 95% HPD are indicated on nodes with grey boxes. French sequences (including those generated in the present study), European, Asian, and North American sequences are in red, orange, blue, and green fonts, respectively.</figcaption></figure>
The Asian lineage of BCoV was predominated by strains collected in South Korea during 2002–2010 and both S gene partition phylogenies suggest that these strains were derived from strains similar to those reported from North America during 1994–1998. Estimation of tMRCA of the South Korean lineages showed that the mean tMRCA of S gene loci-1 was significantly more recent during 1995 (mean tMRCA: 1995; highest posterior density (HPD): 1992.5–1998.7) than in data partition-2 which constituted two distinct lineages that diverged during 1986 (mean tMRCA: 1986; HPD: 1980–1991), however the tMRCAs of the two co-circulating lineages overlapped with the tMRCA of partition-1 highlighting the co-circulation of multiple lineages in South Korea (Table 2). In addition, two strains collected recently from cattle in China exhibited differences in origin between the S data partitions: their first 1000 nucleotides appear to be derived independently from strains similar to those collected in North America during 1994–1998, whereas the second but longer partition was derived directly from South Korean specimens collected during 2002–2010. These results suggest historical movement of BCoV between the North America and Eurasia landmasses, as well as between Asian countries.

<table><caption>**Table 2.** Evolution rates and TMRCA for N and S gene partitions.</caption><thead><tr><th>Data</th><th>Gene</th><th>Demographic Model <sup>a</sup></th><th>Mean nt Sub Rate (95% HPD) × 10<sup>-3</sup> <sup>b</sup></th><th>Mean tMRCA (95% HPD) in Years before 2015</th></tr></thead><tbody><tr><td rowspan="3">Global</td><td>S partition 1</td><td>Constant</td><td>0.8 (0.5–1)</td><td>63.5 (53.9–75.9)</td></tr><tr><td>S partition 2</td><td>Constant</td><td>0.1 (0.7–1)</td><td>60.94 (51.7–70.4)</td></tr><tr><td>N</td><td>Constant</td><td>0.2 (0.1–0.4)</td><td>77.9 (46.2–123.2)</td></tr><tr><td rowspan="8">Europe</td><td>S partition 1</td><td>Constant</td><td>0.7 (0.5–0.9)</td><td>41.6 (37.5–53.5)</td></tr><tr><td rowspan="4">S partition 2</td><td>Bayesian skyline</td><td>0.7 (0.5–0.9)</td><td>41.9 (36.4–52.2)</td></tr><tr><td>Constant</td><td>0.8 (0.6–0.9)</td><td>39.7 (35.3–54.4)</td></tr><tr><td>Bayesian skyline</td><td>0.8 (0.6–1.0)</td><td>37.5 (35.0–47.9)</td></tr><tr><td>Bayesian skyline</td><td>0.5 (0.3–0.7)</td><td>35.8 (35.0–45.4)</td></tr><tr><td rowspan="2">N</td><td>BSP</td><td>0.5 (0.3–0.7)</td><td>35.6 (35.0–46.0)</td></tr><tr><td></td><td></td><td></td><td></td></tr></tbody></table>

The European lineage (EU-lineage) was represented by BCoV collected from several countries in Europe since the 1990s, including France, Denmark, Germany, Italy, and Sweden; however, excluding France, most BCoV sequence data in Europe derives from 2002–2010. There is also evidence that the EU-lineage evolved through recombination as the ancestral relationship of these lineages to M80844/Geissen-Germany/1989 differed between the two loci (Figures 1 and 2), indicating that the recombination event (breakpoint at nt position 1100) may be at the origin of the distinct Asian and European lineages.

BCoV sequences collected from calves in France between 2013–2014 revealed a greater genetic heterogeneity in S gene loci-1 (1–1000 nt) with evidence of co-circulation of two lineages that diverged as early as 2000 (mean tMRCA: Aug 2000; 95% HPD: 1998.5–2009.8; Figure 1), whereas loci-2 exhibited fewer differences between contemporary viruses and revealed a mean divergence of approximately four years (mean tMRCA: Mar 2009; 95% HPD: 2007.5–2010.8). When looking specifically at the south-western French BCoV sequences, BCoV/France/ICSA21L3/2014 differed from its counterparts on locus 1 (being closely related to a recent north-western French sequence, BCoV/France/EPI-Caen/2013) but clustered with its south-western counterparts for the locus 2 (Figure 2). These results highlight a complex evolutionary history that is determined by migration and pervasive recombination.

Our N gene dataset comprised lower sampling outside of France (Table S2). All N gene sequences of BCoV strains from France, isolated since 2003, formed a monophyletic group, with the mean tMRCA of the contemporary BCoV (2013–2014) isolates estimated during 2004 (95% HPD: 1999.2–2009.2), indicating a result that overlaps with the tMRCA of both S gene partitions.

We estimated the evolutionary rate for the S and N gene datasets using the relaxed clock model, and found a relatively lower rates of nucleotide substitution in BCoV datasets, ranging from a mean rate of 0.81 × 10⁻³ (95% HPD: 0.5–0.8 × 10⁻³) for S gene data partition-1, 1.0 × 10⁻³ (HPD: 0.7–1.3 × 10⁻³)

<footnote><sup>a</sup> Bayesian skyline model is a flexible coalescent prior that allows variation in population size over time. <sup>b</sup> estimated using the uncorrelated log-normal relaxed clock model and a constant population size model.</footnote>
for S gene data partition-2 and a lower but with wide-confidence intervals for the N gene dataset (mean rate: $0.2 \times 10^{-3}$; 95% HPD: 0.1–0.4 $\times 10^{-3}$). These estimates were similar to those estimated for other coronaviruses, for example the mean nucleotide substitution rate of MERS was estimated at $1.1 \times 10^{-3}$ (95% CI: 0.8–1.4 $\times 10^{-3}$) [69]. Estimates of the coefficient of variance of the relaxed clock model ranged from a mean of 0.61 for the S gene partition-1, to 0.7 for S partition-2 and the N gene datasets, and HPD intervals of all datasets were lower than 1, indicating that there was significant variation in nucleotide substitution rates along the BCoV branches, and also suggesting that the relaxed clock model was more appropriate than using a strict clock model.

The tMRCA of globally sampled BCoV for the N and S genes was estimated during 1940 and 1964, suggesting the common ancestor for known BCoV was about 51 to 75 years prior to the last included sequence date (2015). Whereas the most common recent ancestor of the EU-lineage was 1978–1981 (Figures 1 and 2). The recombination event that likely led to the two main lineages with clear geographic distinction (US/Asia and European lineages) co-circulating today may thus have occurred in the 1960s–1970s.

### 3.4. Selection Pressure

The selection analysis for the S gene was conducted using two partitions. The average ratio between the mean number of non-synonymous (dN) and synonymous (dS) nucleotide substitution sites (dN/dS) was estimated to be 0.310 for the S gene, with 28 positively selected amino acids (from SLAC, FEL and MEME; Table 3). Amino acids 35, 113, 115, 447, 499, 501, 525, 716, 893, and 1239 were detected as positively selected in at least two of the three used models (Table 3).

<table><caption>**Table 3.** Amino acid sites detected to be under positive selection (with *p*-values < 0.1) through either of the selection methods (SLAC, FEL, and MEME).</caption><thead><tr><th>Codon</th><th>SLAC p-Value</th><th>FEL p-Value</th><th>MEME p-Value</th></tr></thead><tbody><tr><td>35</td><td>0.0600</td><td>0.0700</td><td><strong>0.0130</strong></td></tr><tr><td><u>113</u></td><td><strong>0.0340</strong></td><td><strong>0.0220</strong></td><td>0.0700</td></tr><tr><td><u>115</u></td><td>0.0670</td><td><strong>0.0210</strong></td><td>#N/A</td></tr><tr><td>155</td><td>#N/A</td><td>#N/A</td><td>0.0500</td></tr><tr><td>179</td><td>#N/A</td><td>#N/A</td><td>0.0680</td></tr><tr><td>257</td><td>#N/A</td><td>0.0820</td><td>#N/A</td></tr><tr><td>304</td><td>#N/A</td><td>#N/A</td><td><strong>0.0160</strong></td></tr><tr><td><u>447</u></td><td><strong>0.0260</strong></td><td><strong>0.0260</strong></td><td>0.0660</td></tr><tr><td>484</td><td>#N/A</td><td>#N/A</td><td>0.0980</td></tr><tr><td><u>499</u></td><td><strong>0.0230</strong></td><td><strong>0.0090</strong></td><td>#N/A</td></tr><tr><td><u>501</u></td><td><strong>0.0010</strong></td><td><strong>0.0040</strong></td><td><strong>0.0060</strong></td></tr><tr><td>509</td><td>#N/A</td><td><strong>0.0410</strong></td><td>#N/A</td></tr><tr><td><u>525</u></td><td><strong>0.0220</strong></td><td><strong>0.0280</strong></td><td>0.0910</td></tr><tr><td>543</td><td>#N/A</td><td>#N/A</td><td>0.0760</td></tr><tr><td>689</td><td>#N/A</td><td>#N/A</td><td><strong>0.0380</strong></td></tr><tr><td>716</td><td>0.0590</td><td><strong>0.0490</strong></td><td>#N/A</td></tr><tr><td><u>744</u></td><td>#N/A</td><td>#N/A</td><td><strong>0.0060</strong></td></tr><tr><td>767</td><td>#N/A</td><td>#N/A</td><td><strong>0.0000</strong></td></tr><tr><td>769</td><td>#N/A</td><td>#N/A</td><td>0.0940</td></tr><tr><td>892</td><td>#N/A</td><td>#N/A</td><td>0.0710</td></tr><tr><td><u>893</u></td><td>#N/A</td><td><strong>0.0410</strong></td><td><strong>0.0010</strong></td></tr><tr><td>1085</td><td>#N/A</td><td>#N/A</td><td><strong>0.0030</strong></td></tr><tr><td>1188</td><td>#N/A</td><td>#N/A</td><td>0.0040</td></tr><tr><td>1206</td><td>#N/A</td><td><strong>0.0450</strong></td><td>#N/A</td></tr><tr><td>1237</td><td>#N/A</td><td>#N/A</td><td><strong>0.0380</strong></td></tr><tr><td>1239</td><td>#N/A</td><td><strong>0.0110</strong></td><td>0.0960</td></tr><tr><td>1269</td><td>#N/A</td><td>#N/A</td><td><strong>0.0000</strong></td></tr><tr><td>1362</td><td>#N/A</td><td><strong>0.0150</strong></td><td>#N/A</td></tr></tbody></table>

Sites with *p*-values < 0.05 are shown in bold. Codons positively selected as shown by at least two of the three selection methods are underlined. #N/A: no positive selection detected.
# 4. Discussion

## 4.1. Prevalence of BCoV in Europe

Bovine coronavirus is considered widespread in Europe and worldwide in both its enteric and respiratory forms. Furthermore, BCoV is considered a main pathogen associated with pneumonia in feedlots in North America [15,16,70,71]. Although BCoV is no less a very common pathogen in Europe, limited data is available on its epidemiology. Norwegian studies have revealed a high seroprevalence in cattle with 39.3% and 80.7% in calves and herds, respectively, in 2009, and 72.2% in herds tested in 2016 [21,72]. In Italy, four respiratory disease outbreaks during 2006 were associated with BCoV without concurrent detection of other known respiratory pathogen but with the detection of BCoV in the intestines in two of the four outbreaks [19]. In Belgium, seroconversion for BCoV was detected in 30% of veal calves herds [73]. In a retrospective study in Ireland (2008–2012), BCoV was the most frequently detected virus in cattle nasal swabs from respiratory disease outbreaks [74]. In France, BCoV was detected in 17%–63% of respiratory samples collected in two independent studies [75]. Taken together, literature on BCoV in Europe is limited but highlights the high circulation of the pathogen at the continental scale.

## 4.2. No obvious Link between Virus Tropism and Genetic Markers

It is so far unknown whether a unique or distinct BCoV are responsible for respiratory and enteric infections [18,23,25,76]. As recently observed in Austria and Slovakia [23], we did not succeed in identifying genetic determinants of tissue tropism, nor between enteric and respiratory BCoV strains (based on published sequences), nor between strains present in the upper and lower respiratory tracts (based on the sequences generated within the framework of the present study). In the absence of viral genetic determinants of tissue tropism, the immune status of the animals at the time of the BCoV infection, infectious dose, and route of inoculation (oro-fecal versus aerosol) may play a role in determining sites of infection [77]. In addition, the animals sampled in the present study suffered from respiratory disease, but as both BCoV and co-infecting pathogens were detected, a direct link between BCoV and respiratory signs cannot be made. Further studies are clearly warranted to understand the mechanisms behind the differential tissue tropism of BCoV.

## 4.3. Little Exchanges of BCoV between Europe and the Other Continents

Among available BCoV complete genomes, a minority are of European origins. In this study we therefore characterized respiratory strains of BCoV, from both upper and lower respiratory tracts, at the genetic level. Despite recent efforts, there are notable gaps in our understanding of bovine coronavirus dynamics and diversity. Our results show that the French strains described here cluster with all other Europeans strains, thus indicating a homogeneous evolution of BCoV in Europe, whereas a greater mixing was observed between isolates collected in North America and Europe. This pattern of geographical and temporal distinction between strains was previously described for BRSV, for which distinct phylogenetic groups with continent-specific strains were identified [78–80]. Similarly, BVDV also harbors a geographic clustering [81]. This peculiar evolutionary process could be explained by the relative isolation of Europe as far as cattle traffic is concerned. According to FAOstat [82], cattle traffic involving Europe is mainly unidirectional with export of live animals from France (for the most part) to other continents with a negligible importation activity. On the other hand, the import/export activity between North America and East Asia is much more fluid. In addition, Kin et al. recently highlighted this important trade constraint on BCoV spread at a global scale [28]. It is however unclear why/how BCoV, BRSV, or BVDV would spread differently than other bovine RNA viruses. Influenza D viruses (IDV), which primarily infects bovine hosts, seem for example to have emerged recently but genetically and anti-genetically similar viruses circulate in Europe and North America [29,83–86]. One may therefore hypothesize that different routes of transmission may be used for the different pathogens. In addition, human coronaviruses (OC43 and 229E) were shown to
induce short-term immunity [87]. If BCoV led to similar responses in cattle, herds with animals of multiple ages (as seen more often in France than in North America) might face re-infections with BCoV and might maintain virus over long periods of time, which would impact the virus evolution.

Recombination likely drove the different evolutionary process taken by the Asian and North American strains and the Europeans strains; in the first partition of the S gene dataset, the strains M80844 Giessen/1989 clustered along with French strain D00731/1979, US strain AF058942 LY138 Utah/1965, South Koreans EU401988 SUN5/1994 and EU401987 A3/1994, and North Americans AF220295 Quebec/1972 and U00735 Mebus/1972, suggesting a common ancestor of this cluster. Whereas, in the partition-2 of the S dataset, we noticed a change in the phylogenetic scheme: the M80844 Giessen/1989 clustered with all Europeans strains, while the more recent North Korean (year 2002) and recent Chinese KM985631 HLJ/T4/2014 and KU886219 AKS/01/2015 clustered with the US strains of the 1990s, suggesting a common origin. The old BCoV strains (pre 1994) seemed to have the same origin irrespective of their geographic origin until an event of recombination led the European strains to evolve separately.

### 4.4. *Evolution of BCOV*

Using the different models (constant population under strict and relaxed clock for all the strains and the constant population and BSP under relaxed clock for the European strains), the mean substitution rates calculated for the spike genes ranged between 3.9 and $9.8 \times 10^{-4}$ nucleotide substitution per site per year, which is a very similar range as previously published ($2.1-10.5 \times 10^{-4}$ [88,89]). The slower evolution of N than of S gene is also in line with previous reports from the literature [89]. Substitution rates were in the same range for other coronaviruses with $0.6 \times 10^{-3}$, $0.7 \pm 2 \times 10^{-3}$ for HCoV-OC43 and TGEV, respectively [90,91]. We observed a slightly faster evolution for European BCoV than for the whole dataset. An important bias in our analyses comes from the lack of sequence data and a timing difference between European and American BCoV, such that most North American specimens were at least ten years older than European strains. In addition, the production systems (different numbers of animals per farms, feedlots in North America but not in Europe), the virus prevalence, and the vaccine coverages are likely different, therefore the selection pressures must be different. Our selection pressure calculations unfortunately did not allow for highlighting any difference in selection pressure between the continents. The 10 positively selected amino acids detected here (Table 3) were spread out over the whole spike protein and could not be linked to a biological function due to the limited research so far carried out on BCoV; no spike structure has so far been resolved (comparison is only possible with the SARS spike), and no clear antigenic sites have been detected to the best of our knowledge. Further studies are thus warranted to understand the role of these positions in BCoV evolution.

As far as the tMRCA is concerned, our analysis with the new sequence data confirmed the values previously estimated: 1978 (95% CI: 1910–1963) or 1944 (95% CI: 1910–1963) [88,89]. Here again, more sequences available throughout the last few decades would allow for an increased precision in the tMRCA. We showed here for the first time that a recombination event may be at the origin of the viruses that evolved separately in Europe on the one side and in America/Asia on the other side. This event likely occurred in the 1960s–1970s, when cattle exchanges between Europe and the other continents was already very limited.

### 5. Conclusions

Taken together, our results showed that the BCoV ancestor emerged in the 1940s, and that two geographically distinct lineages diverged from the 1960s–1970s with no genetic mixing detected thus far. Further studies are warranted to understand the mechanisms behind BCoV tissue tropism. Improved surveillance and sequencing are needed to precisely estimate the evolutionary and epidemiological parameters of BCoV circulation and to identify the recombinations in genomic evolution.
**Supplementary Materials:** The following are available online at http://www.mdpi.com/1999-4915/12/5/534/s1, Table S1: Details of the strains used for S gene sequences analyses, Table S2: Details of the strains used for N gene sequences analyses.

**Author Contributions:** E.S., G.M., and M.F.D. designed the study and interpreted the results. E.S., V.D., G.M., and M.F.D. drafted the manuscript. E.S., H.C., and G.M. collected the samples. E.S., V.D., B.H., S.M., and M.F.D. performed the analyses. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the 'RESPICARE' grant of the Institut Carnot Santé Animale (ICSA), and by the French Ministry of Agriculture. Elias Salem was supported by a PhD scholarship of the Lebanese University.

**Acknowledgments:** We are grateful to the Genotoul bioinformatics platform Toulouse Midi-Pyrénées and Sigenae group for providing help and/or computing and/or storage resources on the Galaxy interface http://sigenae-workbench.toulouse.inra.fr.

**Conflicts of Interest:** The present study was part of the PhD thesis of Elias Salem entitled "Infectious bovine bronchopneumonia: deciphering the complex microbiomes and evolutionary and clinical features of yet unfamiliar viruses" [92].

**References**

1. Grissett, G.P.; White, B.J.; Larson, R.L. Structured Literature Review of Responses of Cattle to Viral and Bacterial Pathogens Causing Bovine Respiratory Disease Complex. *J. Vet. Intern. Med.* **2015**, *29*, 770–780.
[CrossRef] [PubMed]

2. Caswell, J.L. Failure of respiratory defenses in the pathogenesis of bacterial pneumonia of cattle. *Vet. Pathol.* **2014**, *51*, 393–409. [CrossRef] [PubMed]

3. Hilton, W.M. BRD in 2014: Where have we been, where are we now, and where do we want to go? *Anim. Health Res. Rev.* **2014**, *15*, 120–122. [CrossRef]

4. Hodgson, P.D.; Aich, P.; Stookey, J.; Popowych, Y.; Potter, A.; Babiuk, L.; Griebel, P.J. Stress significantly increases mortality following a secondary bacterial respiratory infection. *Vet. Res.* **2012**, *43*, 21. [CrossRef]

5. Hodgson, P.D.; Aich, P.; Manuja, A.; Hokamp, K.; Roche, F.M.; Brinkman, F.S.L.; Potter, A.; Babiuk, L.A.; Griebel, P.J. Effect of stress on viral-bacterial synergy in bovine respiratory disease: Novel mechanisms to regulate inflammation. *Comp. Funct. Genomics* **2005**, *6*, 244–250. [CrossRef] [PubMed]

6. Risalde, M.A.; Molina, V.; Sánchez-Cordón, P.J.; Romero-Palomo, F.; Pedrera, M.; Garfia, B.; Gómez-Villamandos, J.C. Pathogenic mechanisms implicated in the intravascular coagulation in the lungs of BVDV-infected calves challenged with BHV-1. *Vet. Res.* **2013**, *44*, 1–13. [CrossRef]

7. Valarcher, J.-F.; Taylor, G. Bovine respiratory syncytial virus infection. *Vet. Res.* **2007**, *38*, 153–180. [CrossRef]

8. Boukahil, I.; Czuprynski, C.J. Mannheimia haemolytica biofilm formation on bovine respiratory epithelial cells. *Vet. Microbiol.* **2016**, *197*, 129–136. [CrossRef]

9. Elswaifi, S.F.; Scarratt, W.K.; Inzana, T.J. The role of lipooligosaccharide phosphorylcholine in colonization and pathogenesis of Histophilus somni in cattle. *Vet. Res.* **2012**, *43*. [CrossRef]

10. Sandal, I.; Shao, J.Q.; Annadata, S.; Apicella, M.A.; Boye, M.; Jensen, T.K.; Saunders, G.K.; Inzana, T.J. Histophilus somni biofilm formation in cardiopulmonary tissue of the bovine host following respiratory challenge. *Microbes Infect.* **2009**, *11*, 254–263. [CrossRef]

11. Zecchinon, L.; Fett, T.; Desmecht, D. How Mannheimia haemolytica defeats host defence through a kiss of death mechanism. *Vet. Res.* **2005**, *36*, 133–156. [CrossRef] [PubMed]

12. Gomez, D.E.; Arroyo, L.G.; Poljak, Z.; Viel, L.; Weese, J.S. Detection of Bovine Coronavirus in Healthy and Diarrheic Dairy Calves. *J. Vet. Intern. Med.* **2017**. [CrossRef] [PubMed]

13. Ribeiro, J.; Lorenzetti, E.; Alfieri, A.F.; Alfieri, A.A. Molecular detection of bovine coronavirus in a diarrhea outbreak in pasture-feeding Nellore steers in southern Brazil. *Trop. Anim. Health Prod.* **2016**, *48*, 649–653. [CrossRef] [PubMed]

14. Alkan, F.; Ozkul, A.; Bilge-Dagalp, S.; Karaoglu, T.; Oguzoglu, T.C.; Caliskan, E.; Burgu, I. The Detection and Genetic Characterization Based on the S1 Gene Region of BCoVs from Respiratory and Enteric Infections in Turkey. *Transbound. Emerg. Dis.* **2011**, *58*, 179–185. [CrossRef]

15. Hasoksuz, M.; Hoet, A.E.; Loerch, S.C.; Wittum, T.E.; Nielsen, P.R.; Saif, L.J. Detection of respiratory and enteric shedding of bovine coronaviruses in cattle in an Ohio feedlot. *J. Vet. Diagn. Investig.* **2002**, *14*, 308–313. [CrossRef]
16. Fulton, R.W.; Step, D.L.; Wahrmund, J.; Burge, L.J.; Payton, M.E.; Cook, B.J.; Burken, D.; Richards, C.J.; Confer, A.W. Bovine coronavirus (BCV) infections in transported commingled beef cattle and sole-source ranch calves. *Can. J. Vet. Res.* **2011**, *75*, 191–199.

17. Fulton, R.W.; Ridpath, J.F.; Burge, L.J. Bovine coronaviruses from the respiratory tract: Antigenic and genetic diversity. *Vaccine* **2013**, *31*, 886–892. [CrossRef]

18. Saif, L.J. Bovine respiratory coronavirus. *Vet. Clin. North Am.—Food Anim. Pract.* **2010**, *26*, 349–364. [CrossRef]

19. Decaro, N.; Campolo, M.; Desario, C.; Cirone, F.; D'Abramo, M.; Lorusso, E.; Greco, G.; Mari, V.; Colaianni, M.L.; Elia, G.; et al. Respiratory disease associated with bovine coronavirus infection in cattle herds in Southern Italy. *J. Vet. Diagnostic Investig.* **2008**, *20*, 28–32. [CrossRef]

20. Bok, M.; Miño, S.; Rodriguez, D.; Badaracco, A.; Nuñes, I.; Souza, S.P.; Bilbao, G.; Louge Uriarte, E.; Galarza, R.; Vega, C.; et al. Molecular and antigenic characterization of bovine Coronavirus circulating in Argentinean cattle during 1994-2010. *Vet. Microbiol.* **2015**, *181*, 221–229. [CrossRef]

21. Toftaker, I.; Sanchez, J.; Stokstad, M.; Nødtvedt, A. Bovine respiratory syncytial virus and bovine coronavirus antibodies in bulk tank milk - risk factors and spatial analysis. *Prev. Vet. Med.* **2016**, *133*, 73–83. [CrossRef] [PubMed]

22. Alekseev, K.P.; Vlasova, A.N.; Jung, K.; Hasoksuz, M.; Zhang, X.; Halpin, R.; Wang, S.; Ghedin, E.; Spiro, D.; Saif, L.J. Bovine-like coronaviruses isolated from four species of captive wild ruminants are homologous to bovine coronaviruses, based on complete genomic sequences. *J. Virol.* **2008**, *82*, 12422–12431. [CrossRef] [PubMed]

23. Vilček, S.; Jacková, A.; Kolesárová, M.; Vlasáková, M. Genetic variability of the S1 subunit of enteric and respiratory bovine coronavirus isolates. *Acta Virol.* **2017**, *61*, 212–216. [CrossRef] [PubMed]

24. Hasoksuz, M.; Lathrop, S.; Al-dubaib, M.A.; Lewis, P.; Saif, L.J. Antigenic variation among bovine enteric coronaviruses (BECV) and bovine respiratory coronaviruses (BRCV) detected using monoclonal antibodies. *Arch. Virol.* **1999**, *144*, 2441–2447. [CrossRef]

25. Gélinas, A.M.; Boutin, M.; Sasseville, A.M.J.; Dea, S. Bovine coronaviruses associated with enteric and respiratory diseases in Canadian dairy cattle display different reactivities to anti-HE monoclonal antibodies and distinct amino acid changes in their HE, S and ns4.9 protein. *Virus Res.* **2001**, *76*, 43–57. [CrossRef]

26. Zhang, X.; Hasoksuz, M.; Spiro, D.; Halpin, R.; Wang, S.; Vlasova, A.; Janies, D.; Jones, L.R.; Ghedin, E.; Saif, L.J. Quasispecies of bovine enteric and respiratory coronaviruses based on complete genome sequences and genetic changes after tissue culture adaptation. *Virology* **2007**, *363*, 1–10. [CrossRef]

27. Gunn, L.; Collins, P.J.; O'Connell, M.J.; O'Shea, H. Phylogenetic investigation of enteric bovine coronavirus in Ireland reveals partitioning between European and global strains. *Ir. Vet. J.* **2015**, *68*, 31. [CrossRef]

28. Kin, N.; Miszczak, F.; Diancourt, L.; Caro, V.; Moutou, F.; Vabret, A.; Ar Gouilh, M. Comparative molecular epidemiology of two closely related coronaviruses, bovine coronavirus (BCoV) and human coronavirus OC43 (HCoV-OC43), reveals a different evolutionary pattern. *Infect. Genet. Evol.* **2016**, *40*, 186–191. [CrossRef]

29. Suzuki, T.; Otake, Y.; Uchimoto, S.; Hasebe, A.; Goto, Y. Genomic characterization and phylogenetic classification of bovine coronaviruses through whole genome sequence analysis. *Viruses* **2020**, *12*, 183. [CrossRef]

30. Kanno, T.; Hatama, S.; Ishihara, R.; Uchida, I. Molecular analysis of the S glycoprotein gene of bovine coronaviruses isolated in Japan from 1999 to 2006. *J. Gen. Virol.* **2007**, *88*, 1218–1224. [CrossRef]

31. Vijgen, L.; Keyaerts, E.; Lemey, P.; Moës, E.; Li, S.; Vandamme, A.M.; Van Ranst, M. Circulation of genetically distinct contemporary human coronavirus OC43 strains. *Virology* **2005**, *337*, 85–92. [CrossRef] [PubMed]

32. Fields, M.N.; Masters, P.; Perlman, S. Coronaviridae. In *Fields Virology*, 6th ed.; Knipe, D.M., Ed.; Wolters Kluwer Health/Lippincott Williams & Wilkins: Philadelphia, PA, USA, 2013; Volume 1, pp. 825–857. ISBN 9781451105636.

33. Wertheim, J.O.; Chu, D.K.W.; Peiris, J.S.M.; Kosakovsky Pond, S.L.; Poon, L.L.M. A Case for the Ancient Origin of Coronaviruses. *J. Virol.* **2013**, *87*, 7039–7045. [CrossRef] [PubMed]

34. Fehr, A.R.; Perlman, S. Coronaviruses: An Overview of Their Replication and Pathogenesis. *Coronaviruses* **2015**, *1282*, 1–23.

35. Cui, J.; Li, F.; Shi, Z.L. Origin and evolution of pathogenic coronaviruses. *Nat. Rev. Microbiol.* **2019**, *17*, 181–192. [CrossRef] [PubMed]
36. Gralinski, L.E.; Menachery, V.D. Return of the coronavirus: 2019-nCoV. *Viruses* **2020**, *12*, 135. [CrossRef] [PubMed]

37. Andersen, K.G.; Rambaut, A.; Lipkin, W.I.; Holmes, E.C.; Garry, R.F. The proximal origin of SARS-CoV-2. *Nat. Med.* **2020**, *89*, 44–48. [CrossRef]

38. Nemoto, M.; Kanno, T.; Bannai, H.; Tsujimura, K.; Yamanaka, T.; Kokado, H. Antibody response to equine coronavirus in horses inoculated with a bovine coronavirus vaccine. *J. Vet. Med. Sci.* **2017**, *79*, 1889–1891. [CrossRef]

39. Decaro, N.; Buonavoglia, C. Canine Coronavirus: Not Only an Enteric Pathogen. *Vet. Clin. North Am.—Small Anim. Pract.* **2011**, *41*, 1121–1132. [CrossRef]

40. Dong, B.; Lu, H.; Zhao, K.; Liu, W.; Gao, W.; Lan, Y.; Zhao, J.; Tang, B.; Song, D.; He, W.; et al. Identification and genetic characterization of porcine hemagglutinating encephalomyelitis virus from domestic piglets in China. *Arch. Virol.* **2014**, *159*, 2329–2337. [CrossRef]

41. Van Boheemen, S.; de Graaf, M.; Lauber, C.; Bestebroer, T.M.; Raj, V.S.; Zaki, A.M.; Osterhaus, A.D.M.E.; Haagmans, B.L.; Gorbalenya, A.E.; Snijder, E.J.; et al. Genomic characterization of a newly discovered coronavirus associated with acute respiratory distress syndrome in humans. *MBio* **2012**, *3*, e00473-12. [CrossRef]

42. Lau, S.K.P.; Li, K.S.M.; Tsang, A.K.L.; Lam, C.S.F.; Ahmed, S.; Chen, H.; Chan, K.-H.; Woo, P.C.Y.; Yuen, K.-Y. Genetic Characterization of Betacoronavirus Lineage C Viruses in Bats Reveals Marked Sequence Divergence in the Spike Protein of Pipistrellus Bat Coronavirus HKU5 in Japanese Pipistrelle: Implications for the Origin of the Novel Middle East Respiratory Sy. *J. Virol.* **2013**, *87*, 8638–8650. [CrossRef] [PubMed]

43. Yang, Y.; Du, L.; Liu, C.; Wang, L.; Ma, C.; Tang, J.; Baric, R.S.; Jiang, S.; Li, F. Receptor usage and cell entry of bat coronavirus HKU4 provide insight into bat-to-human transmission of MERS coronavirus. *Proc. Natl. Acad. Sci. USA* **2014**, *111*, 12516–12521. [CrossRef] [PubMed]

44. Lau, S.K.P.; Poon, R.W.S.; Wong, B.H.L.; Wang, M.; Huang, Y.; Xu, H.; Guo, R.; Li, K.S.M.; Gao, K.; Chan, K.-H.; et al. Coexistence of Different Genotypes in the Same Bat and Serological Characterization of Rousettus Bat Coronavirus HKU9 Belonging to a Novel Betacoronavirus Subgroup. *J. Virol.* **2010**, *84*, 11385–11394. [CrossRef] [PubMed]

45. Denison, M.R.; Graham, R.L.; Donaldson, E.F.; Eckerle, L.D.; Baric, R.S. Coronaviruses. *RNA Biol.* **2011**, *8*, 270–279. [CrossRef] [PubMed]

46. Drosten, C. Identification of a novel coronavirus in patients with severe acute respiratory syndrome. *N. Engl. J. Med.* **2003**, *348*, 1967–1976. [CrossRef] [PubMed]

47. Zaki, A.M.; van Boheemen, S.; Bestebroer, T.M.; Osterhaus, A.D.M.E.; Fouchier, R.A.M. Isolation of a Novel Coronavirus from a Man with Pneumonia in Saudi Arabia. *N. Engl. J. Med.* **2012**, *367*, 1814–1820. [CrossRef]

48. Wang, Q.; Qi, J.; Yuan, Y.; Xuan, Y.; Han, P.; Wan, Y.; Ji, W.; Li, Y.; Wu, Y.; Wang, J.; et al. Bat origins of MERS-CoV supported by bat Coronavirus HKU4 usage of human receptor CD26. *Cell Host Microbe* **2014**, *16*, 328–337. [CrossRef]

49. Han, H.-J.; Yu, H.; Yu, X.-J. Evidence for zoonotic origins of Middle East respiratory syndrome coronavirus. *J. Gen. Virol.* **2016**, *97*, 274–280. [CrossRef]

50. Mohd, H.A.; Al-Tawfiq, J.A.; Memish, Z.A. Middle East Respiratory Syndrome Coronavirus (MERS-CoV) origin and animal reservoir. *Virol. J.* **2016**, *13*, 87. [CrossRef]

51. Lau, S.; Wong, A.; Lau, T.; Woo, P. Molecular Evolution of MERS Coronavirus: Dromedaries as a Recent Intermediate Host or Long-Time Animal Reservoir? *Int. J. Mol. Sci.* **2017**, *18*, 2138. [CrossRef]

52. Lu, G.; Wang, Q.; Gao, G.F. Bat-to-human: Spike features determining ‘host jump’ of coronaviruses SARS-CoV, MERS-CoV, and beyond. *Trends Microbiol.* **2015**, *23*, 1–11. [CrossRef] [PubMed]

53. Lefkowitz, E.J.; Dempsey, D.M.; Hendrickson, R.C.; Orton, R.J.; Siddell, S.G.; Smith, D.B. Virus taxonomy: The database of the International Committee on Taxonomy of Viruses (ICTV). *Nucleic Acids Res.* **2017**, *46*, D708–D717. [CrossRef] [PubMed]

54. Vijgen, L.; Keyaerts, E.; Moës, E.; Thoelen, I.; Lemey, P.; Vandamme, A.; Ranst, V.; Moe, E.; Wollants, E.; Ranst, M.V. Complete Genomic Sequence of Human Coronavirus OC43: Molecular Clock Analysis Suggests a Relatively Recent Zoonotic Coronavirus Transmission Event. *J. Virol.* **2005**, *79*, 1595–1604. [PubMed]

55. Liais, E.; Croville, G.; Mariette, J.; Delverdier, M.; Lucas, M.N.; Klopp, C.; Lluch, J.; Donnadieu, C.; Guy, J.S.; Corrand, L.; et al. Novel avian coronavirus and fulminating disease in Guinea Fowl, France. *Emerg. Infect. Dis.* **2014**, *20*, 105–108. [CrossRef]
56. Angly, F.E.; Willner, D.; Prieto-Davo, A.; Edwards, R.A.; Schmieder, R.; Vega-Thurber, R.; Antonopoulos, D.A.; Barott, K.; Cottrell, M.T.; Desnues, C.; et al. The GAAS metagenomic tool and its estimations of viral and microbial average genome size in four major biomes. *PLoS Comput. Biol.* **2009**, *5*, [CrossRef]

57. Goujon, M.; McWilliam, H.; Li, W.; Valentin, F.; Squizzato, S.; Paern, J.; Lopez, R. A new bioinformatics analysis tools framework at EMBL-EBI. *Nucleic Acids Res.* **2010**, *38*, W695–W699. [CrossRef]

58. Sievers, F.; Wilm, A.; Dineen, D.; Gibson, T.J.; Karplus, K.; Li, W.; Lopez, R.; McWilliam, H.; Remmert, M.; Söding, J.; et al. Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. *Mol. Syst. Biol.* **2011**, *7*, 539. [CrossRef]

59. Hall, T.A. BioEdit: A user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. *Nucleic Acids Symp. Ser.* **1999**, *41*, 95–98.

60. Tamura, K.; Stecher, G.; Peterson, D.; Filipski, A.; Kumar, S. MEGA6: Molecular Evolutionary Genetics Analysis Version 6.0. *Mol. Biol. Evol.* **2013**, *30*, 2725–2729. [CrossRef]

61. Rambaut, A.; Lam, T.T.; Max Carvalho, L.; Pybus, O.G. Exploring the temporal structure of heterochronous sequences using TempEst (formerly Path-O-Gen). *Virus Evol.* **2016**, *2*, vew007. [CrossRef]

62. Kosakovsky Pond, S.L.; Posada, D.; Gravenor, M.B.; Woelk, C.H.; Frost, S.D.W. GARD: A genetic algorithm for recombination detection. *Bioinformatics* **2006**, *22*, 3096–3098. [CrossRef] [PubMed]

63. Delport, W.; Poon, A.F.Y.; Frost, S.D.W.; Kosakovsky Pond, S.L. Datamonkey 2010: A suite of phylogenetic analysis tools for evolutionary biology. *Bioinformatics* **2010**, *26*, 2455–2457. [CrossRef] [PubMed]

64. Drummond, A.J.; Rambaut, A. BEAST: Bayesian evolutionary analysis by sampling trees. *BMC Evol. Biol.* **2007**, *7*, 214. [CrossRef] [PubMed]

65. Drummond, A.J.; Ho, S.Y.W.; Phillips, M.J.; Rambaut, A. Relaxed phylogenetics and dating with confidence. *PLoS Biol.* **2006**, *4*, 699–710. [CrossRef]

66. Ren, L.; Zhang, Y.; Li, J.; Xiao, Y.; Zhang, J.; Wang, Y.; Chen, L.; Paranhos-Baccalà, G.; Wang, J. Genetic drift of human coronavirus OC43 spike gene during adaptive evolution. *Sci. Rep.* **2015**, *5*, 11451. [CrossRef]

67. Kosakovsky Pond, S.L.; Frost, S.D.W. Not So Different After All: A Comparison of Methods for Detecting Amino Acid Sites Under Selection. *Mol. Biol. Evol.* **2005**, *22*, 1208–1222. [CrossRef]

68. Pond, S.L.K.; Frost, S.D.W. Datamonkey: Rapid detection of selective pressure on individual sites of codon alignments. *Bioinformatics* **2005**, *21*, 2531–2533. [CrossRef]

69. Cotten, M.; Watson, S.J.; Zumla, A.I.; Makhdoom, H.Q.; Palser, A.L.; Ong, S.H.; Al Rabeeah, A.A.; Alhakeem, R.F.; Assiri, A.; Al-Tawfiq, J.A.; et al. Spread, circulation, and evolution of the Middle East respiratory syndrome coronavirus. *MBio* **2014**, *5*, e01062-13. [CrossRef]

70. Storz, J.; Purdy, C.W.; Lin, X.; Burrell, M.; Truax, R.E.; Briggs, R.E.; Frank, G.H.; Loan, R.W. Isolation of respiratory bovine coronavirus, other cytocidal viruses, and Pasteurella spp from cattle involved in two natural outbreaks of shipping fever. *J. Am. Vet. Med. Assoc.* **2000**, *216*, 1599–1604. [CrossRef]

71. Cho, K.O.; Hoet, A.E.; Loerch, S.C.; Wittum, T.E.; Saif, L.J. Evaluation of concurrent shedding of bovine coronavirus via the respiratory tract and enteric route in feedlot cattle. *Am. J. Vet. Res.* **2001**, *62*, 1436–1441. [CrossRef]

72. Gulliksen, S.M.; Jor, E.; Lie, K.I.; Løken, T.; Åkerstedt, J.; Østerås, O. Respiratory infections in Norwegian dairy calves. *J. Dairy Sci.* **2009**, *92*, 5139–5146. [CrossRef]

73. Pardon, B.; De Bleecker, K.; Dewulf, J.; Callens, J.; Boyen, F.; Catry, B.; Deprez, P. Prevalence of respiratory pathogens in diseased, non-vaccinated, routinely medicated veal calves. *Vet. Rec.* **2011**, *169*, 278. [CrossRef] [PubMed]

74. O'Neill, R.; Mooney, J.; Connaghan, E.; Furphy, C.; Graham, D.A. Patterns of detection of respiratory viruses in nasal swabs from calves in Ireland: A retrospective study. *Vet. Rec.* **2014**, *175*, 351. [CrossRef] [PubMed]

75. Meyer, G.; Pelletier, C.; Herman, N.; Ducatez, M.; Cassard, H.; Salem, E. Vers une identification de nouveaux virus respiratoires bovins. *NEVA* **2015**, *8*, 80–86.

76. Hasoksuz, M.; Lathrop, S.L.; Gadfield, K.L.; Saif, L.J. Isolation of bovine respiratory coronaviruses from feedlot cattle and comparison of their biological and antigenic properties with bovine enteric coronaviruses. *Am. J. Vet. Res.* **1999**, *60*, 1227–1233. [PubMed]

77. Behura, S.K.; Tizioto, P.C.; Kim, J.; Grupioni, N.V.; Seabury, C.M.; Schnabel, R.D.; Gershwin, L.J.; Van Eenennaam, A.L.; Toaff-Rosenstein, R.; Neibergs, H.L.; et al. Tissue Tropism in Host Transcriptional Response to Members of the Bovine Respiratory Disease Complex. *Sci. Rep.* **2017**, *7*, 17938. [CrossRef]
78. Valentova, V.; Antonis, A.F.G.; Kovarcik, K. Restriction enzyme analysis of RT-PCR amplicons as a rapid method for detection of genetic diversity among bovine respiratory syncytial virus isolates. *Vet. Microbiol.* **2005**, *108*, 1–12. [CrossRef]

79. Sarmiento-Silva, R.E.; Nakamura-Lopez, Y.; Vaughan, G. Epidemiology, molecular epidemiology and evolution of bovine respiratory syncytial virus. *Viruses* **2012**, *4*, 3452–3467. [CrossRef]

80. Valarcher, J.F.; Schelcher, F.; Bourhy, H. Evolution of bovine respiratory syncytial virus. *J. Virol.* **2000**, *74*, 10714–10728. [CrossRef]

81. Vilcek, S.; Durkovic, B.; Kolesarova, M.; Paton, D.J. Genetic diversity of BVDV: Consequences for classification and molecular epidemiology. In *Proceedings of the Preventive Veterinary Medicine*; Elsevier: Amsterdam, The Netherlands, 2005; Volume 72, pp. 31–35.

82. FAO statistics. Available online: http://www.fao.org/faostat/en/#data/TM (accessed on 11 April 2019).

83. Collin, E.a.; Sheng, Z.; Lang, Y.; Ma, W.; Hause, B.M.; Li, F. Cocirculation of Two Distinct Genetic and Antigenic Lineages of Proposed Influenza D Virus in Cattle. *J. Virol.* **2015**, *89*, 1036–1042. [CrossRef]

84. Ducatez, M.F.; Pelletier, C.; Meyer, G. Influenza D Virus in Cattle, France, 2011–2014. *EID J. 2015*, *21*, 2011–2014. [CrossRef] [PubMed]

85. Flynn, O.; Gallagher, C.; Mooney, J.; Irvine, C.; Ducatez, M.; Hause, B.; McGrath, G.; Ryan, E. Influenza D virus in cattle, Ireland. *Emerg. Infect. Dis.* **2018**, *24*, 389–391. [CrossRef] [PubMed]

86. Foni, E.; Chiapponi, C.; Baioni, L.; Zanni, I.; Merenda, M.; Rosignoli, C.; Kyriakis, C.S.; Luini, M.V.; Mandola, M.L.; Nigrelli, A.D.; et al. Influenza D in Italy: Towards a better understanding of an emerging viral infection in swine. *Sci. Rep.* **2017**, *7*, 1–7. [CrossRef] [PubMed]

87. Callow, K.A.; Parry, H.F.; Sergeant, M.; Tyrrell, D.A. The time course of the immune response to experimental coronavirus infection of man. *Epidemiol. Infect.* **1990**, *105*, 435–446. [CrossRef] [PubMed]

88. Bidokhti, M.R.M.; Tråvén, M.; Krishna, N.K.; Munir, M.; Belák, S.; Alenius, S.; Cortey, M. Evolutionary dynamics of bovine coronaviruses: Natural selection pattern of the spike gene implies adaptive evolution of the strains. *J. Gen. Virol.* **2013**, *94*, 2036–2049. [CrossRef]

89. Vijgen, L.; Keyaerts, E.; Lemey, P.; Maes, P.; Reeth, K.V.; Nauwynck, H.; Pensaert, M.; Ranst, M.V. Evolutionary History of the Closely Related Group 2 Coronaviruses: Porcine Hemagglutinating Encephalomyelitis Virus, Bovine Coronavirus, and Human Coronavirus OC43. *J. Virol.* **2006**, *80*, 7270–7274. [CrossRef]

90. Sánchez, C.M.; Gebauer, F.; Suñé, C.; Mendez, A.; Dopazo, J.; Enjuanes, L. Genetic evolution and tropism of transmissible gastroenteritis coronaviruses. *Virology* **1992**, *190*, 92–105. [CrossRef]

91. Lau, S.K.P.; Lee, P.; Tsang, A.K.L.; Yip, C.C.Y.; Tse, H.; Lee, R.A.; So, L.-Y.; Lau, Y.-L.; Chan, K.-H.; Woo, P.C.Y.; et al. Molecular Epidemiology of Human Coronavirus OC43 Reveals Evolution of Different Genotypes over Time and Recent Emergence of a Novel Genotype due to Natural Recombination. *J. Virol.* **2011**, *85*, 11325–11337. [CrossRef]

92. Salem, E. Bronchopneumonies Infectieuses des Jeunes Bovins: De la Complexité du Microbiome aux Particularités évolutives et Cliniques de virus Respiratoires Encore Méconnus. Ph.D. Thesis, INP Toulouse, Toulouse, France, 24 October 2018.

<figure><img src="image_6.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).