High-Throughput Analysis of Human Cytomegalovirus Genome Diversity Highlights the Widespread Occurrence of Gene-Disrupting Mutations and Pervasive Recombination

High-throughput sequencing of complete genomes from clinical HCMV isolates.To efficiently characterize the genetic diversity of a large set of complete genomic sequences derived from clinical HCMV isolates, we recently described a method that combines limited virus culturing and virion DNA purification with multiple-displacement amplification and NGS (36). We showed that this procedure was able to generate highly pure viral DNA suitable for NGS analysis and validated that strain consensus sequences were representative of the original virus populations in the clinical isolates. Here, we implemented this method to characterize complete genomes in a group of 100 clinical HCMV isolates (see Table S1 in the supplemental material). These isolates were collected from Belgian (n = 94) and Czech (n = 6) individuals infected with HCMV, including healthy adults, immunosuppressed patients, and congenitally infected infants.

DNA sequence reads were generated by using a combination of the 454 GS FLX and Illumina NGS platforms. In Table S1 in the supplemental material, some basic genome assembly statistics are listed. We successfully derived a full-genome consensus sequence for 96/100 strains, with the average read depth ranging from 35 to 3,315 (median, 2,031) and the proportion of reads mapping to the consensus ranging from 1% to 98% (median, 89%). While eight isolates were successfully sequenced at read percentages of <10%, sample purity was >50% in the majority of isolates (77/96). For three strains, we could not determine the full-genome consensus because coverage was too low or unevenly distributed, leaving too many sequence gaps for finishing through Sanger sequencing. For one strain, coverage was adequate, but the isolate clearly consisted of multiple genome variants, and a single, meaningful consensus sequence could not be obtained, nor was it possible to segregate the constituting variants, since NGS data provide no connection between variants at different variable loci. The predominance of a single genome variant along the entire genome for all other strains suggests that these sequences constitute a contiguous genome.

As reported previously, genome assembly consisted of de novo assembly, scaffolding on a reference sequence, and subsequent iterative mapping of NGS reads on the genome scaffold (36, 62). Terminal repeats were omitted from the scaffold since these repeats are identical to the internal repeats. This approach was mostly successful, but assembly of the internal repeat regions (IRL/IRS) generally required additional consideration, since de novo contigs tended to break at these regions. Because of the high sequence variability of the internal repeat regions (Fig. 1), reference-assisted iterative mapping usually did not solve this issue. Therefore, contigs were split, sequence reads were assembled on both contigs separately, and contigs were enlarged by the 50% read overhangs at contig ends. This process was reiterated until contigs could be joined. The transition between unique and repeat regions needed to be determined manually to correctly add the terminal repeat sequences at both genome ends. Transitions can be recognized from NGS read assemblies through the concomitant mapping of transition-crossing reads from different genome isomers with inverted UL and US directions. While the position of the US-IRS junction is stable, presumably because it is located inside the TRS1/IRS1 reading frame, the UL-IRL junction is located in a noncoding region and is positioned differently in separate strains. Strain BE/5/2010 had an unusual layout for the IRL/IRS region, with the IRL repeat starting after the IRS repeat and being completely encompassed by it.

• Open in new tab
• Download powerpoint

FIG 1

Diversity and evolution of the HCMV genome. Shown is an overview of genetic diversity and evolutionary pressure along the HCMV genome. The genome is divided into four panels. Each panel consists of four separate tracks. In the top track, nucleotide diversity is calculated in a sliding window of 500 nt with a step size of 100 nt. As only ungapped residues are included in each window, the distance between two data points may vary in areas with many indels. In the second track, the extents of recombination and positive selection are assessed for each gene. Displayed above the center of the appropriate gene (bottom track), dark and light blue bars represent recombination breakpoint density and the percentage of codons under positive selection, respectively. For optimal resolution, values were cut off at 10 breakpoints (brp)/kb and 4% codons under positive selection. Green and red bars in the third track indicate the genome positions of conserved and variable tandem repeats. The bottom track annotates genes and other genome elements in four layers. The first two layers show genes carried on the forward strand, and the last two layers show genes carried on the reverse strand. Spliced exons are connected with thin black lines. Genomic inverted repeats (TRL, IRL/IRS, and TRS) and long noncoding RNAs are represented in black; genes are colored on a scale from green to red, indicating the frequency of ORF-disrupting mutations in separate clinical isolates.

HCMV displays the highest level of genetic diversity of all human herpesviruses.The overall genetic diversity of the HCMV genome was assessed by aligning our 96 genomes with 28 previously reported sequences (see Table S1 in the supplemental material). This alignment contains 255,248 sites: 223,991 sites are without gaps, and 31,528 of these nongapped sites (14%) are polymorphic. The interstrain nucleotide diversity, π, was estimated to be 0.021, and the average number of nucleotide differences between two genomes is 4,734. This number is higher than the recent estimate of π of 0.015 for murine cytomegalovirus (MCMV), based on 11 complete genomes (63). There is a clear discrepancy with the values for intrastrain nucleotide diversity that were reported for congenitally infected infants (π = 0.18 to 0.25) (64). Care should be taken in directly comparing these estimates, as interstrain diversity is estimated from measuring polymorphisms in separate consensus sequences, while intrastrain diversity is derived from characterizing polymorphisms of a single virus population. However, the large discrepancy suggests that many of the variants that were identified in these intrahost populations are deleterious and are not passed on. This notion is supported by the apparent stability of genotype sequences in patients (30, 65–69).

