Mapping the Landscape of Host-Pathogen Coevolution: HLA Class I Binding and Its Relationship with Evolutionary Conservation in Human and Viral Proteins

Conservation scores for an analyzed protein were computed using the ConSeq server (http://conseq.tau.ac.il/), which estimates the conservation score C(i) for each protein site i using a phylogenetic tree built from a set of homologous sequences. The tree is used to infer the evolutionary rates (log probabilities of substitution) for each site along the given protein (2) (Fig. 1A). Conservation scores were computed only for proteins which had at least five homologous proteins in the UniRef100 database (release 11.0). For computational efficiency, no more than 100 homologues were used. Homologous proteins were aligned using the MUSCLE program, version 3.6. The site conservation scores are computed in the context of the entire protein as ConSeq uses aligned homologues to compute the phylogeny in which the targeted protein resides. This analysis also incorporated adjustment for variables that may affect P values, including various numbers of protein homologues in conservation score computations, as well as various protein length distributions. Comparisons of the effects of using inter- versus intraspecies homologous proteins for estimating the evolutionary rates are provided for HIV in Table S6 in the supplemental material.

The binding scores are based on experimental measurements characterizing individual HLA-peptide interactions, as catalogued in the Immune Epitope Database (IEDB) (48, 50), as well as known HLA-peptide binding configurations (32). Binding energies of HLA-peptide complexes were systematically estimated using the adaptive double-threading (ADT) structure-based approach (32) for estimating the binding energy of a major histocompatibility complex (MHC)-peptide complex. The method estimates the 50% inhibitory concentrations ([IC₅₀s] a measure of the binding affinity) after threading both target peptides and HLA proteins (in particular, the known contact residues shown in red on the HLA structure in Fig. 1A) onto solved HLA-peptide complex structures. Here, we have focused on 9-mer peptide targets as the vast majority of known HLA class I epitopes are of this length.

The model parameters were fit to log IC₅₀s obtained from the IEDB for ∼34,000 experiments covering 35 HLA-A and HLA-B alleles. We excluded all HIV epitopes from these data for training in order to avoid a possible bias in the analyses of HIV viral load data. The ADT model can provide estimates for HLA molecules other than the limited number on which it was trained by threading the arbitrary HLA sequence onto the structure of another similar HLA protein and using the estimated model parameters, which generalize for the entire HLA allelic family. We analyzed 95 HLA-A and HLA-B alleles from a Caucasian population in Australia (43), a cohort which provided a total coverage of 95% of HLA-A alleles and 90% of HLA-B alleles in Europe. We note that of these 95 alleles, empirical binding data were provided for only 35 alleles. Due to lack of sufficient experimental measurements for HLA-C alleles in the IEDB, we did not consider these alleles in the current study. The binding score at a given position along the protein is the sum of binding energies for nine overlapping peptides, which measures the probability that the site will be visible to immune surveillance. The binding energy model provides an energy estimate, E_a(e) = E_a(e₁,e₂,…,e₉), where a denotes the index of an HLA allele, and e is (e₁,e₂,…,e₉) the 9-mer peptide. The model is fit to the logarithm of the IC₅₀ measurements for different allele-peptide combinations. Therefore, the probability of peptide presentation is proportional to e^−Ea(e), and high energy indicates low presentation probability and vice versa. In order to estimate the log probability of presentation of a single site in a protein, presentation probabilities of all peptides straddling that site need to be considered. As an estimate that is robust to prediction errors, we define the binding score B(i) for the i-th amino acid in the sequence s = (s₁,s₂,…,s_N) of an arbitrary N-long protein whose segments may be presented by an HLA molecule as follows: \( \(B_{a}(i)\ {=}\ {-}\ {{\sum}_{j{=}i{-}8}^{i}}E_{a}(s_{j},\ s_{j{+}1},\ {\ldots},\ s_{j{+}8})\) \). The predicted binding energies are highly correlated with the true experimentally measured binding energies (Spearman correlation of >0.75), and for some alleles the accuracy of prediction for our method and for other prediction methods (see, for example, references 3, 7, 16, 31, 44, 45, and 49) is believed to be within the accuracy of the IC₅₀ measurement error. A recent analysis of various prediction methods can be found in Nielsen et al. (46).

The allele efficiency score r is defined as the Spearman correlation coefficient of the binding score and the conservation score for a given protein. A positive score indicates preferential targeting of conserved regions, and a negative score indicates preferential targeting of variable regions.

Statistical comparisons were corrected for the potential lack of independence of measurements resulting from the hierarchical structure of the groups, such as correlations generated by similarities within viral families and within HLA supertype classes. The sampling procedure for these analyses, along with a detailed analysis of potential sources of bias in our analysis and a detailed discussion of the statistical power of the methods used, is available in the supplemental material.

We obtained binding data for training the HLA binding predictor from the IEDB resource, consisting of experimentally measured binding affinities (IC₅₀s) for ∼34,000 HLA-peptide pairs, spanning 35 HLA-A and HLA-B alleles. For the viral targeting analysis, we collected the sequences of the most commonly studied human and plant viruses from the NCBI virus database. A full list is provided in Table S2 in the supplemental material. NCBI GenBank accession number information for individual viral proteins is available in Table S8 in the supplemental material. We based our analysis of the human proteins with disease-associated mutations on a version of the Online Mendelian Inheritance in Man (OMIM) database (McKusick-Nathans Institute of Genetic Medicine, Johns Hopkins University, Baltimore, MD, and NCBI, National Library of Medicine, Bethesda, MD, as of February 2007 [http://www.ncbi.nlm.nih.gov/omim]). A control group of randomly sampled proteins was obtained from the UniRef50 database, release 11.0 (http://www.pir.uniprot.org/search/textSearch_NR5.shtml), which contains representatives of natural protein clusters with >50% interclass similarity. The random sample did not include any human proteins. BLAST queries were performed against the UniRef100 database.

A total of 200 Western Australian HIV Cohort participants provided samples with consent for HLA typing. HLA class I genotyping was resolved to the four-digit level based on the sequence of exon 2 and 3 (exon 2/3) and using standard sequence-based typing (SBT). Allele and heterozygote ambiguities between alleles not identical in exon 2/3 were resolved using alternative primers.

In order to test the hypothesis that there is a correlation between conservation and binding patterns in proteins known to interact with HLA molecules, we examined the relationship between the efficiency scores of 95 HLA-A and HLA-B alleles (which account for more than 90% of all HLA alleles in European Caucasians) for target proteins spanning 4,761 human proteins derived from the OMIM database and 52 human viruses (see Table S2A in the supplemental material). We also examined these same correlations in 70 plant viruses (see Table S2B), as well as a control sample of 3,800 random proteins obtained from the UniRef50 database expecting that we would observe significant, but weaker, efficiency scores because of the sharing of conservation patterns between the proteins that typically do not interact with the immune system and the human proteins and human viruses that do. These control groups were also used to ascertain the significance of targeting efficiencies and variations in targeting efficiencies on human proteins and human viruses.

In order to verify that our conclusions were not dependent on the use of a specific predictor, we also used the NetMHCpan online prediction method (3) to replicate our results on Gag efficiency scores and HIV viral load correlations. Negative correlations with viral load were also obtained using this method (data not shown). In addition, we also conducted an analysis of HIV efficiency scores based on experimentally determined T-cell epitope maps and compared this analysis to the one obtained using predicted epitope maps. We found more evidence for HLA targeting of conserved regions when we used the measured epitope maps than when using the predicted ones, indicating that our analysis may underestimate the extent of this phenomenon (see supplemental material and Table S7).

We introduce a novel score termed “targeting efficiency” to quantify the relationship between HLA binding and conservation of target peptides. The HLA allele targeting efficiency score is defined to be the Spearman rank correlation coefficient between binding scores and conservation scores for amino acids along a given protein. Positive scores denote a preference for binding conserved regions while negative scores indicated a preference to bind to variable regions. The process of calculating a targeting efficiency score for the example pair of the herpesvirus-1 capsid triplex subunit 1 protein and the HLA-A*2402 molecule is illustrated in Fig. 1B and C. In addition to calculating targeting efficiency scores for individual proteins, the approach also allowed an efficiency score to be calculated for an entire proteome by concatenating all protein sequences of the given pathogen for which data were available.

In order to examine whether HLA molecules preferentially target conserved targets on human proteins, we analyzed a set of 4,761 proteins spanning the entire OMIM database of human-disease associated proteins (http://www.ncbi.nlm.nih.gov/omim/). This analysis showed that both HLA-A and HLA-B molecules tend to preferentially bind to conserved areas of human proteins (Fig. 2), as first proposed by Hughes and Hughes in 1995 in an analysis of 34 oncoproteins (23). These results are also in agreement with those reported by Yeager et al. (64). The average targeting efficiencies were significantly higher for both HLA-A (P < 10⁻⁷⁶) and HLA-B alleles (P < 10⁻⁸) for OMIM proteins than for a random sample of 3,800 proteins (UniRef50 database) (Fig. 3). Similar targeting preferences were obtained from a smaller comparison using 300 randomly selected non-disease-related human proteins and 300 randomly selected UniRef proteins (see Table S3 in the supplemental material), verifying that this property is not specific to disease-associated human proteins. Importantly, the reported differences in efficiency distributions are not a consequence of the differences in the raw scores of the two parameters used to compute them as neither the binding scores (P = 0.96) nor the conservation scores (P = 0.99) were distributed differently in human and UniRef random protein sets. Additionally, significantly higher binding scores were observed at disease-associated mutation sites (P < 10⁻¹⁰ over 2,825,558 sites), providing further evidence of HLA coevolution with its targets.

We were also interested in exploring the relationship between efficiency scores and the assignment of HLA alleles to loci and to supertype groups (described in reference 29; see also Table S1 in the supplemental material). Therefore, we examined the efficiency scores of HLA-A versus HLA-B alleles for both the human and random protein sets described above. We found a marked preference for HLA-A to target conserved regions of human proteins (HLA-A ≫ HLA-B, P < 10⁻³⁰⁰) compared with the random UniRef50 sample (HLA-B > HLA-A, P = 0.00007). However, further analysis indicated that while the distribution of efficiency scores appears more uniformly positive for the HLA-A allele families (Fig. 4), a limited number of HLA-B alleles also preferentially bind conserved human protein sequences. For example, the B58 supertype (e.g., HLA-B57 and -B58 alleles) and a subgroup of HLA-B7 supertype alleles (HLA-B55 and -B56) have higher efficiency scores than most other B alleles for human proteins (Fig. 4). Interestingly, all of these alleles form a distinct HLA cluster closely associated with chimpanzee Patr-B alleles (see Fig. S2 in the supplemental material). Conversely, HLA-B alleles from the B44, B27, and B7 supertype families (other than B55 and B56) have low or negative average allele efficiencies for human proteins (Fig. 4). These results indicate that grouping HLA alleles by their abilities to target conserved regions of human proteins does not strictly follow their supertype classification, but it does reflect their phylogenetic history as alleles that are more likely to target conserved areas of human proteins cluster with chimpanzee Patr-B alleles (41).

To determine if binding preferences for viruses follow a trend similar to the one observed for human proteins, we computed the HLA efficiency scores for 52 human viruses and 70 plant viruses (listed in Tables S2A and S2B in the supplemental material). Figure 5 (see also Fig. S3) shows the distribution of efficiency scores as heat maps in which the viruses (x axis) are grouped by their families, and the HLA allelic variants (y axis) are organized into HLA supertype families characterized by their HLA peptide binding preferences (56). Figure 1D also presents the distribution of efficiency scores for HLA-A and HLA-B alleles for a selection of human viruses. Importantly, neither the binding scores (P = 0.87) nor the conservation scores (P = 0.59) were distributed differently in the two viral groups (human- and plant-infecting), consistent with previous reports (30). However, differences emerge when allele targeting efficiency is considered, suggesting that the HLA system has been optimized through coevolution with viruses to recognize functionally important protein regions that are relevant to pathogen threats. Most interestingly, these groupings revealed patterns of efficiency variation over different HLA alleles and different viral families (Fig. 5), indicating a possible functional importance of the existence of long-lasting allelic lineages for the HLA-A but not HLA-B allele locus (24, 41) and possible evidence of specialization of the HLA loci. To assess this further, we investigated the distribution of efficiency scores for HLA-A and HLA-B loci according to viral genome composition and viral species, as described in the supplemental material (22, 28, 29, 62).

In our analysis of DNA viruses, we found that both HLA-A and HLA-B loci demonstrated a preference for targeting evolutionarily conserved regions of human DNA viruses relative to plant DNA virus species (P < 0.01) (Fig. 6). Moreover, there was also a striking general preference for HLA-A alleles to target conserved regions of human DNA viruses compared with HLA-B alleles (P < 10⁻¹³). This was most notable for the Herpesviridae and Adenoviridae (Table 1), with significant differences between HLA-A and HLA-B loci for 8 of 10 viral proteomes assessed. Since these viruses have emerged via distinct lineages through vertebrate evolution (13, 42), these observations are not readily explained by sequence similarity between these viral families. Yet the results are consistent with coevolutionary relationships between herpesviruses and adenoviruses and their human and ancestral primate hosts (13). Moreover, we found a linear correlation between HLA allele efficiencies for the analyzed human proteins and these DNA viruses (as shown for cytomegalovirus in Fig. 7A), as well as strong correlations among overall efficiency scores for the proteomes of DNA viruses (Fig. 7C). These observations are consistent with previous evidence that these viruses exploit host peptide mimicry as a means of host immune evasion (39, 60). We suggest that herpesviruses, which establish persistent but nonprogressive infection in the vast majority of human hosts, may evade the immune system (18) by exploiting the “holes” in the T-cell repertoire that are created by negative thymic selection. The similarity of herpesvirus and adenovirus proteome sequences to self-peptides would also be anticipated to induce more specific, less cross-reactive T-cell responses (10, 18), which may contribute to the progressive inflation of herpesvirus-specific T cells with a restricted T-cell receptor (TCR) repertoire during human ageing (33).

In contrast to results obtained for human protein and DNA viral protein targeting by HLA class I, we found that HLA-B alleles had higher efficiency scores for RNA viruses (Fig. 6B) (P < 10⁻¹⁶), with preferential targeting of evolutionarily conserved viral proteins by HLA-B noted for 23/34 (68%) of RNA viruses assessed (Table 1). Further scrutiny revealed a spectrum of targeting efficiency profiles across the range of HLA-B-virus pairs. At one end of the spectrum we found a general trend toward positive HLA-B efficiency scores for the human-adapted Paramyxoviridae (including respiratory viruses such as respiratory syncytial virus, parainfluenza virus, and metapneumovirus, as well as measles and mumps viruses) and the Picornaviridae (predominantly rhinovirus and enterovirus species). These RNA viruses exhibit a diverse range of “high-efficiency” HLA-B specificities (Fig. 5), consistent with the existence of a host-pathogen evolutionary relationship that is relatively specific for the HLA-B locus. Preferential and efficient targeting of these highly infectious, but nonpersistent, RNA viruses by the highly polymorphic HLA-B locus could also provide a potential mechanism for the observed rapid evolution of these viral species, characterized by the emergence of transient viral mutations (which may provide an HLA context-specific selection advantage) that are then purged by purifying selection (22, 52).

At the other end of this spectrum, the arboviral Flaviviridae are a dramatic counter-example to the general observation that HLA molecules preferentially target evolutionarily conserved proteomic regions. As shown in Fig. 5, most HLA-A and HLA-B molecules preferentially targeted nonconserved protein sequences from arboviral flaviviruses.

Targeting efficiency scores for two different protein groups are often correlated across different HLA alleles. For example, we found a strong negative correlation between allele efficiency scores for flaviviruses and the efficiency scores for human proteins or double-stranded DNA (dsDNA) viruses (Fig. 7A and C). The finding that HLA-A and HLA-B alleles that target conserved regions of human (self) proteins tend to target nonconserved regions of the dengue virus suggests that these viruses seek to reduce similarity to self-peptides to the extent that they are ignored by the immune system (18), thereby exploiting positive T-cell selection to their advantage. This evolutionary strategy also focuses HLA binding on relatively nonconserved viral proteome regions that are likely to have functional and genomic plasticity. We also find correlations between unrelated viruses (Fig. 7C), suggesting that these pathogens share common strategies to evade HLA-restricted immune surveillance. In contrast, negative correlations between unrelated viruses may indicate that these pathogens exploit the weaknesses in the immune surveillance created by overadaptation of the immune system to one virus or the other.

We found that these correlations cannot be explained by the sequence (dis)similarities among viruses: when the same correlation factors are computed for simulated HLA binding preferences, they are not as strong as they are for the 95 HLA alleles analyzed here, which cover true binding properties for over 90% of the European population. This again suggests a prominent role for HLA targeting efficiency in shaping HLA and pathogen coevolution (Fig. 7B; see also the supplemental material).

The findings presented thus far point to coadaptive relationships between HLA allelic diversity and human viruses but do not address the functional relevance of these observations. It is apparent from these analyses that relationships between virus species and host HLA diversity are highly specific, indicating specialized roles for HLA loci and for HLA allelic variants within these loci. We thus considered whether HLA targeting efficiencies for a given pathogen (with subsequent possible consequences for HLA-restricted immune responses) are associated with altered infectious disease outcomes. We focused on the specific examples of HIV disease progression, HIV viral load, and the incidence of dengue hemorrhagic fever (DHF).

Dengue hemorrhagic fever is a severe clinical manifestation of a secondary dengue flavivirus infection. Previous studies have suggested that the pathogenesis of this syndrome involves cross-reactive T-cell responses (15), which may be enriched in the context of low-affinity interactions between HLA class I and conserved viral peptides (as described for dengue viruses in Fig. 5; see also Fig. S7 in the supplemental material) (10, 18). In support of this model, we found that HLA alleles known to be associated with susceptibility to dengue hemorrhagic fever (11) appeared to be less likely to target evolutionarily conserved proteomic regions than alleles that confer resistance to this disease (P = 0.05) (Fig. 8).

We next analyzed the effect of HLA targeting efficiency for specific HIV-1 proteins on HIV-1 viral load in a study population of 191 HIV clade B-infected, treatment-naive individuals from the Western Australian HIV cohort for whom viral loads (plasma HIV RNA level) and full HLA typing were available (43). We computed HLA targeting efficiency for each individual by averaging binding scores over each patient's specific HLA-A and HLA-B repertoire, thus approximating the aggregate ability of patient-specific HLA alleles to differentiate between conserved and variable targets.

While the above analysis of viral targeting involved full proteome targeting efficiencies (Fig. 5), the underlying analysis of residue-specific binding and conservation scores calculated for individual viral proteins and for specific HLA alleles allows for a more focused exploration of the landscape of HLA-pathogen interactions based on targeting of particular individual proteins (see, for example, Fig. S4 to S6 in the supplemental material). Analysis of protein targeting efficiencies has the potential for broader use in the investigation of viral immunity relevant to natural infection as well as vaccine design and evaluation, particularly when indicators of functional immunity are available, as is the case in HIV-1 infection.

While overall HLA allele efficiencies toward HIV proteins correlated negatively (but not significantly) with log viral load (P = 0.24), certain individual protein targeting efficiencies showed significant correlations with viral load. For example, HLA-B locus efficiency in targeting Gag protein alone is more strongly correlated with viral load (r = −0.19; P = 0.009), consistent with experimental evidence that HLA-B-restricted cytotoxic T lymphocyte (CTL) responses to Gag epitopes play a significant role in determining the natural history of HIV infection (6, 35, 36, 54). We then analyzed the distribution of HLA allele efficiency scores for HIV-1 Gag according to their known associations with HIV progression (9, 19) and found that protective HLA alleles tend to rank more highly in targeting efficiency of conserved Gag proteomic regions (Fig. 9).

We also investigated the combined effect of efficiency scores for individual proteins and HLA loci on log viral load using multivariate regression. This analysis used the efficiency scores of all nine HIV proteins, taking into account proteasomal cleavage, and was performed on multiple test/train splits of the data (see supplemental material). Here, we found that HIV efficiencies alone account for 7.0% of the log viral load variance (P = 3.67 × 10⁻⁴; correlation, 0.27), gender alone explains 7.0%, and the patient's ethnic group alone shows no significant explanatory power. The HLA efficiencies and gender combined explain 11.2% of log viral load variance. The explanatory power of efficiency scores was not attributable to a single HLA allele or HLA supertype effect nor to confounding demographic effects such as gender and race (described in the supplemental material).

HLA binding is not the sole determinant of potential immune targets. The processing of intracellular antigens relies on relatively monomorphic and evolutionarily conserved proteins to optimize peptide cleavage prior to HLA binding (47), thus providing an additional mechanism for ligand selection. We therefore examined the potential influence of proteasomal cleavage on targeting efficiency by using the NetChop algorithm (46, 55). We analyzed targeting efficiency across human proteins and viral species using a restricted data set of peptides with a high probability of appropriate C terminus cleavage. We found that proteasomal cleavage restriction has indeed coevolved with HLA binding, and cleavage is also directed toward conserved targets. As shown in Fig. S7 in the supplemental material, the distribution of HLA targeting efficiencies remained similar to those identified in Fig. 5, indicating that relationships between HLA binding preference and evolutionary conservation are preserved among peptide targets that are selected via the antigen-processing complex.