key: cord-0001923-7pbfxcmp
authors: Nasrullah, Izza; Butt, Azeem M; Tahir, Shifa; Idrees, Muhammad; Tong, Yigang
title: Genomic analysis of codon usage shows influence of mutation pressure, natural selection, and host features on Marburg virus evolution
date: 2015-08-26
journal: BMC Evol Biol
DOI: 10.1186/s12862-015-0456-4
sha: f2e06c174008ed721d0a244309601c9923548901
doc_id: 1923
cord_uid: 7pbfxcmp

BACKGROUND: The Marburg virus (MARV) has a negative-sense single-stranded RNA genome, belongs to the family Filoviridae, and is responsible for several outbreaks of highly fatal hemorrhagic fever. Codon usage patterns of viruses reflect a series of evolutionary changes that enable viruses to shape their survival rates and fitness toward the external environment and, most importantly, their hosts. To understand the evolution of MARV at the codon level, we report a comprehensive analysis of synonymous codon usage patterns in MARV genomes. Multiple codon analysis approaches and statistical methods were performed to determine overall codon usage patterns, biases in codon usage, and influence of various factors, including mutation pressure, natural selection, and its two hosts, Homo sapiens and Rousettus aegyptiacus. RESULTS: Nucleotide composition and relative synonymous codon usage (RSCU) analysis revealed that MARV shows mutation bias and prefers U- and A-ended codons to code amino acids. Effective number of codons analysis indicated that overall codon usage among MARV genomes is slightly biased. The Parity Rule 2 plot analysis showed that GC and AU nucleotides were not used proportionally which accounts for the presence of natural selection. Codon usage patterns of MARV were also found to be influenced by its hosts. This indicates that MARV have evolved codon usage patterns that are specific to both of its hosts. Moreover, selection pressure from R. aegyptiacus on the MARV RSCU patterns was found to be dominant compared with that from H. sapiens. Overall, mutation pressure was found to be the most important and dominant force that shapes codon usage patterns in MARV. CONCLUSIONS: To our knowledge, this is the first detailed codon usage analysis of MARV and extends our understanding of the mechanisms that contribute to codon usage and evolution of MARV. ELECTRONIC SUPPLEMENTARY MATERIAL: The online version of this article (doi:10.1186/s12862-015-0456-4) contains supplementary material, which is available to authorized users.

The Marburg virus (MARV) is a negative-sense singlestranded RNA virus with a genome size of 19 kb that encodes seven genes in a linear order. MARV belongs to family Filoviridae, which also includes the highly pathogenic Ebola virus (EBOV). The first documented evidence of MARV was in 1967 in laboratory workers and scientists at facilities in Germany and the former Yugoslavia via infected monkeys that were imported from north-western Uganda [1] . MARV is a zoonotic virus and has been detected in both infected and healthy Egyptian fruit bats (Rousettus aegyptiacus) in endemic areas in Africa; therefore, R. aegyptiacus are considered as its natural host. This is most likely the reason that MARV outbreaks have been mostly associated with individuals such as mine workers or tourists in the regions that these bats inhabit [2] . The typical symptoms include general malaise, acute fever, abdominal cramping, bleeding disorders, and shock [3] . Similar to its highly pathogenic cousin EBOV that is the cause of a recent ongoing outbreak, MARV also causes fatal viral hemorrhagic fever in humans and non-human primates with a fatality rate of up to 90 %. Therefore, there is a need for a detailed understanding of replication and evolution of this virus [3] [4] [5] .

It is known that the genetic code shows redundancy and most of the amino acids can be translated by more than one codon. This redundancy represents a key step in modulating the efficiency and accuracy of protein production while maintaining the same amino acid sequence of the protein. Alternative codons within the same group that codes for the same amino acid are often termed 'synonymous' codons, although their corresponding tRNAs might differ in relative abundance in cells and in the speed by which they are recognized by the ribosome. However, the synonymous codons are not randomly chosen within and between genomes, which is referred to as codon usage bias [6, 7] . This phenomenon has been observed in a wide range of organisms, from prokaryotes to eukaryotes and viruses. Studies on codon usage have identified several factors that could influence codon usage patterns, including mutation pressure, natural or translational selection, secondary protein structure, replication and selective transcription, hydrophobicity and hydrophilicity of the protein, and the external environment [8] [9] [10] [11] [12] [13] . Moreover, considering the virus's genome size and other viral features, such as dependence on host's machinery for key processes, including replication, protein synthesis, and transmission, compared with prokaryotic and eukaryotic genomes, the interplay of codon usage among viruses and their hosts is expected to affect overall viral survival, fitness, evasion from host's immune system, and evolution [11, 14] . Therefore, knowledge of codon usage in viruses can reveal information about molecular evolution as well as improve our understanding of the regulation of viral gene expression and aid in vaccine design, where efficient expression of viral proteins may be required to generate immunity.