To put these data in perspective toward the other human herpesviruses, we estimated the overall mean distance and transition/transversion ratio for alignments of all available complete genomic sequences of all nine human herpesviruses (Table 1). Apart from HCMV, only for HSV-1, VZV, and EBV were the numbers of complete genomic sequences adequate for a meaningful estimation of overall diversity. With >0.02 substitutions/site, versus <0.01 substitutions/site, HCMV stands out as being significantly more diverse than these alpha- and gammaherpesviruses (P = 0.012 for comparisons with HSV-1, VZV, and EBV only, or P = 2.5e−06 for comparisons with all human herpesviruses [determined by one-sample t test]). This is not an artifact of the higher number of sequences available for HCMV, since analysis of the overall mean distance in five random, separate subsets gave highly similar results (0.025 to 0.028 substitutions/site). In fact, considering the larger geographical diversity of HSV-1 and VZV strains included in this analysis, the current estimate for HCMV could even be too low (see “HCMV evolution has been shaped by pervasive recombination,” below). There is also a great discrepancy in the estimated transition/transversion ratios between different herpesvirus species. HCMV has a relatively high ratio of 2.53. This is probably a consequence of strong purifying selection removing transversions, which result in more nonsynonymous mutations (see “Positively selected residues provide a genetic fingerprint of the evolutionary arms race between virus and host,” below) (70).

In Fig. 1 (top), π is represented in a sliding window along the HCMV genome. This clearly delineates several diversity hot spots, isolated by long stretches of conserved sequence. To further analyze the heterogeneity of divergence, ratios of nonsynonymous substitutions per nonsynonymous site (dN) were calculated for all 170 genes as a measure of the divergence of the encoded proteins (see Table S4 in the supplemental material). A list of the most divergent genes is provided in Table 2. All 30 genes that have a dN value of >0.025 are listed, along with their gene family, conservation over herpesvirus subfamilies, their confirmed/proposed function(s) (71), and studies that have previously characterized and classified the diversity of these genes. Only 1 of these 30 genes is conserved in all mammalian herpesviruses (UL73, encoding glycoprotein N), while 3 others are conserved within the subfamily Betaherpesvirinae. For 4 of these 30 genes, no previous studies have analyzed diversity in clinical isolates. For 9 others (including 5 of the 6 most variable genes), analyses were limited to partial genomic sequences of <10 clinical isolates (7). This demonstrates the added value of a high-throughput, comprehensive analysis of divergent genome regions. Reciprocally, genes such as UL55 and UL75, encoding glycoproteins B and H, respectively, have been sequenced extensively because of their known functional roles, but they are not in this broad group of the most divergent genes. On the other end of the diversity spectrum, the 25 most conserved genes (dN of <0.002) are listed in Table 3. The majority of these genes (19/25) are conserved in all mammalian herpesviruses (14/25), in betaherpesviruses (3/25), or in beta- and gammaherpesviruses (2/25) (7).

Relationships of gene diversity with gene family, gene conservation across herpesviruses, and encoded function(s) were further explored and are visualized for gene families in Fig. 2 (top). Clear and significant differences in gene diversity between separate gene families and gene conservation groups were found (P = 4.8e−06 and P = 2.0e−08 [determined by KWt]) (Fig. 2). Because of the small number of genes in several gene families, pairwise comparisons were significant only for the RL11 family versus the US6 (P = 0.037 by WRSt), US12 (P = 0.0047 by WRSt), and US22 (P = 0.0044 by WRSt) gene families. The RL11 gene family indeed truly stands out for the high number of variable genes. While variability in RL11 genes was described previously, analyses of the most variable members, RL5A, RL6, RL12, and RL13, were limited to a comparison of seven strains (Table 2). Especially for the RL13 gene, it would be of interest to study the functional behavior of different variants, as this gene has been implicated as a growth temperance factor (72) and in immunomodulation (73). Considering the latter function of RL13, we found that the endocytic YxxL motif essential for internalization of IgGs was 100% conserved among clinical isolates. Cytomegalovirus-specific genes were significantly more diverse than genes conserved in all mammalian herpesviruses (P = 8.4e−08 by WRSt), betaherpesvirus genes (P = 0.0024 by WRSt), and genes conserved between beta- and gammaherpesviruses (P = 0.0025 [determined by WRSt]). The UL73 (encoding glycoprotein N) and UL74 (encoding glycoprotein O) genes are outliers within core and betaherpesvirus genes, respectively. Variability in these genes has been widely studied (Table 2). Similarly, we found statistically significant differences in the genetic diversity of genes classified according to the function of the encoded product (P = 9.8e−06 by KWt) (based on the functional classification reported in reference 71). Pairwise comparisons were significant only for immunomodulation genes versus genes encoding assembly (P = 0.0017 by WRSt), gene regulation (P = 0.0031 by WRSt), and replication (P = 0.011 by WRSt) functions. Generally, diverse genes are involved in interactions with the host (immunomodulation, entry, spread, cell tropism, and virion proteins, which include surface glycoproteins), while conserved genes perform core viral functions such as replication, assembly, modulation of the host cell cycle and proteins, gene regulation, cellular trafficking, nucleotide repair, virion stability, latency, and viral growth.

• Open in new tab
• Download powerpoint

FIG 2

Variability, recombination, and selection in HCMV gene families. Gene diversity (dN), recombination breakpoint density, and the percentages of codons under positive and negative selection are indicated for HCMV genes within gene families; each dot represents a gene. Only genes belonging to specific gene families are represented. Group averages are designated with horizontal lines.

