Horizontal gene transfer and nucleotide compositional anomaly in large DNA viruses

Compositionally anomalous genes in LDVs

We next analyzed compositional heterogeneities across different genes encoded in individual LDV genomes. Genes with nucleotide compositions substantially deviating from the average in a genome tend to be under distinct selection pressures or have particular evolutionary histories. We denote such genes as compositionally anomalous (cA) genes. To identify cA genes with a robust statistical support, we used a method [26] based on Markov modeling and Bayesian inference originally developed for the identification of horizontally transferred genes within prokaryotic genomes. Composition-based approaches have a clear limitation in detecting HGTs; they detect only a subset of HGTs (i.e. recent horizontal gene acquisitions by a recipient genome from a donor genome with nucleotide compositions that are different from the recipient genome). Our aim here is to examine the nature of cA genes in LDVs in the light of previous observations of cA genes for prokaryotic genomes.

We identified cA genes in all the analyzed LDV genomes. In many LDV genomes, the cA genes were more A+T rich than the remaining genes (for 43 out of 67 LDVs; binomial test, p < 0.05). The proportion of cA genes per genome was highly variable across the 67 LDVs, ranging from 1% to 53% (Table 1, Fig. 4). In contrast, the proportion obtained by the same method is less variable for prokaryotes (0% to 25%) [9]. Herpesviruses exhibited the widest range of cA gene proportion (6% to 53%). All but two of the analyzed NCLDVs, including mimivirus, had a relatively low level of cA proportion (1%–17%). Two exceptions were Paramecium bursaria chlorella virus 1 (28%) and molluscum contagiosum virus (31%). Phages showed a relatively low proportion (4% to 21%). The cA gene proportions for other viruses were 5% to 26% for four baculoviruses, 39% for nimavirus and 24% for Heliothis zea virus 1. The very high cA proportions in some betaherpesviruses appear to be linked with genomic islands harboring betaherpesvirus-specific genes, which we will discuss below (in the "Co-localization of cA genes" section).

Nakamura et al. reported a significant positive correlation between the cA gene proportion and the total number of genes in a genome for prokaryotes [26]; larger prokaryotic genomes tend to more rapidly acquire foreign genes (of "anomalous nucleotide compositions") than smaller ones. In contrast, we found no such correlation for the 67 LDVs (R² = 0.03) (Additional file 3).

We observed a weak but significant positive correlation (R² = 0.53, p < 0.005) between the proportion of cA genes and genomic G+C content (Additional file 4). This correlation was not due to a bias induced by a given viral family or subfamily. For example, when we excluded herpesviruses, many of which exhibit both a high cA gene proportion and a high G+C content, the correlation was reduced but still remained significant (R² = 0.24, p ≤ 0.005). Such a significant correlation was observed also within a single subfamily (i.e. alphaherpesviruses, R² = 0.66, p ≤ 0.005). A potential variation in gene prediction quality (which could be dependent on genomic G+C compositions) does not appear to explain these correlations (Monier et al., unpublished data). We also examined a possible relationship between the cA gene proportion and chromosome topology. There was no difference in the cA gene proportions between linear LDV genomes (15.5% on average) and circular LDV genomes (17.3% on average; t-test, p = 0.63).

Our global genomic signature analysis revealed remarkably different nucleotide compositions between LDVs and their hosts. To investigate further the discrepancy of nucleotide compositions, we examined if the identified cA genes are enriched in "host-like" sequences in terms of their nucleotide compositions. We determined the genomic signature distances (di-nucleotides for word length) from individual cA genes to the viral genome where they originate as well as to the host genomic sequences. As shown in the Additional file 5, we observed a marked enrichment of genes with host-like signatures among the identified cA genes (19%) relative to non-cA genes (3.3%). The host-like signatures in these cA genes could be due to a large variation in their nucleotide compositions, being unrelated to possible host origins for some of the cA genes. To test this hypothesis, we performed the same analysis by randomizing the pairings between viruses and their hosts. A comparable fraction (23% versus 19%) of the cA genes indeed exhibited smaller distances to the genomic signatures of randomly chosen hosts than to the viral genomes where they originate. Thus the enrichment of genes with host-like signatures in cA genes may not suggest their horizontal acquisitions from hosts.

Correlation between cA genes and replication-associated strand asymmetry