In the present study, we report genome-wide comprehensive analyses of codon usage and various factors that have contributed to the molecular evolution in MARV.

The complete genome sequences of 63 MARV strains were obtained from the National Center for Biotechnology (NCBI) GenBank database (http://www.ncbi.nlm.nih.gov). The accession numbers and demographics of the selected MARV genomes are provided in Additional file 1: Table S1 .

Identification of potential recombinant events in the MARV genomes were determined with the Recombination Detection Program (RDP) 4 Beta (version 4.27) software suite [15] , which incorporates several phylogeneticsubstitution and distance-based methods, including GENECONV [16] , RDP [17] , MaxChi [18] , Chimaera [19] , Bootscan [20] , SiScan [21] , 3Seq [22] , and LARD [23] . The P-value cut-off was set to 0.05 in all analyses, and the Bonferroni correction was applied. The default settings were used for all analyses.

The following compositional properties were calculated for the coding sequences of MARV genomes: (i) overall frequency of occurrence of the nucleotides (A%, C%, T/U%, and G%); (ii) frequency of each nucleotide at the third site of the synonymous codons (A 3S %, C 3S %, U 3S %, and G 3S %); (iii) frequencies of occurrence of nucleotides G + C at the first (GC 1S ), second (GC 2S ), and third synonymous codon positions (GC 3S ); (iv) mean frequencies of nucleotides G + C at the first and second position (GC 1,2S ); and (v) overall GC and AU content. The codons AUG and UGG are the only codons for Met and Trp, respectively, and the termination codons UAA, UAG, and UGA do not encode any amino acids. Therefore, these five codons are not expected to exhibit any usage bias and were therefore excluded from the analysis.

The RSCU values for all of the coding sequences of MARV genomes were calculated to determine the characteristics of synonymous codon usage without the confounding influence of amino acid composition and coding sequence size of different gene samples following a previously described method [24] . The RSCU index was calculated as follows:

where g ij is the observed number of the ith codon for the jth amino acid, which has n i kinds of synonymous codons. RSCU values represent the ratio between the observed usage frequency of one codon in a gene sample and the expected usage frequency in the synonymous codon family, given that all codons for the particular amino acid are used equally. The synonymous codons with RSCU values > 1.0 have positive codon usage bias and were defined as abundant codons, whereas those with RSCU values < 1.0 have negative codon usage bias and were defined as less-abundant codons. When the RSCU value is 1.0, it means there is no codon usage bias for that amino acid and the codons are chosen equally or randomly [25] . Moreover, the synonymous codons with RSCU values > 1.6 and < 0.6 were treated as over-represented and under-represented codons, respectively [26] .

CAI analysis is a quantitative method that predicts the expression level of a gene based on its coding sequence. CAI values range from 0 to 1. The most frequent codons have the highest relative adaptiveness towards its host, and sequences with higher CAIs are suggested to be preferred over those with lower CAIs [27] . The CAI analysis for MARV genes was performed using CAIcal server [28] . The synonymous codon usage patterns of the viral hosts (H. sapiens and R. aegyptiacus) were used as references. Non-synonymous codons and termination codons were excluded from the calculation. The reference datasets for H. sapiens and R. aegyptiacus were obtained from the Codon Usage Database [29] . The correlation analysis between CAI and ENC values was performed to determine the relative influence of mutation and selection. If selection is preferred over mutation, the correlation (r) between the two quantities should be very high (r → −1). In contrast, if mutation force is more important, r should approach 0 (no correlation) [30, 31] .

ENC analysis was used to quantify the absolute codon usage bias by evaluating the degree of codon usage bias exhibited by the MARV coding sequences, regardless of gene length and the number of amino acids. ENC values range from 20, which indicates extreme codon usage bias using only one of the possible synonymous codons for the corresponding amino acid, to 61, which indicates no bias using all possible synonymous codons equally for the corresponding amino acid. The larger the extent of codon preference in a gene, the smaller the ENC value. It is also generally accepted that genes have a significant codon bias when the ENC value is less than or equal to 35 [32, 33] . ENC was calculated using the following formula:

Where F k (k = 2, 3, 4, 6) is the mean of F k values for the k-fold degenerate amino acids, which is estimated using the following formula:

where n is the total number of occurrences of the codons for that amino acid and

where n i is the total number of occurrences of the ith codon for that amino acid. Genes for which the codon choice is only constrained by a mutation bias will lie on or just below the curve of the expected ENC values. Therefore, to elucidate the relationship between GC 3S and ENC values, the expected ENC values for different GC 3S were calculated as follows:

where s represents the given GC 3S % [32] .

PCA is a multivariate statistical method that is used to explore the relationships between variables and samples.

In the present study, PCA was used to analyze the major trends in codon usage patterns among MARVs coding sequences. PCA involves a mathematical procedure that transforms correlated variables (RSCU values) into a smaller number of uncorrelated variables called principal components. To minimize the effect of amino acid composition on codon usage, each coding sequence was represented as a 59 dimensional vector, and each dimension corresponded to the RSCU value of each sense codon, which only included synonymous codons for a particular amino acid excluding the codons AUG, UGG, and the three stop codons.

The neutrality plot or neutral evolution analysis was performed to determine and compare the extent of influence of mutation pressure and natural selection on the codon usage patterns of MARV by plotting the P 12 (GC 1,2S ) values of the synonymous codons against the P 3 (GC 3S ) values. In the plot, the regression coefficient against P 3 is regarded as the mutation-selection equilibrium coefficient and the evolutionary speed of the mutation pressure and natural selection pressure is expressed as the slope of a regression line. If all of the points lie along the diagonal distribution, no significant difference exists at the three codon positions, and there is no or weak external selection pressure. Alternatively, if the regression curve tends to be sloped or parallel to the horizontal axis, it means that the variation correlation between GC 1,2S and GC 3S is very low. Therefore, the regression curve effectively measures the degree of neutrality when selecting the effect that dominates evolution.

The Parity rule 2 (PR2) plot analysis was performed to investigate the impact of mutation and selection pressure on codon usage of genes. PR2 is a plot of AU-bias [A 3 /(A 3 + U 3 )] as the ordinate and GC-bias [G 3 /(G 3 + C 3 )] as the abscissa at the third codon position of the four-codon amino acids of entire genes. In this plot, the center of the plot, where both coordinates are 0.5, is the place where A = U and G = C (PR2), with no biasness between influence of mutation and selection rates (substitution rates) [34, 35] .

The RSCU values of MARVs and its two hosts, H. sapiens and R. aegyptiacus, were compared to determine influence of the host. The codon usage data of H. sapiens and R. aegyptiacus were obtained from the Codon Usage Database [29] . Furthermore, the influence of the overall codon usage patterns of hosts on the formation of the overall codon usage of viruses, defined as the similarity index D(A,B) [36] , was calculated as follows:

where R(A,B) is defined as a cosine value of an included angle between A and B spatial vectors and represents the degree of similarity between MARV and host overall codon usage pattern. a i is defined as the RSCU value for a specific codon among 59 synonymous codons of MARV coding sequence. b i is termed as the RSCU value for the same codon of the host. D(A,B) represents the potential effect of the overall codon usage of the host on that of MARV, and its value ranges from 0 to 1.0 [36] .

The relative abundance of dinucleotides in the coding regions of MARV genomes was calculated using a previously described method [37] . A comparison of actual and expected dinucleotide frequencies of the 16 dinucleotides in coding regions of MARV was also undertaken. The odds ratio was calculated using the following formula:

where f x denotes the frequency of nucleotide X, f y denotes the frequency of nucleotide Y, f y f x denotes the expected frequency of the dinucleotide XY, and f xy the frequency of the dinucleotide XY. This was calculated for each dinucleotide. As a conservative criterion, for P xy > 1.23 or < 0.78, the XY pair is considered to be over-represented or under-represented, respectively, in terms of relative abundance compared with a random association of mononucleotides.

Correlation analysis was carried out to identify the relationships between nucleotide composition, PCA, and codon usage patterns of MARV using Spearman's rank correlation analysis. All statistical analyses were carried out using SPSS 17 (SPSS Inc., Chicago, IL, USA) for Windows.

It has been previously shown that occurrence of recombination events at either gene or genome levels can influence codon usage bias patterns. For example, recombination can influence the effect of natural selection on codon usage [38] [39] [40] [41] . Therefore, to avoid influence of recombination on codon analysis, we first performed recombination analysis on the 63 MARV genomes. No evidence of recombination was found among MARV genomes. Therefore, coding sequences of all 63 of the initially selected MARV genomes were included in codon usage analysis as discussed in following sections.

To determine the potential impact of nucleotide constraints on codon usage, nucleotide composition analysis was performed. It was found that the A and U nucleotides were most abundant in MARV coding sequences with a mean composition of 31.9 and 27.7 %, respectively, compared with C (20.8 %) and G (19.6 %). The nucleotide composition at the third position of synonymous codons (A 3S , U 3S , G 3S , C 3S ) showed that the mean A 3S (31.3 %) and U 3S (33.0 %) were also highest compared with G 3S (17.7 %) and C 3S (18 %) ( Table 1 ). The mean AU and GC compositions were determined to be 59.6 and 40.4 %, respectively, highlighting that there is an AU-rich composition of MARV coding sequences. The analysis of nucleotide composition at first, second, and third positions of synonymous codons showed that GC 1S values ranged from 46.3 to 46.7 %, with a mean of 46.5 % and standard deviation (SD) of 0.12. GC 2S values ranged from 38.8 to 39.1%, with a mean of 39.0 % and an SD 0.08. GC 1,2S values ranged from 42.6 to 42.8 %, with an average of 42.7 % and SD of 0.05. In the case of GC 3S , the values ranged from 35.1 to 36.3 %, with a mean of 35.67 % and SD of 0.36; alternatively, the AU 3S values ranged from 63.6 to 65.0 %, with a mean of 64.33 % and an SD of 0.36. These data confirmed that a substantial portion of MARV coding sequences are composed of A and U nucleotides (Table 1) .

The patterns of synonymous codon usage in MARV coding sequences were assessed by RSCU analysis. All of the 18 most abundantly used codons for their corresponding amino acids in MARV coding sequences were A/U-ended and exhibited an equal distribution of A and U (A-ended: 9; U-ended: 9) ( Table 2 , Additional file 2: Figure S1 ). Analysis of over-and under-representation of codons showed that four out of 18 preferred codons had RSCU values >1.6. These are UUA(L), ACA(T), UCA(S), and AGA(R), whereas the RSCU values of the remaining preferred codons were also found to be >0.6 and <1.6. However, the under-represented (RSCU <0.6) and nonpreferred codons were all G/C-ended (Table 2) . Nucleotide composition and RSCU analyses showed that selection of the preferred codons has been mostly influenced by compositional constraints (A and U in this case), which accounts for the presence of mutation pressure.

To estimate the degree of codon usage bias within coding sequences of different isolates of MARV, ENC were calculated. The ENC values among MARV coding sequences ranged from 53.4 to 54.9, with a mean of 54.2 (ENC > 40) and an SD of 0.35 (Table 1) , indicating a relatively stable and conserved genomic composition among different MARV coding sequences.

To determine the trends in codon usage variation among coding sequences of different MARV isolates, we performed PCA on the RSCU values, which were examined as a single dataset (Fig. 1a) . The first principal axis (f' 1 ) accounted for 65.55 % of the total variation, and the next three axes (f' 2 -f' 4 ) accounted for 14.17, 10.48, and 2.36 % of the total variation in synonymous codon usage, respectively. Next, we plotted principal axes based on geographical locations of MARV isolates (Fig. 1b) . Three separate clusters were observed. Cluster A, which formed the largest cluster, consisted of isolates from the Demographic Republic of Congo (DRC), Uganda, and a single isolate from South Africa. The DRC isolates formed the majority of Cluster A. Cluster B was dominated by isolates from Uganda, Kenya, and a single isolate from DRC, whereas Cluster C consisted of all of the isolates from Angola and single isolates from Germany, the Netherlands, Uganda, and DRC.

To determine whether the biased codon choices are restricted to highly biased genes, the relation between A and U content and G and C content in four-fold degenerate codon families (alanine, arginine, glycine, leucine, proline, serine, threonine and valine) were analyzed by PR2 plot (Fig. 2) . It was found that A and U were used more frequently than G and C in MARV four fold degenerate codon families. This shows that preference towards codon To determine whether the patterns of codon usage have been influenced by mutation pressure, ENC-plot, correlation, linear regression and neutrality plot analyses were performed. In case of ENC-plot, The GC 3S values were plotted against ENC, which showed that all spots clustered slightly below on the left side of the expected curve (Fig. 3) . This indicates that mutational pressure has dominated in shaping codon usage patterns of MARVs.

In the next step, correlation analysis among the nucleotide compositions, codon compositions, and ENC values was performed (Table 3) . Several strong and significant correlations (P < 0.01; P < 0.05) were observed between nucleotide compositions and codon compositions. GC and GC 1,2S were also compared with GC 3S and highly significant positive correlations (GC 1,2S versus GC 3S : r = 0.747, P < 0.01; GC versus GC 3S : r = 0.828, P < 0.01) were observed. Furthermore, a significant positive correlation between GC 3S and ENC values (r = 0.424, P < 0.001) as well as a significantly negative correlation between AU 3S and In addition to correlation analysis, linear regression analysis was also performed to determine correlations between f' 1 and f' 2 and nucleotide constraints of MARV genomes (Table 4 ). In agreement to above findings, significant correlations were observed between both axes and compositional quantities indicating that mutation pressure has played a major role in shaping the dynamics of codon usage patterns within MARV genomes.

A neutrality plot was constructed between P 12 (GC 1,2S ) and P 3 (GC 3S ) values to determine the extent of variation between mutation pressure and natural selection (Fig. 4) . A significant positive correlation (r = 0.747, P < 0.01) was found between P 12 and P 3 values with a correlation coefficient of 0.926 ± 0.067, suggesting that the effect of 

To determine the potential influence of natural selection, linear regression analysis was performed between General average hydropathicity (GRAVY) and aromaticity (ARO) values and the f' 1 , f' 2 , ENC, GC, and GC 3S values to investigate the influence of natural selection on MARV codon usage patterns. The correlations of both GRAVY and ARO with f' 1 were non-significant, whereas GRAVY and ARO showed significant positive and negative correlations with f' 2 , respectively. Furthermore, it was found that GRAVY had significantly positive correlations with ENC, GC 3S , and GC values, whereas ARO had significant negative correlations with GC 3S and GC and non-significant correlations with the ENC value ( Table 5 ). The non-significant correlation of both GRAVY and ARO with f' 1 , which accounts for 65.55 % of the total variation, shows that natural selection has contributed to some extent; however, it is not the most substantial influencing factor on MARV codon usage patterns.

In order to determine codon usage optimization and adaptation of MARV to its hosts, CAI analysis was performed. A mean CAI of 0.712 was obtained for MARV genes in relation to H. sapiens, while a mean CAI of 0.534 was obtained in relation to R. aegyptiacus (Additional file 3: Table S2 ). There was a trend for a lower CAI values for MARV in relation to R. aegyptiacus, with the consequent lower efficiency of protein synthesis in R. aegyptiacus. Furthermore, correlation was investigated between CAI and ENC values to examine the relative influence of mutation pressure and natural selection. The CAI values of MARV genes in relation to H. sapiens and R. aegyptiacus were found to be 

To study the possible effect of dinucleotides on codon usage, we calculated the relative abundances of the 16 dinucleotides from the MARV coding sequences. The occurrence of dinucleotides was found to be non-random, and only CpU was present at the expected frequencies (i.e., 1.0) ( Table 6 ). Furthermore, only CpA was overrepresented and showed marginal over-representation (1.23 ± 0.01). CpG (mean ± SD = 0.51 ± 0.01) and GpC (mean ± SD = 0.90 ± 0.03) were both under-represented.

The analysis of RSCU values of both CpG-containing codons (CCG, GCG, UCG, ACG, CGC, CGG, CGU, and CGA) and GpC-containing codons (GCU, GCC, GCA, UGC, AGC, and GGC) showed that all codons were also under-represented (RSCU < 0.6) and were not preferred codons for their respective amino acids ( Table 2 ). Similar to CpG and GpC, the relative abundance of UpA also deviated from the "normal range" (mean ± SD = 0.69 ± 0.01) and was under-represented. Except for UUA (RSCU = 1.72), which is a preferred codon for the amino acid leucine, the remaining five UpA containing codons (CUA, GUA, UAU, UAC, and AUA) were under-represented (RSCU < 0.6) and not preferred codons. Five (UCA, ACA, GCA, CAA, and CAU) out of eight codons that contain CpA (CCA, CAG, and CAC) were also over-represented and preferred codons compared with the rest of the codons for their respective amino acids (Table 2) . Correlation between the relative abundance of dinucleotides with the f' 1 and f' 2 was also investigated. Fourteen and 12 out of 16 dinucleotides Table 3 Summary of correlation analysis between nucleotide constraints in MARV genomes The numbers in the each column represent correlation coefficient "r" values, which are calculated in each correlation analysis NS non-significant (P > 0.05) *represents 0.01 < P < 0.05 **represents P < 0.01 The numbers in the each column represents correlation coefficient "r" values, which are calculated in each correlation analysis NS non-significant (P > 0.05) *represents 0.01 < P < 0.05 **represents P < 0.01 (Table 6 ).

To determine the influence of host on MARV codon usage patterns, the codon usage of MARV isolates was compared with that of its two hosts, H. sapiens and R. aegyptiacus, via comparison of RSCU values. The results showed that the codon usage patterns and selection of preferred codons in MARV genomes is antagonist to both H. sapiens and R. aegyptiacus for majority of codons ( Table 2 ). The only exception was codon AGA, which was a preferred codon for the amino acid arginine in MARV and H. sapiens but not in R. aegyptiacus. Moreover, analysis of RSCU values in a heatmap also showed that the MARV RSCU values did not cluster along any of its hosts RSCU values (Additional file 4: Figure S2 ).

To determine how the overall codon usages of MARV's hosts have contributed to evolution of virus codon usage patterns, similarity index analysis was conducted. It was found that R. aegyptiacus exerted a more dominant effect on shaping MARV codon usage compared with that of H. sapiens, as the similarity index was found to be higher in R. aegyptiacus (Fig. 5 ).

In the present study, we analyzed synonymous codon usage in coding sequences from 63 MARV genomes to understand its molecular evolution under the influence of multiple viral, host, and environmental factors. It has previously been shown that codon usage bias, or preference for one type of codon over another, can be greatly influenced by overall genomic composition [42] . Nucleotide composition analysis showed that A and U nucleotides constitute the majority of overall nucleotide composition in MARV genomes ( Table 1) . The RSCU analysis also showed that MARV genomes exhibit greater codon usage bias toward A-and U-ended codons (Table 2) . Therefore, once it is established that there is codon bias toward A-and U-ended codons in MARV genomes, we next determined the extent of this bias within and in between different MARV isolates. This was accomplished by ENC analyses. In the case of MARV, the mean ENC value was found to be 54.2 in MARV coding sequences, which indicates slightly biased, relatively stable and conserved ARO aromaticity; NS non-significant (P > 0.05) *represents 0.01 < P < 0.05 **represents P < 0.01 Table 6 Summary of correlation analysis between the first two principal axes and relative abundance of dinucleotides in MARV genomes [46] , and West Nile virus (ENC: 53.81) [11] . A possible explanation given for this is that the low codon bias of RNA viruses might be advantageous for efficient replication in host cells by reducing the synthesis machinery competition between the virus and host with potentially distinct codon preferences. Whether the same holds true for MARVs as well, warrants further investigations. However, this can be attributed to the fact that MARVs maintain low yet surviving replication rate within in its natural host, R. aegyptiacus without causing any disease conditions [47] . Therefore, it seems that evolution of low codon bias within MARV coding sequences have enabled it to successfully maintain its survival cycle within both of its hosts each of which possess distinct codon usage preferences from that of MARV (Table 2) .

Considering the multivariate nature of codon usage, we next performed PCA analysis on RSCU values to determine the trends of codon usage variations that showed that f' 1 accounted for the major portion of codon usage variation followed by f' 2 . Moreover, MARV isolates formed three separate PCA clusters following distribution on principal axes. The clustering of diverse MARV lineages that are separated by thousands of miles within a single cluster as well as clustering of closely related lineages into different clusters highlights an important role of the mobility of MARV's natural host, R. aegyptiacus. This also indicates that the isolates of MARV might have independently evolved in three clusters after diverging from a common ancestor that potentially originated from DRC, based on inclusion of DRC isolates in all three clusters. Moreover, it appears that the geographical diversity and associated factors, such as presence of natural host within region of infection, climatic features, and host susceptibility, have contributed to shaping codon usage in MARV genomes.

Our initial analysis indicated influence of nucleotide constraints on MARV codon usage patterns. However, it has previously been shown that, although overall RSCU could reveal the codon usage pattern for genomes, it may hide the codon usage variation among different genes in a genome [48] , thereby indicating that composition frequencies of nucleotides are not always the only factor associated with codon usage patterns. The ENC plot is widely used to determine codon usage variation among genes in different organisms. It has been postulated that the ENC-plot of genes for which codon choice is constrained only by compositional constraints or mutation pressure will lie on the continuous curve of the predicted ENC values [32] . When the ENC and GC 3S values of MARVs were plotted, it was found that although none of the isolates fell on the expected continuous curve but clustered closely below the curve, thereby showing major influence of mutation pressure on MARV codon usage patterns and of natural selection to some extent. Besides that, it has also been reported previously that both mutation pressure and natural selection can influence the overall ENC and it might not be a robust index to show the relative contribution of mutation and selection on structuring codon usage patterns. Moreover, the codon usage bias of base composition of the genes of a species with A/U biased genomes will behave differentially than those species with G/C biased genomes, as such, ENC-GC 3S plot might be potentially misleading. In contrast, CAI is suggested as the most robust index for showing the influence of natural selection on codon usage patterns of such genes [27, 30, 31, 49] . CAI is regarded as a measure of gene expression and can be used to assess the adaptation of viral genes to their hosts. It has been postulated that the highly expressed genes exhibit a stronger bias for particular codons compared with genes that are less expressed. Compared with ENC, which is another way of calculating codon usage bias and measures deviation from a uniform bias (null hypothesis), CAI measures the deviation of a given protein coding gene sequence with respect to a reference set of genes [27] . If CAI value is high, then codon usage bias is extremely high and the influence of natural selection is prevailing. CAI values were calculated for MARV genes separately for both of its hosts. MARV genes showed higher CAI values for H. sapiens (0.712) as compared to R. aegyptiacus (0.534) indicating that natural selection from both hosts have influenced the codon usage patterns of MARV. Furthermore, comparative analysis of CAI values between MARV and its hosts suggests that MARV genes have optimized their codon usage patterns to utilize the translational resources of H. sapiens more efficiently than that of R. aegyptiacus. A higher CAI values of MARV genes for H. sapiens represents an interesting evolutionary step which might have supported MARV to turn out to be a highly pathogenic virus for H. sapiens and at the same time remaining completely harmless for its natural host, R. aegyptiacus.

Among multiple influencing factors, mutation pressure and natural selection are considered the two major factors that shape codon usage patterns [50] . A general mutation pressure that affects the whole genome would certainly account for the majority of the codon usage among certain RNA viruses [42] . The ENC and CAI analyses highlighted the influence of both mutation pressure and natural selection on codon usage patterns of MARV genes. In order to determine the share of each factor on evolution of MARV codon usage patterns, a neutrality plot analysis was performed which showed that influence of mutation pressure dominates over natural selection. Furthermore, we also examined the influence of mutation pressure on MARV codon usage via correlation and linear regression analyses between different nucleotide compositional constraints, ENC, and principal axes. Strong and significant correlations were observed, which indicates a dominant influence of mutation pressure. This was further supported when these indices were plotted against the first two principal axes via PCA, and significant strong correlations were observed. However, in the case of MARV genomes, involvement of factors other than mutation pressure such as natural selection cannot be ignored because nucleotide base compositions showed variation, distribution of MARV isolates were although close to but still below the expected curve on ENC plot, and there was a weak codon bias. A weak codon usage bias may be caused by natural selection when the viruses try to adapt to the host cell [51] [52] [53] . It has been suggested that, if synonymous codon usage bias is affected by mutation pressure alone, then the frequency of nucleotides A and U/T should be equal to that of C and G at the synonymous codon third position [53] . To test this phenomena in MARV, PR2 biasness plot analysis was performed on four fold degenerate codons. The occurrence frequency of AU and GC nucleotides at the synonymous codon third position was not found to be equal and AU biased preference was observed in four fold degenerate codons of MARV genes which indicates the potential influence of natural selection on codon usage patterns of MARV genes. In addition to this, correlation analysis between principal axes and GRAVY and ARO also revealed that, although natural selection has influenced MARV codon usage patterns to some extent, it is much weaker compared with mutation pressure.

Dinucleotide abundance has been reported to influence overall codon usage bias in several organisms, including DNA and RNA viruses [37, 54, 55] . Toll-like receptor 9 (TLR9), which is a type of intracellular pattern recognition receptor (PRR), recognizes unmethylated CpGs, which leads to activation of several immune response pathways [56] . The vertebrate immune system relies on unmethylated CpG recognition in DNA molecules as a signature of infection, and CpG underrepresentation in RNA viruses is exclusively observed in vertebrate viruses; therefore, it is reasonable to suggest that a TLR9-like mechanism exists in the vertebrate immune system that recognizes CpGs when in an RNA context (such as in the genomes of RNA viruses) and triggers immune responses [57] . In contrast to the CpG usage of + ssRNA viruses that are greatly influenced by their hosts and because of which + ssRNA viruses mimic their hosts' CpG usage, -ssRNA viruses do not produce DNA intermediates during the replication of their genome. As a result, CpGs are under-represented, independent of the infected host or their phylogenetic relationship. The under-representation of CpG in -ssRNA viruses is therefore due to the U/A mutation bias in overall genomic composition that further indicates a dominating effect of mutation pressure [54, 58] . In the case of MARV, none of the dinucleotides were found at the expected frequencies and were also markedly underrepresented. As inferred from the RSCU analysis, codons containing CpG and UpA dinucleotides were also underrepresented and were not preferred codons for their respective amino acids within MARV genomes. These results collectively indicate that, although dinucleotide representation has some influence over the codon usage of MARVs, the overall influence is not strong because of exceptions among dinucleotide frequencies and selection of preferred codons in MARV genomes.

It has been previously shown that, among many other factors, the codon usage patterns of viruses are also affected by its hosts [59] . For example, the codon usage pattern of poliovirus is reported to be mostly coincident with that of its host [60] , whereas the codon usage pattern of hepatitis A was reported to be antagonistic to that of its host [61] . Alternatively, some viruses exhibit a mix of both coincidence and antagonism, such as the HCV [46] , enterovirus 71 [9] , and CHIKV [44] . However, MARV showed almost complete antagonism to both of its hosts, as inferred from the RSCU analysis with the exception for a common preferred codon for arginine between MARV and H. sapiens (Table 2) . A recent codon usage analysis in the EBOV, which is also from the same family as MARV and spread via bats, reported similar antagonism of RSCU toward its host, H. sapiens [43] . It has been proposed that the coincident portions of codon usage among viruses and their hosts could enable the corresponding amino acids to be efficiently translated, whereas the antagonistic portions of codon usage may enable viral proteins to be properly folded, although the translation efficiency of the corresponding amino acids might decrease [46] . Whether the same holds true for MARV warrants further investigations. In addition to the RSCU comparison analysis, we also performed a similarity index analysis to determine which of the MARV hosts have a dominant influence over its RSCU patterns. Evidence of selection pressure from both hosts was detected which is in agreement to the CAI analysis; however, level of pressure was significantly different. Compared with H. sapiens, R. aegyptiacus have a more profound effect on shaping MARV RSCU patterns, as inferred from the similarity index analysis. As R. aegyptiacus is consider as the natural reservoir and host for MARV, it makes sense that virus has evolved its genomic features to a stable level in order to better adapt to its primary host's environment. It has also been recently suggested that flight, a factor common to all bats but to no other mammals, provides an intensive selective force for coexistence with viral parasites through a daily cycle that elevates metabolism and body temperature analogous to the febrile response in other mammals. On an evolutionary scale, this hostvirus interaction might have resulted in the large diversity of zoonotic viruses in bats, possibly through bat viruses adapting to be more tolerant of the fever response and less virulent to their natural hosts [47] .

The aetiology of an unknown human infection transmitted by monkeys (preliminary communication)

Marburg virus infection detected in a common African bat

Marburgvirus genomics and association with a large hemorrhagic fever outbreak in Angola

Marburg virus, a filovirus: messenger RNAs, gene order, and regulatory elements of the replication cycle

Forty-five years of Marburg virus research

Codon catalog usage and the genome hypothesis

Variation in G + C-content and codon choice: differences among synonymous codon groups in vertebrate genes

Analysis of synonymous codon usage in SARS Coronavirus and other viruses in the Nidovirales

The characteristics of the synonymous codon usage in enterovirus 71 virus and the effects of host on the virus in codon usage pattern

Analysis of synonymous codon usage in dengue viruses

A detailed comparative analysis on the overall codon usage patterns in West Nile virus

Codon usage patterns in Escherichia coli, Bacillus subtilis, Saccharomyces cerevisiae, Schizosaccharomyces pombe, Drosophila melanogaster and Homo sapiens

a review of the considerable within-species diversity

Analysis of synonymous codon usage in classical swine fever virus. Virus genes

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

RDP3: a flexible and fast computer program for analyzing recombination

Possible emergence of new geminiviruses by frequent recombination

RDP: detection of recombination amongst aligned sequences

Analyzing the mosaic structure of genes

Evaluation of methods for detecting recombination from DNA sequences: computer simulations

A modified bootscan algorithm for automated identification of recombinant sequences and recombination breakpoints

Sister-scanning: a Monte Carlo procedure for assessing signals in recombinant sequences

An exact nonparametric method for inferring mosaic structure in sequence triplets

Phylogenetic evidence for recombination in dengue virus

Codon usage in regulatory genes in Escherichia coli does not reflect selection for 'rare' codons

An evolutionary perspective on synonymous codon usage in unicellular organisms

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

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

Codon usage tabulated from international DNA sequence databases: status for the year 2000

Codon usage in twelve species of Drosophila

Analysis of synonymous codon usage patterns in duck hepatitis A virus: a comparison on the roles of mutual pressure and natural selection

The 'effective number of codons' used in a gene

An evaluation of measures of synonymous codon usage bias

Intrastrand parity rules of DNA base composition and usage biases of synonymous codons

Translation-coupled violation of Parity Rule 2 in human genes is not the cause of heterogeneity of the DNA G + C content of third codon position

The distribution of synonymous codon choice in the translation initiation region of dengue virus

Dinucleotide relative abundance extremes: a genomic signature

Does recombination improve selection on codon usage? Lessons from nematode and fly complete genomes

Codon usage bias: causative factors, quantification methods and genome-wide patterns: with emphasis on insect genomes

Extensive homologous recombination in classical swine fever virus: A re-evaluation of homologous recombination events in the strain AF407339

Codon usage bias and recombination events for neuraminidase and hemagglutinin genes in Chinese isolates of influenza A virus subtype H9N2

The extent of codon usage bias in human RNA viruses and its evolutionary origin

Genome-wide analysis of codon usage bias in Ebolavirus

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

Analysis of codon usage in bovine viral diarrhea virus

The characteristic of codon usage pattern and its evolution of hepatitis C virus

Bat flight and zoonotic viruses

Synonymous codon usage analysis of thirty two mycobacteriophage genomes

A comparison of synonymous codon usage bias patterns in DNA and RNA virus genomes: quantifying the relative importance of mutational pressure and natural selection

GC3 biology in corn, rice, sorghum and other grasses

Selective pressure dominates the synonymous codon usage in parvoviridae

Optimization of codon usage of poxvirus genes allows for improved transient expression in mammalian cells

Analysis of synonymous codon usage patterns in torque teno sus virus 1 (TTSuV1)

CpG usage in RNA viruses: data and hypotheses

Second codon positions of genes and the secondary structures of proteins. Relationships and implications for the origin of the genetic code

Clinical application of CpG-, non-CpG-, and antisense oligodeoxynucleotides as immunomodulators

Virushost coevolution: common patterns of nucleotide motif usage in Flaviviridae and their hosts

Patterns of evolution and host gene mimicry in influenza and other RNA viruses

Heterogeneity in codon usages of sobemovirus genes

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

Genome variability and capsid structural constraints of hepatitis a virus

In summary, this study showed that overall codon usage bias within MARVs is slightly biased, and the major factor that has contributed to shaping codon usage is mutation pressure followed by influence of hosts. In addition, natural selection, environment, geographical conditions, and dinucleotides have also been determined to influence codon usage. The evolution of MARV probably reflects a dynamic process of mutation and natural selection to adapt its codon usage to different environments and hosts.

The data sets supporting the results of this article are available in the Dryad digital repository http://dx.doi.org/ 10.5061/dryad.3hc5t.

Additional file 1: Table S1 . Demographics of MARV genomes that were analyzed in the present study. (DOCX 23 kb) Additional file 2: Figure S1 .