Tandem repeats in the HCMV genome.Another important source of sequence variation is the heterogeneity in the copy numbers of adjacently repeated elements or tandem repeats (TRs), caused by recombination or strand slippage replication. Variation in TRs is associated with phenotypic variability, regulation of gene expression, and genetic evolvability in both prokaryotes and eukaryotes (74–76). Furthermore, several studies have found evidence that TR variations may impact strain functionality and pathogenicity in viruses (77–82). The presence of TRs in HCMV was described previously, and TR polymorphisms could be used as epidemiological markers to distinguish clinical isolates (83–85). A comparative analysis of TRs in several members of the family Herpesviridae, based on a single genome sequence for each species, found the highest TR content in the alphaherpesvirus pseudorabies virus (18% of total nucleotides), followed by HSV-1 (9%), EBV (7%), KSHV (4.5%), and VZV and HCMV (3%) (86). To assess the total set of repeats in HCMV genomes, we identified all homopolymers (repeats with a period size of 1 nt), microsatellites (period size of 2 to 9 nt), and minisatellites (period size of >9 nt) in the genome of reference strain Merlin and subsequently searched for orthologous TRs in 123 other HCMV genome sequences. The total set of identified repeats is reported in Table S5 in the supplemental material. In total, 779 TRs were found in the genome of strain Merlin, 23 of which are duplicated or triplicated in TRL and/or TRS inverted repeats. For 683/779 TRs (88%), an orthologous repeat could be found in the majority of the other HCMV strains. These 683 TRs constitute 3.9% of the total nucleotides and are annotated in the HCMV genome in Fig. 1. Only 51 of these orthologous repeat sets (7%) were classified as variable (<50% conservation of repeat sequence and copy number) (Fig. 1). While 81% of the total nucleotides are within protein-coding regions, only 65% of TR nucleotides are found inside genes (Fig. 3A). Reciprocally, there is a clear overrepresentation of TRs in noncoding regions, including the 4 long noncoding RNAs, introns, and intergenic regions. Likewise, the internal repeat regions (IRL-IRS) that make up only 1% of the trimmed Merlin genome contain 8% of TR nucleotides (Fig. 3B). Overall, the level of TR conservation is higher in coding than in noncoding regions (P = 9.9e−09 by Fisher's exact test [FEt]) (Fig. 3C). When different repeat types were analyzed separately, this held true for homopolymers (P = 1.3e−06 by FEt) and minisatellites (P = 0.017 by FEt) but not for microsatellites (P = 0.56 by FEt). Recently, TRs of HSV-1 were analyzed based on a collection of 26 complete genomes by using a similar approach (48). The authors of that study found 584 orthologous TRs in this data set, corresponding to 4.3 TRs/kb (5.4% of nucleotides), which is higher than the 2.9 TRs/kb (3.9% of nucleotides) that we found for HCMV. HSV-1 also has more variable TRs than does HCMV (17% versus 7%). HSV-1 shows a similar overrepresentation of TRs in noncoding regions and in the genomic inverted repeats. Likewise, most variable TRs are located in noncoding regions in both viruses. While HSV-1 has a higher proportion of genes containing TRs (92%, versus 79% for HCMV), this is caused mostly by the large number of homopolymers in HSV-1 genomes (3.4 TRs/kb, versus 1.2 TRs/kb for HCMV). When only the proportion of genes containing micro- and minisatellites is calculated (thus excluding homopolymers), 68% of HCMV genes still contain TRs, while this proportion is decreased to 39% for HSV-1.

• Open in new tab
• Download powerpoint

FIG 3

Tandem repeats (TRs) in the HCMV genome. TRs in reference strain Merlin were identified, and orthologous repeats were searched for in a data set of 124 complete HCMV genomes. Only TRs with orthologs in >50% of strains were included in the analysis. (A) TR nucleotide content (9,008 nt) in coding and noncoding (intergenic, intron, and ncRNA) regions compared to the distribution of total nucleotides (231,784 nt) over these regions. (B) Similar to panel A, where TR and total nucleotide distributions over unique long (UL), unique short (US), and internal repeat long and short (IRL-IRS) genome regions are compared. Percentages in panels A and B do not add up to 100% because of rounding errors. (C) Conservation of TRs between coding and noncoding TRs and between different TR types (homopolymers, microsatellites, and minisatellites). TRs were reported to be conserved if >50% of strains had identical repeat sequences and copy numbers.

TR polymorphisms in noncoding regions might have a profound impact on gene regulation and expression by altering binding sites for regulatory proteins, chromatin structure, transcript stability and transcription, splicing, or translation efficiency (87–90). However, we can only speculate about their effects based on sequence data alone. Overall, selection seems to have constrained the presence of unstable TR elements inside coding regions, illustrated by the discrepancy of TR frequencies between coding and noncoding regions (Fig. 3A) and the higher level of conservation of TRs in coding regions (Fig. 3C). Therefore, it is conceivable that some TR variations in coding regions might have specific functions or provide the virus with greater adaptability because of their intrinsic instability (so-called “evolutionary tuning knobs” [91]). Mutation rates in TRs can be up to 100,000 times higher than those in other parts of the genome. Therefore, we assessed the potential impact of TR variation in coding regions on the encoded proteins (see Table S5 in the supplemental material). Nine out of the 13 variable TRs inside protein-coding regions constitute variability in the longer and nonperfect TRs, as determined by TRF analysis. Diversity in these minisatellites comprises mostly variations in repeat sequence and period length, caused by nucleotide divergence in these areas. The UL50 and UL111A genes contain variable trinucleotide microsatellites, causing amino acid stretches of various lengths (Table 4). Homopolymer length variation is present in the RL12 and UL1 genes, but it does not cause frameshifts and results only in amino acid divergence. In fact, none of these 13 variable TRs in protein-coding regions cause frameshifts that disrupt ORF integrity. Because TRs cataloged as conserved (>50% of strains with conserved repeat sequence and copy number) could also contain variation in a minority of strains, we analyzed an additional set of 53 protein-encoding TRs with variations in period length or copy number below the 50% threshold (see Table S5 in the supplemental material). These TRs comprised 26 homopolymers, 5 microsatellites, and 22 minisatellites. While most of these TR variations either were conserved or led to amino acid variations and indels without a clear repetitive character at the protein level, seven TRs caused the occurrence of repetitive single-amino-acid stretches of various lengths (Table 4). These stretches contain mostly small and hydrophilic amino acids, suggesting selective constraints toward these residues in coding TRs (92). It is assumed that TRs inside coding regions form flexible, unstructured, and hydrophilic loops (75). These amino acid loops might be involved in protein-protein interactions that could be altered by changes in loop length. For example, the variable proline stretch in UL50 might have functional consequences for the efficiency of nuclear egress of HCMV capsids. Together with pUL53, pUL50 forms the nuclear egress complex (NEC). A random screen for dominant negative mutants of M50 (the MCMV homolog of UL50) identified this proline-rich motif to be essential for nuclear egress, and this finding was confirmed for HCMV UL50 (93). The authors of that study suggested that this motif likely controlled a binding site for a NEC interaction partner. The latter could be HCMV pUL97, a protein kinase that was recently found to phosphorylate S216, a site neighboring the proline stretch. This phosphorylation modulates NEC localization and nuclear egress (94, 95). It might be of interest to assess the effect of the large length heterogeneity (4 to 12 residues) in the proline motif on nuclear egress efficiency. Finally, homopolymer length variation in the UL111A and UL133 genes caused frameshifts that resulted in premature ORF termination in strains BE/16/2010 and BE/17/2010 (UL111A) and in strain BE/2/2012 (UL133) (see Tables S2 and S5 in the supplemental material).

Wild-type HCMV strains contain ORF-disrupting mutations in a wide range of nonessential genes.The accumulation of gene-disrupting mutations in cell culture-passaged HCMV strains is a well-described phenomenon (35, 96). As first suggested for isolates of koi herpesvirus, some disrupting mutations might also occur in vivo (97). Recently, a few studies have indicated that some HCMV mutants may indeed be present in clinical isolates prior to culture passage (36, 62, 98). In particular, strains JP (62) and BE/21/2010 (36) were sequenced directly from clinical material and displayed disruptive mutations in the RL5A and UL111A genes and in the RL5A, UL9, and UL150 genes, respectively. Furthermore, we showed that RL5A, UL1, UL9, and UL111A mutants identified in three additional strains after limited passaging were also present in the original clinical isolate (36). With our current data set providing full-genome information for 96 clinical isolates, supplemented with 32 previously reported sequences, we have an ideal opportunity to further assess the occurrence of ORF-disrupting mutations in clinical HCMV isolates and provide a more detailed estimate of the mutation frequency in different genes.

An overview of all genes that contain ORF-disrupting mutations in one or more isolates is presented in Table 5, and more detailed information is provided in Table S2 in the supplemental material. In total, 26 of 170 genes (15%) have a disrupted ORF in at least one clinical isolate. Unsurprisingly, none of these genes are essential for growth on fibroblast cells, although UL30 was found to be growth augmenting (99, 100). Looking at the distribution of mutants over different clinical isolates (Fig. 4), only 28 of 124 isolates (23%) have the complete set of 170 intact genes, with the other isolates having 1 (33%), 2 (27%), 3 (13%), or 4 (3%) mutated genes. With these data, we show that only 1 out of 4 clinical isolates is genetically intact and that gene-disrupting mutations are extremely common in a defined set of nonessential genes. We cannot rule out the possibility that even more genes are mutated in our isolates, as discussed in Text S1 in the supplemental material. Some mutations affect only a small proportion of the ORF (e.g., in the US9 and US27 genes) (Table 5). These cases could be better described as variants than as mutants if the encoded gene products are not affected by these N- or C-terminal deletions, but functional experiments are needed to evaluate this.

• Open in new tab
• Download powerpoint

FIG 4

Distribution of ORF-disrupting mutations. Shown is a graphical representation of the distribution of ORF-disrupting mutations over 124 clinical isolates. Rows represent different isolates, and columns represent all 26 genes containing disruptive mutations. Disrupted genes are represented in red, and intact genes are represented in green.

As our strains were minimally passaged in fibroblast cell culture, some mutations might be artifacts of culture adaptation. Therefore, the presence of mutations in the original clinical material was assessed by direct PCR sequencing of these samples, if available. All mutations that were confirmed with this procedure are shown in Table 5 and in Table S2 in the supplemental material. Importantly, all but one of the mutations that were characterized in the original clinical material could be confirmed to be unrelated to culture amplification. Mutations in the RL1, RL5A, RL6, UL1, UL9, UL30, UL40, UL111A, UL142, UL150, US7, US9, and US27 genes are thus indeed present in clinical isolates. Moreover, the occurrence of identical mutations in unrelated and geographically distinct isolates confirms their circulation in natural populations. The only mutation identified as an artifact of culture passage was a substitution in the first splice donor site of the UL128 gene in strain BE/11/2011; the specific case of mutations in the RL13 and UL128 genes is discussed in further detail in Text S1 in the supplemental material.

Mutated genes are highlighted on the genome map in Fig. 1 with a color code depicting mutation frequency. The RL11 gene family stands out, with 7 out of 14 members containing disruptive mutations (or 6 if RL13 is omitted). Furthermore, 4 out of 5 genes that are mutated in >10% of isolates are part of this family. RL11 genes share homology with the CR1 domain of the adenovirus E3 genes through their RL11 domain (101). The encoded proteins have similarities to the IgD family, and because of their hypervariability (see “HCMV displays the highest level of genetic diversity of all human herpesviruses,” above), it is believed they could have a function in modulating variable host proteins. Immunomodulatory capacities have recently been proposed for UL7 (102), UL11 (103), and RL12 and RL13 (73). While RL13 mutants are probably culture artifacts and RL12 and UL11 mutants are rare and unconfirmed in clinical material, mutations in RL5A, RL6, UL1, and UL9 are much more common and definitely circulate in the host. Functional data for these genes are limited, with UL1 being implicated as a tropism factor (104) and UL9 being implicated in growth temperance (100). While our analysis shows that almost 35% of strains have a mutation in UL9, a previous study that characterized UL9 sequences in unpassaged and moderately and highly passaged isolates from a diverse geographical background found only 2 mutations in 41 strains (5%) (GenBank accession numbers DQ847465 to DQ847505). It is unclear whether this discrepancy is due to geographical differences in UL9 mutations, the types of patients involved, or the body compartment that was sampled.

Several genes in the RL1 family (2/2), the US6 family (3/6), and the US12 family (2/10) are also affected. For the RL1 family genes RL1 and UL145, functional knowledge is lacking, and in both cases, only one strain contains a mutation. Members of the US6 family have established or tentative immunomodulatory functions, interfering with the major histocompatibility complex class I (MHC-I) antigen-processing pathway (71, 105). The role of US6 in inhibition of antigen peptide transport to MHC-I molecules has been well described (106). Studies on US7 and US9 function are much scarcer (107–109), but immunomodulatory functions have been predicted (71). The natural occurrence of mutants in US6, US7, and US9 may suggest some functional redundancy, as the US2, US3, US10, and US11 genes all target MHC-I antigen presentation (110). The US12 gene family encodes 10 seven-transmembrane proteins, with some members recently being associated with immunomodulation (111, 112) and cell tropism for epithelial and endothelial cells (113). Functions of the US12 and US13 genes have not yet been revealed. These different HCMV gene families exist as clusters of adjacent (apart from the RL1 family) and (distantly) related genes that probably originated from duplication of an ancestral gene (114). As recent experiments with poxviruses have shown, these so-called genomic accordions can rapidly expand under strong selective pressure and provide dsDNA viruses with low mutation rates an alternative mode to evolve more quickly under specific circumstances (115, 116). Subsequently, adaptive mutations in these new copies can accumulate at a much higher rate than in one isolated gene. Finally, when the selective environment changes, genes that did not acquire beneficial functions can be removed, and the accordion contracts. Especially for the RL11 and US6 gene families, where mutant genes were confirmed to be present in clinical isolates, this accordion contraction could be currently ongoing.

Several additional, isolated genes are mutated in the original clinical material. The fact that almost 10% of strains are affected in the viral interleukin-10 (IL-10)-encoding UL111A gene is striking. UL111A encodes separate transcripts during productive (cmvIL-10) and latent (LAcmvIL-10) infection, and both gene products are affected in all 12 mutant isolates (Fig. 5). cmvIL-10 binds and signals through the human IL-10 receptor and mimics its immunomodulatory properties. It has been shown to inhibit the production of proinflammatory cytokines and MHC-I and -II expression in monocytes and to stimulate monocyte differentiation to a phagocytic phenotype and B cell proliferation and differentiation. It is believed that these concerted cmvIL-10 actions have an important impact on the immune system's capacity to control HCMV replication (117). In this regard, cmvIL-10 has been shown to impair cytotrophoblast remodeling of the uterine vasculature, thereby possibly enhancing congenital disease (118). The latency-associated LAcmvIL-10 product cannot signal through human IL-10 receptors in the same fashion but was recently shown to upregulate the expression of cellular IL-10 and CCL8 (119). Another mutated gene with potential implications for the immunomodulatory capacities of isolates is UL40. The UL40 signal peptide is necessary for the cell surface expression of HCMV pUL18 and HLA-E molecules, both of which are natural killer (NK) cell ligands that can inhibit NK cell activation in the absence of normal MHC-I antigen presentation (120). In 5.5% of isolates, UL40 has a substitution in its original start codon, presumably leading to translation initiation from an alternative start codon located 15 nt downstream and truncating the signal peptide (see Table S2 in the supplemental material). Whereas this mutation was previously found only in strain 3157, we identified six more instances of signal peptide truncation in our isolates. By using UL40 from strain 3157 in comparison to wild-type UL40, it was recently demonstrated that this mutation did not affect UL40 translation or pUL18 surface expression but did inhibit the surface expression of HLA-E and thereby sensitized infected cells to NK cell lysis (120). For both the UL111A and UL40 genes, unrelated strains display identical deletions or substitutions, suggesting wide circulation of mutants in the population (Fig. 5; see also Table S2 in the supplemental material). Studying the impact of UL111A and UL40 mutants on the immune-evasive potential of strains and their implications for strain pathogenicity might be worthwhile, as initial findings for UL40 illustrate (120). The remaining genes containing disruptive mutations are discussed further in Text S1 in the supplemental material; mutants of these genes in clinical material were rare and/or unconfirmed.

• Open in new tab
• Download powerpoint

FIG 5

ORF-disrupting mutations in the UL111A gene. Shown is a nucleotide alignment of wild-type UL111A (strain Merlin) and all 12 mutants. Countries of isolation are listed for all strains with the international two-letter code (GB, Great Britain; BE, Belgium; CZ, Czech Republic; DE, Germany). Mutations (deletions, insertions, and substitutions) are highlighted in red, and the predicted stop codons are underlined, with untranslated sequences after stop codons being crossed out. Introns have a gray background, unless they are aberrantly translated because of the deletion of splice donor sites. LAcmvIL-10 transcripts are similar, but the second intron is not spliced, with translation proceeding into it.

HCMV evolution has been shaped by pervasive recombination.An important role for recombination in HCMV evolution has long been suggested. Given the common occurrence of multiple infections in a single host, the potential for recombination is obvious (reviewed in reference 15). Recombination has been identified in individual gene sequences and by the small amount of gene linkage observed (98, 121–124). However, overall recombination on a full-genome scale has not yet been quantified.

Since standard phylogenetic trees cannot account for recombination and are incorrect if recombined sequences are included, we have analyzed phylogenetic relationships by constructing a split-decomposition network (Fig. 6A) (51). The large numbers of reticulate connections that are apparent at the root of the network indicate conflicting evolutionary signals such as those caused by recombination. By using the Phi-test for recombination, it was confirmed that recombination between isolates was detected (P = 0.0) (54). To further visualize recombination events, nine strains were selected from the network, and a BootScan analysis was performed based on a sequence alignment of these strains only (Fig. 6B). The constant shifting of phylogenetic relationships along the genome provides further evidence of numerous recombination events.

• Open in new tab
• Download powerpoint

FIG 6

Widespread recombination between HCMV strains. Recombination between separate HCMV strains was analyzed. (A) Neighbor-net split network of 124 full-genome sequences showing numerous reticulate connections that are indicative of recombination. The Phi-test for recombination gave strong statistical evidence for recombination. Countries of isolation of different strains are represented with the international two-letter country codes at the beginning of strain names (CN, China; KR, South Korea; IT, Italy; US, United States; other codes are defined in the legend to Fig. 5). Asian strains JHC and HAN are highlighted with a gray background. (B) BootScan analysis of 9 strains highlighted in the split network. Strain BE/25/2010 (highlighted in black) was used as a reference strain.

When split networks were constructed for HSV-1, VZV, and EBV, similar evidence for recombination was detected (P values are shown in Fig. 7). Previous studies have reported distinct geographic clusters for HSV-1, VZV, and EBV strains (48, 125–127). Split networks for these viruses indeed show distinct clusters (Fig. 7), while such clusters do not become apparent from our HCMV split network. However, apart from two Asian strains, all isolates were collected from European and North American patients. Interestingly, Asian strains JHC (128) and HAN are neighbors in the split network (highlighted in gray in Fig. 6A), but they do not cluster separately from European/North American strains. Additional full-genome sequences from Asian and African strains will be necessary to investigate the potential existence of separate geographic clusters of HCMV.

• Open in new tab
• Download powerpoint

FIG 7

Diversity and recombination in HSV-1, VZV, and EBV. Using the full-genome sequences listed in Table S1 in the supplemental material, neighbor-net split networks were constructed for HSV-1, VZV, and EBV strains. In all three cases, statistically significant evidence for recombination was detected. Countries of isolation of different strains are represented with the international two-letter country codes at the beginning of strain names (KE, Kenya; JP, Japan; HK, Hong Kong; NG, Nigeria; GH, Ghana; CA, Canada; RU, Russia; NL, the Netherlands; MA, Morocco; MX, Mexico; other codes are defined in the legends to Fig. 5 and Fig. 6). In all three networks, distinct clusters are recognizable. Strains chosen for recombination analysis with RDP3 are underlined.

The amount of recombination was quantified by detecting individual recombination breakpoints with RDP3 (53). To allow comparison of recombination densities in HCMV, HSV-1, VZV, and EBV genomes, we selected a group of nine full-genome sequences for each virus species, i.e., the total number of genomes that was available for EBV. Strains were selected to properly reflect total diversity (Fig. 6A and 7). Of these four human herpesviruses, HCMV clearly has the highest recombination density (Table 6). There is a statistically significant, positive correlation between overall nucleotide diversity (Table 1) and recombination density (Table 6) for these four viruses (P = 0.017 and ρ = 0.98 by Pearson's product-moment correlation). In comparisons of clusterings of isolates in the split networks (Fig. 6A and 7), there is a tentative inverse relation between the existence of distinct clusters and the recombination breakpoint density. VZV has the lowest recombination density, and the divergence of clusters is very clear, with few recombinants between clusters. HSV-1 and EBV have intermediate recombination densities: distinct clusters are still recognizable, but their delineation is less pronounced because of the presence of intermediate recombinants. Finally, HCMV has a markedly higher recombination density, resulting in a star-like phylogeny with no apparent clustering. Hypothesizing that these clusters or “genome types” have evolved during human radiation across the planet, we may now see a fading of geographical genome types because of increased global travel and mixing of populations and viruses. This hypothesis is supported by findings that the geographic separation of VZV clades might be slowly fading because of recent immigration (129). Given their different recombination densities, distinct herpesviruses would currently be at different stages of this process. In comparisons of recombination in HSV-1 and VZV, it has been proposed that inherent disparities in biological characteristics might be at the root of recombination potential (130). Higher recombination rates in HSV-1 were attributed to more frequent reactivation, longer episodes of asymptomatic shedding, and subsequent increased occurrence of multiple infections. In the same fashion, its superior immunomodulatory capabilities and broad cell tropism might explain the even higher capacity for recombination in HCMV.

While it is clear that recombination at the genome level permits the exchange of alleles at separate loci, we also wanted to evaluate the contribution of recombination to the generation of variation at the level of individual genes. Evidence for recombination in codon-aligned gene sequences was assessed by three different methods. The Phi-test for recombination simply determines whether there is statistically significant evidence for recombination, while GARD and RDP3 identify specific recombination breakpoints. An overview of results for all genes is presented in Table S6 in the supplemental material. Recombination was confirmed in 93 out of 170 genes (55%) by at least 2 out of 3 methods. The densities of breakpoints varied enormously between different genes, as illustrated in Fig. 1. Gene-specific recombination densities were grouped in different gene families, gene conservation groups, and gene functions, as described above for diversity (see second panel in Fig. 2 for gene families). There was a statistically significant difference between recombination densities in different gene families (P = 0.0031 by KWt) but not in different conservation categories (P = 0.13 by KWt). For the former, none of the pairwise comparisons were significant. Overall, gene families with higher levels of diversity display higher recombination densities, although there are clearly some exceptions in conserved genes, such as IRS1 and TRS1 of the US22 family (Fig. 2). When divided over different gene functions, there was still a statistically significant difference between recombination densities (P = 0.023 by KWt; no significant pairwise comparisons). The discrepancy between overall KWt results and the results of pairwise comparisons for recombination over different gene families and functions is caused by the rigorous correction for multiple testing in the latter. The numbers of genes assigned to gene families and functions vary widely, reducing the power of the tests. In summary, our data demonstrate the potential for recombination to generate diversity both at the level of the individual gene as well as over the complete genome. The high recombination rate in HCMV compared to that in other human herpesviruses could be one of the most important reasons for its higher level of strain diversity. Reciprocally, the increased diversity could also result in a higher sensitivity for detection of recombination events that are undetectable in more conserved herpesviruses.

Positively selected residues provide a genetic fingerprint of the evolutionary arms race between virus and host.While intermolecular recombination is a very strong mechanism to increase nucleotide diversity, its ability to do so relies on preexisting variation through mutations caused by polymerase infidelity. Obviously, the accumulation of these mutations is dependent on the selective pressure that acts upon them. The persistent nature of HCMV and its constant interplay with host immune components highlight the potential selective pressure that might act on gene variants. To characterize the overall selective pressure acting on viral genes, we calculated the ratio of nonsynonymous substitutions per nonsynonymous site (dN) to synonymous substitutions per synonymous site (dS) for each gene alignment (Fig. 8; see also Table S4 in the supplemental material). Most HCMV genes seem to be under strong evolutionary constraints, as 96% of genes have a dN/dS ratio of <0.5 and 64% of genes have a ratio of <0.2. Recent studies of HSV-1 and MCMV also found a predominance of negative selection, but the proportion of genes with a dN/dS ratio of <0.1 is much higher for HCMV (38%, versus 7% for HSV-1 and 13% for MCMV) (48, 63). A few genes have dN/dS ratios closer to 1, indicating neutral selection or genetic drift. Most of these genes are poorly characterized, but the tegument protein ppUL25 was found to be a major target of anti-CMV antibodies (131). This strong negative/purifying selection is an indication of the excellent adaptation of HCMV to its human host, with most mutations having negative fitness effects and being quickly removed.

• Open in new tab
• Download powerpoint

FIG 8

The majority of HCMV genes are under strong purifying selection. The selection mode acting on genes is represented by calculation of dN/dS ratios. A ratio close to zero indicates strong negative/purifying selection, and a ratio close to 1 indicates neutral selection or genetic drift. A ratio significantly higher than 1 indicates positive/diversifying selection. Genes are binned in groups with similar dN/dS ratios in steps of 0.1.

Notwithstanding the strong purifying selection acting on most HCMV genes, individual residues might be experiencing a different selective regimen. To identify positive/diversifying and negative/purifying selection at the codon level, we made use of four separate algorithms (SLAC, FEL, MEME, and FUBAR) included in the HyPhy package for phylogenetic hypothesis testing (56) (see Table S7 in the supplemental material). Since these methods have different sensitivities and specificities, we retained only sites that were independently confirmed by at least two out of four methods. Of 68,287 codons analyzed, 431 (0.6%) showed evidence of positive selection (in 105 genes), and 7,731 (11.3%) showed evidence of negative selection (in 169 genes), again demonstrating the predominance of negative selection acting on HCMV genes. However, there were clear discrepancies in the distribution of positively selected residues, as shown in Fig. 1. There was a statistically significant difference between gene families (P = 0.0054 by KWt) and gene conservation groups (P = 0.00055 by KWt); genes with higher levels of variability generally display higher percentages of positively selected codons (statistically significant only for cytomegalovirus-specific versus core genes [P = 0.0033 by WRSt]). Clearly, there are multiple exceptions, most notably the MHC family (UL18 and UL142) and the RL1 (RL1 family), US7 (US6 family), IRS1 and TRS1 (US22 family), and US14 and US18 (US12 family) genes, which are subjected to higher levels of positive selection than would be expected from their diversity (Fig. 2, third panel). Their products function in the evasion of both adaptive immune responses (UL18, UL142, and US18) and innate antiviral mechanisms (IRS1 and TRS1), or their functions are not yet characterized. Hence, it could be suggested that these functionally uncharacterized genes might also interact with host antiviral mechanisms. Also, for negative selection, there was a significant difference between gene families (P = 0.00023 by KWt) and gene conservation groups (P = 0.015 by KWt), but this was unrelated to gene diversity (Fig. 2, bottom). No statistically significant differences in selection between gene functions were observed (P = 0.16 for positive selection; P = 0.32 for negative selection [determined by KWt]).

Of the genes experiencing the highest percentages of positive selection, there is a clear predominance of genes modulating host immune and antiviral pathways (Fig. 1; see also Table S7 in the supplemental material). The UL147 gene encodes an α-chemokine (132), UL142 interferes with NK cell activation (133), IRS1 inhibits the protein kinase R antiviral pathway (134), and UL20 is poorly characterized but has also been implicated in immune evasion (135). The role of the UL40 signal peptide in inhibiting NK cells is discussed above (see “Wild-type HCMV strains contain ORF-disrupting mutations in a wide range of nonessential genes”). Interestingly, 6 out of 8 residues under positive selection in pUL40 are located within this signal peptide, and more specifically, 4 out of 8 are located within the nonamer peptide that is processed to bind the viral MHC-I mimic pUL18 and HLA-E. For pUL18, a crystal structure of its interaction with human β-2-microglobulin (MHC-I light chain), an actin peptide, and the inhibitory LIR-1 receptor was reported (136). We have visualized pUL18 residues under positive selection on the pUL18 structure (Fig. 9). All positively selected residues are located on the surface of pUL18, with L265 interacting with β-2-microglobulin and D248, V250, and E252 being situated adjacent to LIR-1 binding sites. Although the overall evolutionary mode of HCMV is strongly shaped by constraints and purifying selection, these examples illustrate that specific protein residues can experience diversifying selective pressure through their interactions with the host.

• Open in new tab
• Download powerpoint

FIG 9

Residues under diversifying selection in pUL18. Codons under positive/diversifying selection in the UL18 gene were determined with the SLAC, FEL, MEME, and FUBAR algorithms of the HyPhy package. Sites that showed significant evidence of positive selection by at least two of four methods are represented in red on the protein structure of pUL18 (green). The structure shows a complex of pUL18 (a viral MHC-I homolog), human β-2-microglobulin (blue) (a MHC-I light chain), and an actin peptide (pink) bound to the inhibitory immunoglobulin receptor LIR-1 (yellow) (136). The three-dimensional structure is visualized from two opposite angles. All selected residues are located at the surface of pUL18.

Concluding remarks.In this study, we have applied our previously described workflow (36) to the characterization of HCMV genomes of 96 distinct clinical isolates from Belgian and Czech patients. Based on a comparative analysis of 124 full-genome sequences, HCMV interstrain nucleotide diversity, π, was estimated to be 0.021, significantly higher than the diversity in other human herpesviruses. Nevertheless, overall levels of purifying selection were very high, reflecting the remarkable adaptation of HCMV to its human host. Given the proofreading capacity of the HCMV DNA polymerase, additional strategies are necessary to generate diversity apart from replication error. Because of the wide range of immune-evasive functions, multiple infections are common, permitting the virus to recombine extensively. Furthermore, gene duplications can effectively enlarge mutational space. Finally, interactions with host components generate positive selective pressure at specific loci, which may help the virus avoid immune recognition.

We have demonstrated the widespread occurrence of gene-disrupting mutations in wild-type HCMV strains, unrelated to culture passage. While some of these mutants might be evidence of ongoing contraction of gene family duplications, others could have implications for the immunomodulatory and, ultimately, the pathogenic potential of the isolates involved. In particular, a closer look at mutants of the UL40 and UL111A genes is warranted. We are currently retrieving clinical data regarding our patient population to explore associations of gene variants and mutants with disease outcome.

Our data set is completely derived from European patients. Currently, only two complete genome sequences from patients outside Europe and North America are publicly available. The availability of more complete genomes, especially from African and Asian patients, will reveal whether there are geographical discrepancies in gene diversity and mutational patterns. Additionally, the majority of our isolates (74/96) were derived from urine samples. Our preliminary results show identical mutations in different body compartments for several genes (including nasopharyngeal isolates, blood, amniotic fluid, and bronchoalveolar lavage fluid), but we cannot exclude the possibility that specific mutations can influence the tissue tropism of the viral strain. In the future, it will be interesting to see whether there are differences in the occurrence of disruptive mutations depending on the body compartment that is sampled.

To our knowledge, this study is the most comprehensive analysis of genetic variability and evolution in HCMV to date, providing both conceptual insights into diversity generation and a large source of sequence information of outstanding value for functional experiments with this important human pathogen.