key: cord-1012082-6706zxf7
authors: Castells, Matías; Victoria, Matías; Colina, Rodney; Musto, Héctor; Cristina, Juan
title: Genome-wide analysis of codon usage bias in Bovine Coronavirus
date: 2017-06-17
journal: Virol J
DOI: 10.1186/s12985-017-0780-y
sha: 62bd8c44c827603bc80a41b75fa8ba219969d48a
doc_id: 1012082
cord_uid: 6706zxf7

BACKGROUND: Bovine coronavirus (BCoV) belong to the genus Betacoronavirus of the family Coronaviridae. BCoV are widespread around the world and cause enteric or respiratory infections among cattle, leading to important economic losses to the beef and dairy industry worldwide. To study the relation of codon usage among viruses and their hosts is essential to understand host-pathogen interaction, evasion from host’s immune system and evolution. METHODS: We performed a comprehensive analysis of codon usage and composition of BCoV. RESULTS: The global codon usage among BCoV strains is similar. Significant differences of codon preferences in BCoV genes in relation to codon usage of Bos taurus host genes were found. Most of the highly frequent codons are U-ending. G + C compositional constraint and dinucleotide composition also plays a role in the overall pattern of BCoV codon usage. CONCLUSIONS: The results of these studies revealed that mutational bias is a leading force shaping codon usage in this virus. Additionally, relative dinucleotide frequencies, geographical distribution, and evolutionary processes also influenced the codon usage pattern. ELECTRONIC SUPPLEMENTARY MATERIAL: The online version of this article (doi:10.1186/s12985-017-0780-y) contains supplementary material, which is available to authorized users.

Coronaviruses belong to the family Coronaviridae and are the largest enveloped single-strand RNA viruses, ranging from 26 to 31 kilobases in genome size [1, 2] . These viruses infect a wide range of avian and mammalian species, and are responsible for enteric or respiratory infections [3] . There is a rising concern about the emergence of two human coronaviruses, Severe acute respiratory syndrome-related coronavirus (SARS-CoV) and Middle-East respiratory syndrome coronavirus (MERS-CoV), who emerged in 2002 and 2012, respectively [4, 5] . Both SARS-CoV and MERS-CoV have a zoonotic origin, revealing the importance of the control of coronaviruses associated with domestic animals in close contact with human populations [6] .

Coronaviruses consists of four genera named Alpha-, Beta-, Gamma-and Deltacoronavirus based on phylogenetic distance of highly conserved domains. In turn, Betacoronavirus genus is divided into four clades, namely A to D. Bovine coronavirus (BCoV) belongs to the Betacoronavirus genus clade A [7, 8] .

BCoV are widespread around the world and cause enteric or respiratory infections among cattle [9, 10] . These viruses are associated with different syndromes in cattle, ranging from neonatal calf diarrhea, winter dysentery in adult cattle, to respiratory infection in cattle of different age groups [11] .

BCoV infection leads to important economic losses to the beef and dairy industry throughout the world, associated with decreased performance, morbidity, mortality, direct cost of treatment of sick animals, and long-term effects on health and productivity of surviving calves [10, 12] .

BCoV is closely related to the Human coronavirus OC43 (HCoV-OC43), isolated in 1967 [13] . Recent studies revealed that HCoV-OC43 resulted from a zoonotic transmission from bovine to human [6] .

The redundancy of the genetic code provides evolution with the opportunity to adjust the efficiency and accuracy of protein production, while preserving the same amino acid sequence [13] . The relation of codon usage among viruses and their hosts may affect viral fitness, evasion from host's immune system and evolution [14] [15] [16] . Synonymous triplets are generally not used randomly, and the main forces that drive this bias from equal usage are natural selection and mutational biases [17, 18] . Therefore, the study of codon usage in viruses can reveal important information about virus evolution, regulation of gene expression and protein synthesis [19] .

