Bovine Coronavirus: Variability, Evolution, and Dispersal Patterns of a No Longer Neglected Betacoronavirus

Abstract

Bovine coronavirus (BoCV) is an important pathogen of cattle, causing severe enteric disease and playing a role in the bovine respiratory disease complex. Similar to other coronaviruses, a remarkable variability characterizes both its genome and biology. Despite their potential relevance, different aspects of the evolution of BoCV remain elusive. The present study reconstructs the history and evolution of BoCV using a phylodynamic approach based on complete genome and spike protein sequences. The results demonstrate high mutation and recombination rates affecting different parts of the viral genome. In the spike gene, this variability undergoes significant selective pressures-particularly episodic pressure-located mainly on the protein surface, suggesting an immune-induced selective pressure. The occurrence of compensatory mutations was also identified. On the contrary, no strong evidence in favor of host and/or tissue tropism affecting viral evolution has been proven. The well-known plasticity is thus ascribable to the innate broad viral tropism rather than mid- or long-term adaptation. The evaluation of the geographic spreading pattern clearly evidenced two clusters: a European cluster and an American-Asian cluster. While a relatively dense and quick migration network was identified in the former, the latter was dominated by the primary role of the United States (US) as a viral exportation source. Since the viral spreading pattern strongly mirrored the cattle trade, the need for more intense monitoring and preventive measures cannot be underestimated as well as the need to enforce the vaccination of young animals before international trade, to reduce not only the clinical impact but also the transferal and mixing of BoCV strains.

Keywords: bovine coronavirus; evolution; host; phylodynamics; phylogeography; selection.

Figure 1

Recombination distribution plot. The significance (p-value) of hot and cold-recombination spots is reported along the genome.

Figure 2

Upper figure: Boxplot (a) and density plot (b) of the evolutionary rate posterior probability. Lower figure: Boxplot (c) and density plot (d) of time to most recent common ancestor (tMRCA) posterior probability. Results are reported for all considered genes. The 95% high posterior density (95HPD) intervals are reported for both figures. The median values are displayed as a bold line in the boxplots.

Figure 3

Relative genetic diversity (Ne x t) of the bovine coronavirus (BoCV) population over time. The mean value is reported in black, while the upper and lower 95HPD values are reported as a shaded area. The mean and upper 95HPD values of the time to Most Recent Common Ancestor (tMRCA) are reported as black dotted lines.

Figure 4

Time-scaled phylogenetic tree based on BoCV spike protein. The estimated ancestral geographic location has been color-coded (different nuances of green, red, and blue have been selected for American, Asian, and European countries), while the posterior probability is depicted as a circle on the corresponding node, whose size is proportional to the estimated value. The statistically supported migration routes have been represented in the upper left insert.

Figure 5

Well-supported migration paths between countries are depicted. The arrows indicate the directionality of the process, while the edge color is proportional to the base-10 logarithm of the migration rate. The location of each country has been matched with its centroid.

Figure 6

Time-scaled phylogenetic tree based on BoCV spike protein. The estimated ancestral host has been color-coded, while its posterior probability is depicted as a circle on the corresponding node, whose size is proportional to the estimated value. The statistically supported migration routes have been represented in the upper left insert.

Figure 7

Upper (left) and lateral (right) view of the quaternary structure of the BoCV spike protein reconstructed using homology modeling. The appearance has been edited to highlight different selective pressure features. In the color-coded monomer, the difference between non-synonymous and synonymous substitution rates (dN–dS), calculated using FUBAR, is reported on the protein surface ranging from purple (dN > dS) to light blue (dN < dS). In the gray-colored monomer, sites under episodic diversifying selection are reported in red. In the remaining monomer, colored in orange, the receptor-binding domain (RBD) is highlighted in green. In the transparency, the ribbon structure is visualized. In the central insert, the RBD is magnified (upper image), and the strength of selective pressures (lower image) acting in that region is represented using the previously described scale. A more detailed representation of the overall protein structure is reported in the Animations (S1 and S2).