Genome sequences of many prokaryotes and vertebrates [33, 34] show strong strand asymmetry in nucleotide composition between the leading and the lagging strands. For most bacteria, compositional strand asymmetry is characterized by an excess of G and T bases in the leading strands relative to the lagging strands [33, 35, 36]. A substantial part of compositional strand asymmetry is independent of gene distribution in the two distinct DNA strands, and is probably due to mutational biases linked to asymmetric mechanisms of DNA replication [37, 38]. Compositional strand asymmetry spanning large genomic segments has been described also for some large DNA viruses, including mimivirus [5, 39] and several herpesviruses [40, 41]. Such replication-associated strand bias can potentially result in two classes of genes with distinct nucleotide compositions, depending on which strand they are located. We thus examined if there is a correlation between the distribution of cA genes and compositional strand asymmetry.

We first generated cumulative GT-excess plots [42] for all of the 67 LDV genomes to assess their compositional strand asymmetry. LDVs presented various GT-excess patterns. Several LDVs, mostly phages, showed strong global strand asymmetry with a monotonous increase (or decrease) of the cumulative GT-excess curve along their entire genome. Some LDVs including mimivirus exhibited a few peaks in their GT-excess profiles. For other LDVs, the GT-excess curves locally fluctuated with no long genomic segments exhibiting a consistent compositional asymmetry. We selected fifteen LDVs presenting long (>10 kb) genomic segments with uniform compositional asymmetry (Additional file 6). After identifying the genomic coordinates where the sign of the nucleotide-skew (G+T versus C+A) changes, we split the genomic sequences into sub-strands. Those sub-strands were classified into two classes: class I consists of sub-strands having a positive nucleotide-skew (i.e. G+T% > C+A%), and class II with a negative skew (i.e. G+T% < C+A%). Genes were then classified into either class I or II according to the sub-strand they originated from. For instance, the plus strand (according to the GenBank entry) of the Pseudomonas phage phiKZ genome is G and T rich along its entire length, thus being classified as class I. This class I strand contains 229 genes (≥ 300 nt), of which 22 were detected as cA genes. The complementary strand is classified as a class II and contains 49 genes, of which 5 were cA genes. No significant correlation was found between the distribution of cA genes and the compositional strand asymmetry (Fisher's exact test, p-value = 0.99). The mimivirus genome shows a clear switch of nucleotide-skew (G+T versus C+A) at position 380,000 nt as previously noted [5]. Thus a part of the plus strand (from 0 to 380,000 nt) and a part of the complementary strand (from 380,000 to 1,118,404 nt) constitute the class I strands. Class II is represented by the remaining strands (Additional file 6). Again, no significant correlation was found between the distribution of cA genes and the compositional strand asymmetry (p-value = 0.63). Of the fifteen LDVs that we analyzed, only three showed a significant correlation between cA gene distribution and strand asymmetries (p-values < 0.05). The three viruses are the fowlpox virus and two strains of the human herpesvirus 5 (the wild type strain merlin and the laboratory strain AD169). These results suggest that replication-associated compositional strand asymmetry accounts for only a small part of the nucleotide compositional biases of the cA genes observed in LDVs.

Co-localization of cA genes

In prokaryotes, horizontally transferred genes are often clustered in the chromosome and make genomic islands [26, 43, 44]. We investigated co-localization of the cA genes along the chromosomes of the LDVs using a Monte-Carlo simulation (see Methods). None of the genomes of entomopoxviruses, baculoviruses, iridoviruses, mimivirus, nimavirus or phycodnaviruses exhibited a significant propensity for cA genes to be clustered along the genome. Of the 67 LDVs, only 12 presented a significant frequency of pairs of co-localized cA genes (p ≤ 0.01; Table 1). Using the same Monte-Carlo simulation, we identified significantly large sized cA gene islands (ranging from 3 up to 32 neighboring cA genes) in 11 LDV genomes (p ≤ 0.01; Additional file 7). Those were 6 betaherpesviruses, 1 alphaherpesvirus, 2 chordopoxviruses, African swine fever virus, and the cyanophage P-SSM4. Betaherpesviruses showing cA islands of significant sizes include murid herpesvirus 2 with an extravagant 53% of cA gene proportion. The abundance of cA gene islands in betaherpesviruses parallels the mosaicity in G+C composition along their genomes [17]. It is known that betaherpesvirus genomes are composed of core genes located in the middle of their genomes and shared by different subfamilies of Herpesviridae, plus betaherpesvirus-specific genes located at the extremities of the genomes [45]. As shown in Fig. 5, many of the beta-family genes were identified as cA genes making genomic islands.

Correlation between cA genes and functional properties

Functional constraint on specific genomic regions can be a cause of atypical nucleotide compositions of genes. We selected five genomic regions to examine if there were correlations between cA genes and certain of their functional attributes. For the first three cases (gamma- and betaherpseviruses), the correlation between the nucleotide compositional biases and gene functions has already been described, and could serve as a positive control for our method. For the last two cases (Emiliania huxleyi virus-86 and mimivirus), no such correlation has been previously reported.

Human herpesvirus 4 (Epstein-Barr virus) of the gammaherpesvirus subfamily has two different life cycle modes: latent and productive. Nine genes are known to be expressed during latency [46]. Those genes are hardly expressed during the productive mode. They are highly A+T rich compared to the remaining genes. Karlin et al. hypothesized that the unique codon usage of those genes helps to minimize the competition with host genes for various host resources during latency [19]. We found that eight of the nine latency genes were indeed detected as cA genes (89%). In contrast, the cA gene proportion was significantly lower for the remaining genes (29/81, 36%, p = 0.003).

The immediate-early transcription region (≈10-kb) of the murid herpesvirus 1 (murine cytomegalovirus, betaherpesvirus subfamily) genome is known to be CpG suppressed [47]. Experimental data suggest that the methylation of the cytosines in CpG dinucleotides in the enhancer/promoter of these regions has a regulatory role in gene expression [48]. The CpG suppressed region of the murid herpesvirus 1 genome encodes 10 genes. We found that the cA gene proportion in this region (6/10 genes, 60%) was twice as important as in the remaining regions (46/146, 32%) of the genome, though their difference might be due to chance (p = 0.085).

The murid herpesvirus 1 genome possesses three (19-kb, 10-kb and 17-kb) regions that are A+T rich relative to other parts of the genome [47]. Those regions are enriched in genes encoding membrane glycoprotein (e.g. m02 and m145 families). Some of those proteins were suggested as responsible for the evasion from natural killer cell-mediated immune surveillance through their interaction with inhibitory natural killer cell receptors [49, 50]. Of 37 genes within the A+T rich regions, 67% (24 genes) were detected as cA genes, which is significantly higher than the cA gene proportion for the remaining genes (23.5%, 28/119 genes, p < 10^-5). It is uncertain if the increased A+T levels of these regions have functional roles. However, it should be noted that these A+T rich regions contain five of the seven genes exhibiting significant sequence polymorphisms between strains of murid herpesvirus 1 [51].

Emiliania huxleyi virus-86 (EhV-86) of the Phycodnaviridae family exhibits two distinct transcription phases during its lytic infection to the host alga, E. huxleyi [52]. The primary phase is characterized by the transcriptions of a group of genes by 1 hour post-infection. The secondary phase involves the transcriptions of other genes between 2 and 4 hours post-infection. Respectively thirty-eight and 253 genes in our data set (ORF length ≥ 100 codons) correspond to the primary and the secondary transcription phases. We found a significantly greater fraction of cA genes in the primary phase (8/38, 21%) than in the secondary phase genes (19/253, 7.5%; p = 0.014). Genes expressed during the primary phase map in the 104-kb central genomic region (bases 200,000 to 304,000), which shows a similar G+C content as the rest of the genome [53]. It has been suggested that promoter-like elements (family A repeats) uniquely found in this region control this early expression pattern of EhV-86 [7]. The functions of the early expression are unknown as most of the transcribed genes lack detectable homologs in the databases. Allen et al. postulated that an ancestral EhV acquired the 104-kb region by horizontal transfer [53].

Mimivirus has unusually well conserved promoter-like AAAATTGA motifs in the upstream regions of half of its predicted genes in the genome [39]. Based on the putative associated gene functions, Suhre et al. predicted that these promoter-like elements regulate gene expression during the early stages of the viral infection, whilst most of the genes contributing to the virus particle are devoid of this motif [54]. Of 402 genes with the promoter-like motifs in our data set, 43 (10.6%) were detected as cA genes. Of 508 genes lacking the motifs, 40 (7.9%) were detected as cA genes. The difference is not statistically significant (p = 0.16).

In summary, we found a significant correlation between cA genes and previously described functional categories of genes for three of the five cases. It should be noted that in the case of EhV, no relationship between expression timing and gene nucleotide composition was previously reported. This suggests the possible use of nucleotide composition analysis to predict expression patterns of viral genes. These results indicate that anomalous nucleotide compositions of some of the cA genes can be due to functional constraints, although the distinction between functional constrains and horizontal transfer events is generally difficult given the known bias in functions of horizontally transferred genes in prokaryotes [24, 26].

Additional criteria 1

- From independent evidences, one can chose a group of sequences, O1, representing an outgroup of I.

- From independent evidences, one can chose a group of sequences, O2, representing an outgroup of I and O1. [Consequently, the tree topology will become (O2, (O1, (V, I))), when rooted by O2.]

Additional criteria 2

- There are only a few species and/or genera in I.

- Other sequences outside the (V, I) group are from a wide variety of genera.

- The branches from the root of (V, I) to V and I are relatively shorter than the branches from the root of (V, I) to other peripheral nodes.