In the present study, we performed comprehensive analyses of codon usage and composition of BCoV strains and explored the possible leading evolutionary determinants of the biases found.

Complete genome sequences for 15 BCoV strains were obtained from GenBank database (available at: http:// www.ncbi.nlm.nih.gov). For strain names and accession numbers see Additional file 1. By concatenation of different genome ORF's sequences, different datasets were constructed: one dataset includes the concatenation of nonstructural region ORFs (ORF1ab), a second one containing the structural region ORFs (hemagglutinin-esterase, spike, envelope, membrane, nucleocapsid) and a third including the concatenation of nonstructural and structural ORFs. For each strain the ORFs were aligned using the MUSCLE program [20] . The alignment of concatenation of nonstructural and structural ORFs is available in Additional file 2.

Codon usage, dinucleotide frequencies, base composition, the relative synonymous codon usage (RSCU) [21] , the effective number of codons (ENC) [22] , aromaticity (Aromo) and hydrophathicity (Gravy) values were calculated using the program CodonW (available at http://sourceforge.net/projects/codonw). The total G + C genomic content and G + C content at first, second and third codon positions were calculated using the Codon W program and EMBOSS Cusp program (available at http://emboss.bioinformatics.nl/ cgi-bin/emboss/cusp). The RSCU values of Bos taurus cells were obtained from Kazusa database (available at: http://www.kazusa.or.jp/codon/). Codon usage preferences in BCoV in relation to the codon usage of Bos taurus were established by means of the use of the codon adaptation index (CAI) [23] . CAI was calculated using the approach of Puigbo et al. [24] . This method permits to compare a given codon usage (BCoV) to a predefined reference set (Bos taurus). In order to show whether the BCoV genes are well adapted to the codon usage of the reference set, as measured by CAI, we constructed a dataset composed of 22 Bos taurus genes selected at random and obtained from ARSA at DNA Database of Japan (available at: http:/http://www.ddbj.nig.ac.jp/ arsa). The RCSU values of those 22 Bos taurus genes were found to be similar to the RCSU values obtained from the Kazusa database (see Additional file 3). Moreover, a strong positive correlation are found among the RSCU of those genes and the RSCU obtained from Kazusa database by Spearmans'rank correlation test [25] (r = 0.956, P < 0.00001).

Statistically significant difference among CAI values was determined by applying a Wilcoxom & Mann-Whitney test [25] . To discern if the statistically significant differences in the CAI values arise from codon preferences, we used e-CAI [26] to calculate the expected value of CAI (eCAI) at the 95% confident interval. A Kolmogorov-Smirnov test for the expected CAI was also performed [26] .

Correspondence analysis (COA) is a type of multivariate analysis that permits a geometrical representation of the sets of rows and columns in a dataset [27, 28] . Each ORF is represented as a 59-dimensional vector and each dimension corresponds to the RSCU value of each codon (excluding AUG, UGG and stop triplets). Major trends within a dataset can be established using measures of relative inertia and genes ordered according to their position along the different axes [29] . COA was performed on the RSCU values by means of the use of the CodonW program. Correlation analysis was performed using Spearman's rank correlation analysis method [25] .

In order to gain insight into the genetic variability and evolution of BCoV, a phylogenetic tree analysis was performed for all BCoV strains enrolled in these studies, using complete genome codes. Sequences were aligned using the MUSCLE program [20] . Once aligned, the FindModel program (available at: https:// www.hiv.lanl.gov/content/sequence/findmodel/findmodel.html) was used to identify the optimal evolutionary model that best fitted our sequence dataset.

Akaike information criteria (AIC) and the log of the likelihood (logL) indicated that the GTR+ Γ model was the most accurate (AIC = 60,864.73, logL = −30,423.36). Using this model, maximum likelihood trees were constructed using software from the MEGA 6 program [30] . As a measure of the robustness of each node, we employed the bootstrap method (500 replicas).

In order to gain insight into the degree of codon usage bias in BCoV, the ENC's values were calculated for the complete genome of all BCoV strains. A mean value of 43.78 ± 0.07 was obtained for BCoV strains included in these studies. Then, a plot of ENC versus GC3S (ENC plotted against G + C content at the third codon position) was constructed. An ENC-GC3S plot of genes whose codon choice is constrained only by a GC3 mutational bias, will lie on or just below the continuous curve of the predicted ENC values [31] . As shown in Fig. 1 , all points lie together under the expected ENC curve, indicating that G + C compositional constraints might play a role in ZIKV codon usage. Additionally, a correlation analysis between ENC and GC3S showed significant results (r = 0.811, p = 0.004).

The aromaticity (Aroma) and hydropathicity (Gravy) values of a given gene product can be indicative of the effect of translation or natural selection [32] . For these reasons, the effect of natural selection on BCoV codon usage was studied by correlation analysis between Gravy and Aroma values and nucleotide compositions at the third codon position and ENC values. No significant correlations between Gravy/Aroma values and nucleotide compositions at third codon position as well as between Gravy/Aroma and ENC values were observed (see Additional file 4).

These results indicate that codon bias in BCoV is related to mutational bias, whereas translational control, may be weak or absent in the reading frames of these viruses.

To study the trends in codon usage variation among different BCoV genomes, we performed a COA analysis. COA was performed on the RSCU values for the ORF1ab of each BCoV strain enrolled in these studies and we examined the distribution of the strains in the plane defined by the first two principal axes of COA. The first axis generated by the analysis accounts for 43.37% of the total variation, while the second axis accounts for 18.96%. The results of these studies are shown in Additional file 5: Figure S1A . Interestingly, different BCoV isolates are located at different positions in the plane, suggesting that trends in variation of BCoV can be observed. Moreover, these trends correlate with the genetic variability observed by phylogenetic analysis (see Additional file 5: Figure S1B ).

To compare the codon usage preferences of BCoV with those of Bos taurus, the RSCU values of the codons in nonstructural, structural and complete BCoV genome ORFs were calculated and compared with those of B. taurus. The results of these studies are shown in Table 1 .

The frequencies of codon usage in BCoV ORFs are significantly different in relation to B. taurus ones. Indeed, highly biased frequencies were found for UUU (Phe), UUA (Leu), UUG (Leu), CUU (Leu), AUU (Ile), GUU (Val), UAU (Tyr), CAU (His), CAA (Gln), AAU (Asn), GAU (Asp), UCU (Ser), CCU (Pro), ACU (Thr), GCU (Ala), UGU (Cys), CGU (Arg), AGU (Ser) and GGU (Gly). As can be seen, most of the highly preferred codons are U-ending and UpU containing codons, which strongly suggests that mutational bias is a main force shaping codon usage in BCoV (see Table 1 ). Moreover, most of the highly decreased codons frequencies with respect to B. taurus cells are C-ending codons, also suggesting a strong mutational bias in the use of these codons (Table 1) . A Wilcoxom & Mann-Whitney test on the frequencies of BCoV codon usage among nonstructural and structural genome regions revealed no significant differences among both regions (T = 1983, p-value = 0.758).

In this study, a CAI metric was used as a measure of relative adaptedness of BCoV codon usage to Bos taurus host. CAI values for all triplets were calculated for the complete genome ORFs of BCoV strains enrolled in these studies, using B. taurus codon usage as the reference set. The results of these studies are shown in Table 2 . A mean value of 0.638 was obtained for BCoV genes in relation to B. taurus; while a mean CAI value of 0.756 was obtained for a Bos taurus sequence dataset in relation to the same reference set (see Table 2 ). In order to evaluate if the differences were statistically significant, we performed a Wilcoxon & Mann-Whitney test. The results of this test revealed that the differences in CAI values are statistically significant (T = 0, p-value <0.001). The CAI value obtained for Bos taurus genes is higher than the one obtained for BCoV in relation to Bos taurus codon usage (see Table 2 ). This reveals that BCoV genes are relatively less adapted to Bos taurus than Bos taurus genes itselves.

In order to discern if the statistically significant differences in CAI values arise from codon preferences [24] , the expected CAI (e-CAI) values were calculated for BCoV complete genome ORFs sequences in relation to B. taurus codon usage reference set. The e-CAI algorithm [26] generated 500 random sequences with the same nucleotide and amino acid composition as the sequences of interest (in this case BCoV sequences). Then, we calculated the CAI values for all of them, and a Kolmogorov-Smirnov test for the e-CAI of these random sequences was performed in order to show if the generated sequences follow a normal distribution. The results of these studies revealed an e-CAI value of 0.656. Kolmogorov-Smirnov test revealed a normal distribution of the generated sequences (Kolmogorov-Smirnov test of e-CAI value of 0.028, which is below the critical value of 0.061). To avoid the effect of extreme compositional constraint and to make sure that CAI is directly correlated with codon usage preferences, Puigbo et al. [26] suggested that if the eCAI value of a gene is higher than its CAI value, it may be considered as evidence of codon usage adaptation. For all BCoV strains, CAI values were found to be lower than their corresponding eCAI values, when compared against Bos taurus (eCAI = 0.656, p < 0.05).

Taking all these results together, our studies revealed that the CAI values for BCoV genes are different from the CAI values obtained for B. taurus sequences and these differences are related to codon usage preferences.

It has been previously shown that dinucleotide biases can play a role in codon usage bias [29] . In order to determine if this is the case in BCoV, the relative abundances of the 16 dinucleotides in BCoV complete genome ORFs were established. The results of these studies are shown in Table 3 . (3) , percentage of G + C genomic content at codon positions 1 through 3, respectively. In all cases, mean ± standard deviation values are shown

As can be seen, the relative abundance of UpU, UpA and ApU showed a strong deviation from the expected frequencies (i.e. 1.0) (a mean of 2.02, 1.44 and 1.47, respectively), while UpC, CpC, ApC and GpC frequencies were markedly underrepresented (a mean of 0.60, 0.44, 0.70 and 0.70, respectively). The relative abundance of CpG and GpC also showed a strong deviation from the expected frequencies (i.e. 1.0) (a mean of 0.25 and 0.70, respectively).

In these studies, we first study the general codon usage pattern in BCoV. When the ENC's values were calculated for the complete genome of all BCoV strains, a mean value of 43.78 ± 0.07 was obtained. Since the ENC values obtained are >40, these results suggest a relatively conserved codon usage bias among different BCoV genomes. This is in agreement with previous reports in other members of the family, like SARS-CoV (mean ENC = 48.99) [33] , the avian coronavirus Infectious bronchitis virus (ENC = 42.79) [34] or the Porcine epidemic diarrhea virus (ENC = 47.91) [35] . The ENC-GC3S plot revealed that all values obtained for BCoV lie bellow the continuous curve of predicted ENC values in absence of selection, revealing that G + C compositional constrain play a role in BCoV codon usage (Fig. 1) .

Moreover, no significant correlations between both Gravy and Aroma values with nucleotide composition at the third codon position, as well as ENC values (see Additional file 3). This indicates that the role of translational selection in BCoV codon usage bias is weak or absent. Taking all together, the results of these studies suggests that mutational bias is a main force shaping codon usage in this virus.

Then, in order to study the trends in codon usage variation across BCoV strains, a COA analysis was performed on the RSCU values for the ORF1ab of each BCoV strain enrolled in these studies. Different BCoV isolates are located at different positions in the plane defined by the first two principal axes of COA, revealing that trends in codon usage variation can be observed (see Additional file 5: Figure S1A ). This is also in agreement with the results found in the phylogenetic analysis of BCoV strains enrolled in these studies (see Additional file 5: Figure S1B ) and with recent studies on genetic variability of coronaviruses, showing that BCoVs strains are distributed on three main sub-clusters named C1, C2, and C3 [6] . Subcluster C1 includes BCoVs from America and Asia, sub-cluster C2 includes BCoVs from Europe and subcluster C3 includes prototype, vaccine, or attenuated BCoV strains [6] . As can be seen in Additional file 5: Figure S1A , BCoV strains Mebus and Quebed, who belong to C3 cluster, are situated in different positions in the plane defined by the first two axes of COA in relation to C1 strains. Moreover, different genetic lineages of C1 sub-cluster are also located at different positions in the plane (Additional file 5: Figure S1A ). This is in agreement with the results found in the phylogenetic analysis of BCoVs enrolled in these studies (Additional file 5: Figure S1B ). Moreover, BCoV strain BCV-AKS-01, who shows a more distant genetic relation with C3 and C1 BCoVs enrolled in the studies (Additional file 5: Figure S1B ), is situated in a different position in the plane defined by the first two axes of COA (Additional file 5: Figure S1A ). The results of these studies revealed that evolutionary processes also influenced the codon usage pattern of BCoV.

In these studies, significant differences in codon usage frequencies were found among BCoV and B. taurus ones (Table 1) . Interestingly, recent studies on human coronaviruses revealed significant biases in nucleotide composition [36] . While the A/G bias is a relatively stable property among coronaviruses, the C/U bias differs significantly in each virus type, with U-counts ranging from 30.7% (SARS-CoV) to 40.3% (HCoV-HKU) and C-counts from 20.3% (MERS-CoV) to 12.9% (HCoV-HKU) [36] . The U-count for BCoV revealed a value of 35.6%, while the C-count shows a value of 15.1%. This is in agreement with the results found in this work and may help to explain the bias Mean values of BCoV strains relative dinucleotide ratios ± standard deviation. b Correlation analysis between the first axis in COA and the sixteen dinucleotides frequencies in BCoV genes is shown found in BCoV genome composition and codon usage. Moreover, the CAI values for BCoV genes resulted to be significantly different from the ones obtained for B. taurus, revealing differences in codon usage preferences ( Table 2 ).

The results of these studies revealed that the relative abundance of CpG and GpC dinucleotides showed a strong deviation from the expected frequencies (Table 3) . The under-representation of CpG might be due to its immunostimulatory properties as recognition of unmethylated CpG by Toll like receptor 9 (TLR9), which leads to activation of several immune response pathways in the host [37] . Moreover, an increase in CpG dinucleotide frequency has been shown to lead to attenuation of replication in RNA viruses [38] . Cytosine deamination and selection against CpG motifs have been proposed as two independent selection forces that shape codon usage bias in coronaviruses [39] , suggesting that immune selection may play a role in the observed BCoV codon usage bias. This is in agreement with the results of this work and indicates that the composition of dinucleotides also determines the variation in synonymous codon usage among BCoV.

Mesoniviridae: a proposed new family in the order Nidovirales formed by a single species of mosquito-borne viruses

Nidovirales: evolving the largest RNA virus genome

Coronavirus diversity, phylogeny and interspecies jumping

Severe acute respiratory syndrome

Isolation of a novel coronavirus from a man with pneumonia in Saudi Arabia

Comparative molecular epidemiology of two closely related coronaviruses, bovine coronavirus (BCoV) and human coronavirus OC43 (HCoV-OC43), reveals a different evolutionary pattern

Ratification vote on taxonomic proposals to the International Committee on Taxonomy of Viruses

Discovery of seven novel Mammalian and avian coronaviruses in the genus deltacoronavirus supports bat coronaviruses as the gene source of alphacoronavirus and betacoronavirus and avian coronaviruses as the gene source of gammacoronavirus and deltacoronavirus

Detection of respiratory and enteric shedding of bovine coronaviruses in cattle in an Ohio feedlot

Bovine coronavirus associated syndromes

Bovine respiratory coronavirus

Development and performance evaluation of calf diarrhea pathogen nucleic acid purification and detection workflow

Synonymous codon usage in Escherichia coli: selection for translational accuracy

Modulation of poliovirus replicative fitness in HeLa cells by deoptimization of synonymous codon usage in the capsid region

Reduction of the rate of poliovirus protein synthesis through large-scale codon deoptimization causes attenuation of viral virulence by lowering specific infectivity

Hepatitis A virus adaptation to cellular shutoff is driven by dynamic adjustments of codon usage and results in the selection of populations with altered capsids

Forces that influence the evolution of codon bias

What we know and what we should know about codon usage

Genome-wide analysis of codon usage and influencing factors in Chikungunya viruses

MUSCLE: a multiple sequence alignment method with reduced time and space complexity

An evolutionary perspective on synonymous codon usage in unicellular organisms

Accounting for background nucleotide composition when measuring codon usage bias

The codon adaptation index -a measure of directional synonymous codon usage bias, and its potential applications

CAIcal: a combined set of tools to assess codon usage adaptation

Free Statistics Software, Office for Research Development and Education, version 1.1.23-r7

a novel server to estimate an expected value of Codon Adaptation Index (eCAI)

Codon usage bias and the evolution of influenza A viruses. Codon Usage Biases of Influenza Virus

Theory and Applications of Correspondence Analysis

Analysis of synonymous codon usage in classical swine fever virus

MEGA6: Molecular Evolutionary Genetics Analysis version 6.0

Analysis of codon usage bias and base compositional constraints in Iridovirus genomes

expressivity and aromaticity are the major trends of amino-acid usage in 999 Escherichia coli chromosomeencoded genes

Analysis of synonymous codon usage in 11 human bocavirus isolates

Avian coronavirus spike glycoprotein ectodomain shows a low codon adaptation to Gallus gallus with virus-exclusive codons in strategic amino acids positions

Characterization of the porcine epidemic diarrhea virus codon usage bias

On the biased nucleotide composition of the human coronavirus RNA genome

Evolutionary basis of codon usage and nucleotide composition bias in vertebrate DNA viruses

RNA virus attenuation by codon pair deoptimisation is an artefact of increases in CpG/ UpA dinucleotide frequencies

Cytosine deamination and selection of CpG suppressed clones are the two major independent biological forces that shape codon usage bias in coronaviruses

Acknowledgements MC acknowledges support from ANII and INIA through PhD fellowships.

This work was found by Instituto Nacional de Investigación Agropecuaria (INIA) through project N-15156, Agencia Nacional de Investigación e Innovación (ANII) through project PE_ALI_2009_1_1603, Fondo MaríaViñas and PEDECIBA, Uruguay. We acknowledge Comisión Sectorial de Investigación Científica (CSIC), Universidad 18 de la República, Uruguay, for encoragement support through Grupos I + D program.

The datasets generated and analyzed during the current study are available in Additional file 2.

The results of these studies revealed significant differences in codon preferences in BCoV genes in relation to B. taurus codon usage. The overall codon usage among BCoV strains is similar. All U-ending codons are highly frequent codons, which strongly suggests that mutational bias is a leading force shaping codon usage in this virus. G + C compositional constraint influences the codon usage of BCoV. Dinucleotide composition also plays a role in the overall pattern of BCoV codon usage. Authors' contributions JC conceived of the study, and participated in its design and coordination. MC have made substantial contributions to the acquisition of data and analysis. MV, RC and HM have been involved in the analysis and in interpretation of the results found in this work, as well as in revising the manuscript critically for important intellectual content. JC wrote the paper. All authors read and approved the final manuscript.

The authors declare that they have no competing interests.

Ethics approval and consent to participate Not applicable.• We accept pre-submission inquiries • Our selector tool helps you to find the most relevant journal Submit your next manuscript to BioMed Central and we will help you at every step: