Introduction
============

Pathogens are suspected to be one of the strongest selective forces in human evolution and the constant exposure to parasites over evolutionary time has likely contributed to the genetic variation found at a large number of genes within and among present day populations ([@msy116-B18]). In this light, the genes of the major histocompatibility complex (MHC) with their exceptional genetic diversity and immune function are a prime candidate to investigate the exact mechanisms through which pathogen-mediated selection has contributed to human evolution.

The MHC is a key component of the adaptive immune system common to all jawed vertebrates ([@msy116-B32]). In humans, it is known as a gene-dense region that spans ∼4 Mb on the short arm of chromosome 6. It comprises over 200 genes, many of which are involved in immunity ([@msy116-B4]). Among these genes, the classical MHC genes (also called human leukocyte antigen, HLA) encode for cell-surface glycoproteins with a key role in adaptive immunity ([@msy116-B28]; [@msy116-B81]). In cells infected by intracellular parasites, MHC class I molecules can present parasite-derived peptides to cytotoxic T lymphocytes (CTL). Upon recognition of these foreign peptides, the infected cells are destroyed. The MHC class II molecules present antigens, mainly derived from extracellular pathogens, on the surface of specialized antigen-presenting cells. The exposed peptides are recognized by helper T lymphocytes (T~H~ cells), leading to a complex cascade of specific immune responses ([@msy116-B28]; [@msy116-B29]; [@msy116-B56]).

The classical MHC genes are among the most polymorphic genes in the human genome and thousands of different alleles have been identified at some of these loci ([@msy116-B32]; [@msy116-B81]). This polymorphism is characterized by a remarkable sequence variation in the peptide-binding grooves of MHC molecules (i.e., the pocket where antigens are bound) ([@msy116-B62]; [@msy116-B67]) as well as an enhanced rate of nonsynonymous substitutions ([@msy116-B26]). MHC polymorphisms are often ancient and allele lineages whose origin predates species divergence are retained across multiple species, an observation described as transspecies polymorphism ([@msy116-B33]). The general action of balancing selection in enhancing both the rate of nonsynonymous substitutions in codons forming the peptide binding groove ([@msy116-B26], [@msy116-B27]) and the persistence of allelic diversity over extremely long time periods is strongly supported ([@msy116-B34][@msy116-B35]). However, the exact mechanisms of balancing selection are still disputed. Accordingly, three main mechanisms of pathogen-mediated selection have been suggested ([@msy116-B77]) which are potentially not mutually exclusive and may interact with one another: heterozygote advantage ([@msy116-B13]), rare-allele advantage ([@msy116-B6]), and fluctuating selection ([@msy116-B24]).

The heterozygote advantage was first proposed by [@msy116-B13]. Heterozygous individuals at MHC loci are assumed to present a broader range of pathogen-derived peptides than homozygotes, thus increasing the probability of triggering a specific immune response. They show increased resistance to pathogens, and are more likely to have higher relative fitness, resulting in an increased persistence of different MHC alleles in the population ([@msy116-B28]; [@msy116-B63]). The heterozygote advantage hypothesis has been further extended by taking into account the sequence level, leading to the idea of a divergent allele advantage ([@msy116-B64]; [@msy116-B83]). The high sequence divergence observed at MHC genes results in structural polymorphism that may impact the functional properties of MHC molecules. Heterozygous individuals with more divergent MHC allele combinations (i.e., larger number of amino acid differences along the sequence of the antigen-binding domains) are thought to encode glycoproteins that differ more in the repertoire of antigens they can bind. Those individuals may thus be able to present a wider array of antigens to immune effector cells, conferring an advantage against pathogen infections. In contrast, alleles more similar at the sequence level presumably exhibit more similar peptide binding specificities, thus leading to recognition of a lower overall number of peptides when co-occurring in a heterozygous individual ([@msy116-B41]).

Because of the extremely high number of pathogen proteins to which each host might be exposed throughout its lifetime, comprehensively measuring the relevant repertoire of MHC-bound peptides is impractical in humans and impossible in nonmodel species. Consequently, the divergent allele advantage hypothesis has been difficult to test. However, different measures of MHC sequence divergence are increasingly being used as a proxy for the potential MHC-bound peptide repertoire diversity, leading to correlative evidence that highlights how selection has favored the evolution of multiple MHC loci with divergent alleles in natural populations ([@msy116-B74]; [@msy116-B38]; [@msy116-B68]; [@msy116-B17]; [@msy116-B57]; [@msy116-B45]; [@msy116-B73]; [@msy116-B43]; [@msy116-B44]). In humans, the development of computational MHC antigen-binding prediction algorithms has enabled a more direct test of the divergent allele advantage. With this approach it has been shown previously that more divergent *HLA-DRB1* allele pairs experience less overlap in the antigenic peptides they can bind, that is, they are able to present a broader range of potential antigens ([@msy116-B41]), thus supporting the divergent allele advantage hypothesis at this locus. The investigation of the *DRB1* locus has been further extended by considering two distinct phylogenetic groups of alleles, denoted as group A and B ([@msy116-B87]). The same pattern of increased pathogen recognition capacity was observed only for those alleles that in the phylogenetic tree cluster together with primate alleles forming a polyphyletic group (group B) ([@msy116-B39]). Recently, a mechanism of joint divergent asymmetric selection acting on *HLA-A* and *B* as a whole was suggested, which has potentially evolved to counter-balance the lack of diversity at individual HLA loci often found in small-sized and isolated human populations ([@msy116-B9]).

In order to evaluate the divergent allele advantage hypothesis more systematically across all key classical human MHC genes, we here investigated the relationship between sequence divergence and peptide binding properties for three class I genes (*HLA-A*, *-B*, *-C*) and two class II genes (*HLA-DRB1* and *-DQB1*). Focusing on "*common*" alleles for each locus, as defined by the CWD catalogue ([@msy116-B49]) ([supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), we evaluated different estimates of amino acid sequence divergence as proxies for the functional divergence among different alleles. Functional divergence was characterized by allele-specific computational binding prediction for a broad range of representative human pathogens ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Considering a larger and more comprehensive data set of pathogen-derived peptides compared with the set of pathogenic peptides that has been used in previous studies, we were also able to describe the functional features of the divergent allele advantage, by investigating the differential pattern of antigenic presentation between MHC class I and class Il loci. Finally, the frequency distribution of HLA allele pools was investigated in several European populations in order to explore ongoing selection for divergent MHC alleles in modern humans.

Results
=======

Functional Characterization of Common Human MHC Alleles
-------------------------------------------------------

The set of 232 proteins from a broad collection of relevant human pathogens (*N* = 27), including macroparasites, bacteria, and viruses, resulted in a total of 118,097 unique pathogen-derived peptides. These peptides are meant to represent a comprehensive repertoire of potential antigens to which humans may have been exposed to throughout their evolutionary history, and which may thus have contributed to the exceptional MHC diversity that we see in present-day human populations.

The number of alleles defined as *common* varied among the different MHC loci: *HLA-A*: 63, *HLA-B*: 123, *HLA-C*: 40, *HLA-DRB1*: 73, *HLA-DQB1*: 21 ([supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), reflecting general differences in allelic diversity among the loci ([@msy116-B81]). The proportion of peptides predicted to be bound by a given allele varied substantially within and among the different loci ([table 1](#msy116-T1){ref-type="table"} and [supplementary fig. S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). For each locus, the proportions of common (shared among different alleles) and private (allele-specific) peptides bound by each allele are reported in [supplementary figure S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online. The overlap in bound peptides between different loci was significantly higher for class I genes (*A* ∩ *B* = 8,253, *B* ∩ *C* = 5,930, *A* ∩ *C* = 5,692, *A* ∩ *B* ∩ *C* = 3,545) than for class II gene (*DRB1* ∩ *DQB1* = 133) (χ^2^ test, *P* \<0.001) ([supplementary fig. S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Interestingly, we also found a remarkable amount of bound peptides that was shared between class I and class II genes (*A* ∩ *DRB1* = 1,377, *A* ∩ *DQB1* = 385, *B* ∩ *DRB1* = 933, *B* ∩ *DQB1* = 403, *C* ∩ *DRB1* = 689, *C* ∩ *DQB1* = 251) ([supplementary fig. S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). However, since the binding prediction algorithms do not account for different antigen processing pathways, it remains to be explored whether in reality these shared peptides are presented by alleles from both classes of MHC. The number of combined peptides bound by any two alleles of a given locus (equivalent to a heterozygous genotype) differed significantly among the genes (Kruskal--Wallis test, *P* \< 0.001) ([supplementary fig. S5](#sup1){ref-type="supplementary-material"}*a*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Generally, class II genes showed a lower median number of peptides bound by allele pairs (*HLA-DRB1* = 587, *HLA-DQB1* = 446) compared with class I genes (*HLA-A* = 3,936, *HLA-B* = 3,753, *HLA-C* = 3,349). Table 1.Proportion of Bound Peptides across the Five Classical MHC Loci.LocusProportion of Peptides Bound By At Least One AlleleProportion of Bound Peptides Per AlleleMinMaxMedian (95% CI)*HLA-A*0.1850.0140.0740.017 (0.0159, 0.0178)*HLA-B*0.1920.0120.0180.017 (0.0168, 0.0172)*HLA-C*0.0790.0140.0190.017 (0.0165, 0.0171)*HLA-DRB1*0.0250.0020.0050.003 (0.0026, 0.0031)*HLA-DQB1*0.0110.0010.0040.002 (0.0015, 0.0027)[^1]

Sequence-Based Divergence Parameters and Functional Differences among MHC Alleles
---------------------------------------------------------------------------------

A growing number of studies are investigating the fitness consequences of MHC allele divergence both in natural populations and in model species by using some estimate of sequence dissimilarity as a proxy for functional divergence among alleles. These estimates range from simple nucleotide differences (only partially relevant at the protein level) to sophisticated methods that take into account different physicochemical properties of amino acids at the protein sequence level.

For instance, the pairwise amino acid p-distance simply counts the relative number of differences along the amino acid sequence, but information on amino acid properties and relationships are not incorporated, and all nonidentical amino acids are treated as equivalent ([@msy116-B23]). However, as the substitution rate usually varies among amino acid site, methods including information about different mutation rates for each amino acid as well as scores that take into account residue-specific properties have been introduced ([@msy116-B53]). DayHoff ([@msy116-B12]) and JTT ([@msy116-B30]) are two examples of the most popular methods mainly used to investigate mutational trajectories and evolutionary distances between amino acids. Additionally, quantitative measures of pairwise distance have been developed, in which the physicochemical properties of the amino acids, and thus the functional similarity between sequences are considered. Among the different physicochemical features, the molecular volumes of amino acid residues might be particularly meaningful for the question whether a peptide fits into the various pockets of the peptide-binding groove of an MHC molecule. Grantham ([@msy116-B22]) and Sandberg ([@msy116-B72]) distances are two examples of sequence-based measures of sequence divergence where the molecular volume of the different amino acids is taken into account.

So far, a comprehensive evaluation has been lacking as to which sequence parameters are most suitable, that is, most strongly correlated with functional divergence. Thus, in order to identify the most relevant sequence-based parameter, the allele-specific functional binding properties at each MHC locus were correlated with five different, commonly used measures of sequence divergence: pairwise amino acid p-distance ([@msy116-B23]), DayHoff ([@msy116-B12]), JTT ([@msy116-B30]), Grantham ([@msy116-B22]) and Sandberg ([@msy116-B72]).

For each possible allele pair at a given HLA locus, we calculated the total number of unique pathogen-derived peptides (obtained by computational binding prediction as described earlier) bound by both alleles together (meant to reflect the MHC-presented antigen repertoire conferred by a heterozygote genotype). The pairwise number of bound peptides was then correlated with the sequence divergence of the two given alleles, estimated by the different measures. As expected, the five measures of genetic distance were highly correlated with each other ([supplementary fig. S6](#sup1){ref-type="supplementary-material"} and [table S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Accordingly, the correlation values between allele divergence and the combined number of bound peptides by all possible allele pairs at each locus under investigation were largely consistent across the different parameters of sequence divergence ([supplementary table S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Nevertheless, despite similar correlation values, a rank analysis across the five human MHC gene (*HLA-A*, *-B*, *-C*, *-DRB1*, and *-DQB1*) revealed that the Grantham distance measure consistently ranked at the top, that is, showed the strongest correlation values ([table 2](#msy116-T2){ref-type="table"}). Therefore, for subsequent analyses, we focused on one parameter only, the Grantham distance. Table 2.Rank analysis between different measures of sequence divergence.TotalExtracellularIntra-ExtraIntracellularAverage TauRankAverage TauRankAverage TauRankAverage TauRankP-distance0.32120.26530.29720.2972Dayhoff0.31330.25740.28550.2972JTT0.30950.25650.28640.2874Grantham0.32710.27210.29810.3051Sandberg0.31040.26920.28930.2825[^2]

Relation between Sequence Divergence and Functional Divergence
--------------------------------------------------------------

The median distance between allele pairs of a given locus differed significantly between the five genes (Kruskal--Wallis test, *P* \< 0.001) with class II genes showing larger median Grantham distances (*HLA-DRB1* = 9, *HLA-DQB1* = 11) compared with class I genes (*HLA-A* = 7, *HLA-B* = 8, *HLA-C* = 5) ([supplementary fig. S5](#sup1){ref-type="supplementary-material"}*b*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

For the sake of completeness, we first used our data to confirm the intuitive assumption of the heterozygote advantage hypothesis that allele pairs (representing heterozygous genotypes) together bind a larger number of peptides than single alleles (representing homozygous genotypes). This was generally true for all five loci (Kruskal--Wallis test, *P* \< 0.001; [supplementary fig. S7](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), even though there were some rare cases where certain alleles alone bound more peptides than certain allele combinations, suggesting interesting variation in peptide promiscuity among alleles. This general result is in line with a large body of empirical studies showing higher pathogen resistance for MHC heterozygotes ([@msy116-B10]; [@msy116-B63]).

Subsequently, we focus all our analyses on allele pairs with two different alleles (reflecting heterozygous genotypes). According to the divergent allele advantage hypothesis, we expect that the number of peptides bound by heterozygote genotypes increases with increasing sequence divergence between the two given HLA alleles. Following this expectation, all five HLA genes revealed a significant positive correlation between the pairwise genetic distance and the combined number of bound peptides across all possible allele pairs ([fig. 1](#msy116-F1){ref-type="fig"} and [table 3](#msy116-T3){ref-type="table"}). Interestingly, for *HLA-A*, *-B*, and *-DRB1*, the rate at which the number of bound peptides increases in response to larger allele divergence appears to slow down after a certain point, seemingly approaching a maximum. This can be explained by the fact that, for some loci, even allele pairs with only intermediate sequence divergence do not share any bound peptides anymore ([supplementary fig. S8](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). As the combined number of bound peptides cannot be larger than the sum of peptides bound be each allele, as soon as zero overlap is reached, any further sequence divergence cannot increase the combined number of bound peptides any further. This suggests that alleles at some loci can diverge functionally with only a small number of sequence changes, possibly at sites located in the peptide binding region. Table 3.Divergent Allele Advantage and Different Pathogen Groups.HLA-AHLA-BHLA-CHLA-DRB1HLA-DQB1Pathogen GroupsTau*P*~adj~Tau*P*~adj~Tau*P*~adj~Tau*P*~adj~Tau*P*~adj~Total0.361\<0.0010.397\<0.0010.507\<0.0010.157\<0.0010.2100.001Extracellular0.345\<0.0010.192\<0.0010.392\<0.0010.130\<0.0010.303\<0.001Intra--Extra0.351\<0.0010.289\<0.0010.475\<0.0010.166\<0.0010.2100.024Intracellular0.293\<0.0010.377\<0.0010.544\<0.0010.137\<0.0010.172\<0.001[^3][^4]

![In silico evidence for divergent allele advantage across five classical human MHC genes. Correlation between pairwise genetic distances reported as Grantham distance (*x* axes) and number of bound peptides (*y* axes) counted for all possible pairs of common HLA alleles. Each dot represents an allele pair. Binding prediction analyses performed on the complete data set of pathogen proteins (*n* = 232). Linear model (red line) and smoothed lowess curve (dashed blue line), describing the association between the combined number of bound peptides and pairwise Grantham sequence divergence. Note the different axis scales.](msy116f1){#msy116-F1}

Furthermore, a significant negative correlation could be observed between pairwise genetic distance and the proportion of shared peptides (peptides bound by both alleles of a given combination; [supplementary fig. S8](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), revealing a decreasing proportion of peptides shared between more divergent alleles. Of note, for the two class II loci, the correlations between genetic distance and peptide sharing were stronger than between genetic distance and the total number of bound peptides. This is owing to the fact that the latter measure includes additional variation from differences in the size of the bound peptide repertoire among HLA class II alleles ([supplementary fig. S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), an allele-specific property that is independent of the divergence between alleles.

Phylogenetic analysis of the *HLA-DRB1* gene has revealed two subgroups of allelic lineages: a human-specific monophyletic group (Group A), and a polyphyletic group with primates (Group B) ([@msy116-B87]). It has been proposed that group A and B allele lineages have evolved with contrasting binding capacity, and only the alleles from the polyphyletic group B showed increased presentation of pathogen peptides with increasing sequence divergence ([@msy116-B39]). In contrast with previous findings, our binding prediction analysis revealed a significant positive correlation for the whole set of *DRB1* alleles ([fig. 1](#msy116-F1){ref-type="fig"} and [table 3](#msy116-T3){ref-type="table"}) as well as for both groups of alleles separately ([supplementary fig. S9](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). This discrepancy to the earlier results might be due to the much larger and more comprehensive data set of pathogen-derived peptides used in our analysis (here 118,097 peptides vs. 265 peptides in the previous study).

Antigen Processing and Different Origins of Pathogenic Peptides
---------------------------------------------------------------

Up to this point, our analysis treated all pathogen peptides as equally likely targets for each given MHC locus. However, in reality some of those peptides will never be in contact with certain MHC molecules, due to the different processing pathways by which eukaryotic cells degrade proteins: the proteasome and lysosomal proteases. Peptides resulting from proteasome degradation, generally derived from intracellular proteins, are presented by MHC class I molecules, whereas peptides presented by MHC class II molecules are usually of extracellular origin and processed through lysosomal protease degradation in antigen-presenting cells ([@msy116-B29]). If the exceptional sequence divergence among MHC alleles evolved at least partly as a consequence of pathogen-mediated selection for divergent alleles, we would expect to see a stronger signature of selection at a given MHC locus when focusing on biologically meaningful pathogens that are likely to actually be encountered by a given MHC molecule. That is, we expect a stronger correlation between sequence divergence and functional divergence when only focusing on peptides originating from the locus-specific antigen-processing pathway. We thus divided the pathogen proteins used for binding prediction analysis into three groups. The three groups were based on their agent's lifestyle in the host: "extracellular" (*n* = 58 proteins), "intracellular" (*n* = 100) and a third group of pathogen proteins belonging to those agents whose life cycle involves both intracellular and extracellular stages inside the host ([@msy116-B75]), here named as "intra-extracellular" (*n* = 75) ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Within each of the three groups of pathogen proteins, and for all five investigated HLA genes, we again observed a significant positive correlation between the pairwise sequence divergence and the combined number of bound peptides across all possible allele pairs. Interestingly however, in some cases, the strength of correlation differed among the three pathogen groups: at two of the MHC class I genes, *HLA-B* and *HLA-C*, a stronger positive correlation was observed for the group of intracellular pathogens, compared with extracellular pathogens, while this bias was not observed at *HLA-A* ([fig. 2](#msy116-F2){ref-type="fig"} and [table 3](#msy116-T3){ref-type="table"}). Conversely, a stronger correlation for peptides derived from extracellular pathogens was detected at one MHC class II locus, *HLA-DQB1*; while correlation values between pairwise sequence divergence and the combined number of bound peptides at the *HLA-DRB1* locus were comparable across the three groups of pathogen proteins ([fig. 2](#msy116-F2){ref-type="fig"} and [table 3](#msy116-T3){ref-type="table"}). In the three cases where a bias across the groups of pathogen proteins was detected, the intra-extracellular proteins showed intermediate correlation values between the intra- and the extracellular proteins ([table 3](#msy116-T3){ref-type="table"}). However, permutation of the proteins among the three groups revealed that only for *HLA-B* and *HLA-C*, the observed difference between extracellular and intracellular correlation values were indeed larger than expected by chance (both *P* \< 0.001; [supplementary fig. S10](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), while the other three loci did not show a statistically significant deviation from random expectations.

![Different origins of pathogenic peptides. Correlation between pairwise genetic distances reported as Grantham distance (*x* axes) and number of bound peptides (*y* axes) calculated for all possible pairs of common HLA alleles. Each dot represents an allele pair. Binding prediction analyses performed considering proteins within three groups of pathogens: extracellular (orange, *n* = 57 proteins), intracellular (green, *n* = 100), and intra-extracellular (purple, *n* = 75).](msy116f2){#msy116-F2}

Our main analyses were performed considering pathogen-derived peptides that are all of the same length (9 aa). However, MHC class II molecules allow binding longer peptides than class I. Thus, for alleles at the two class II loci, binding prediction of 15mer peptides from the same set of pathogen proteins were considered. This analysis showed stronger correlations for *DRB1* and weaker correlations for *DQB1*, compared with the 9mer predictions, but overall support the main conclusions ([supplementary table S5](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Furthermore, to test if our results were sensitive to the choice of the binding threshold, we additionally repeated the primary analysis using a different binding threshold (%rank of 0.5, indicating strong binding). The comparable results suggest that our main conclusions hold across a range of established binding thresholds ([supplementary table S6](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

Distinct Amino Acid Composition among Pathogen Groups
-----------------------------------------------------

The observed bias in the correlation between sequence divergence and functional divergence (peptide binding) toward a specific group of pathogens suggests distinct differences in the peptide repertoires among these groups. This could be due either to certain group-specific peptide sequences or to a more general difference in the amino acid composition of proteins among the pathogen groups. Amino acid usage has changed over evolutionary time in different species, and proteins have evolved in terms of physico-chemical and structural properties, reflecting adaptations to specific environmental conditions ([@msy116-B7]; [@msy116-B79]). Intriguingly, intra- and extracellular environments exhibit significant differences, including different pH value ([@msy116-B11]) and availability of different nutrients ([@msy116-B20]; [@msy116-B60]; [@msy116-B71]). It thus appears plausible that intra- and extracellular pathogens may have evolved proteomes with distinct amino acid compositions.

In order to test if the observed bias in correlation values across the groups of pathogen proteins was the results of group-specific peptide sequences or due to more general differences in the amino acid composition within each group, we created four different data sets of artificial proteins. These four data sets were then analyzed in the same way as above, again assessing the strength of correlation between allele divergence and functional divergence (here based on bound peptides from the artificial proteins) for the three pathogen groups and across the five HLA loci. First, amino acids forming each pathogen protein sequence were randomly shuffled, maintaining the same amino acid composition of a given protein, but changing its actual sequence. If the stronger correlation with intracellular pathogen proteins by *HLA-B* and *HLA-C* were due to group-specific peptide sequences, we would expect this bias to disappear when reshuffling the protein sequences. However, correlation values resulting from shuffled proteins did not differ substantially from the true observed correlation values obtained with the real data set of pathogen proteins ([supplementary table S7](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). For the second set of artificial proteins, we created random protein sequences but maintained amino acid frequencies as they occurred within each group of pathogen proteins. Again, correlation values did not differ substantially from the observed true correlation values obtained in the initial test and the specific bias across the three groups was still observed ([supplementary table S7](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). For the third set of artificial proteins, we again created random protein sequences but this time maintained amino acid frequencies as they occurred in the whole data set of pathogen proteins. While a general positive correlation was also observed in this data set, the specific bias across the groups of pathogen protein was not detected anymore ([supplementary table S7](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Finally, the amino acid composition computed from UniProtKB/Swiss-Prot data bank ([@msy116-B19]; [@msy116-B8]) was used to create the fourth set of artificial proteins. Again, the specific bias across the three groups of pathogen proteins was not detected anymore ([supplementary table S7](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Thus, the observed bias in correlation values persisted only when amino acid frequencies mirrored the specific frequencies observed within each group of pathogen proteins. One of the possible explanations of our results could be that the observed stronger correlation of *HLA-B* and *HLA-C* alleles with peptides from intracellular pathogens is not due to specific peptides but is the results of adaptation of MHC alleles to the differences in amino acid composition between groups of pathogens. To further explore this hypothesis, a nonparametric multivariate analysis of variance was performed to quantify the similarity among proteins with regard to their amino-acid composition. The two groups of intracellular and extracellular pathogen proteins indeed differed significantly in their amino acid composition (PerMANOVA test, *P* \< 0.001), with 9% of the total variance associated with the divisions in intracellular and extracellular proteins (R\^^2^ estimate from PerMANOVA) ([fig. 3](#msy116-F3){ref-type="fig"}). Accordingly, when average amino acid compositions were compared between the two groups of pathogen proteins, significant variations in the mean amino acid composition were observed for specific amino acids (one-way ANOVA, *P* \< 0.05) ([supplementary fig. S11](#sup1){ref-type="supplementary-material"} and [table S8](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). The observed differences in amino acid composition could be linked to glycosylation patterns, which differ between extracellular and intercellular peptides ([@msy116-B51]). N-glycosylation is one of the forms of protein glycosylation in eukaryotic organisms which is mainly targeting extracellular and secreted proteins. It has been shown that N-glycosylation sites are specific to the consensus sequence Asn-Xaa-Ser/Thr and that the presence of proline between Asn and Ser/Thr inhibits N-glycosylation ([@msy116-B3]). Accordingly, in our analysis proteins of extracellular pathogens show low proline concentration which is instead prevalent in proteins of intracellular pathogens. These exploratory analyses suggest that the amino acid composition might be different between the two groups of intracellular and extracellular pathogen proteins and that MHC alleles might have potentially adapted their binding specificities accordingly, at least at the *HLA-B* and *HLA-C* loci. However, further research is necessary to validate this conclusion and to exclude other potential causes, such as taxonomy, driving the observed difference in amino acid composition.

![Multidimensional scaling plot of amino acid composition in intracellular and extracellular pathogen proteins. Multidimensional scaling (MDS) based on amino acid frequencies indicates similarity in amino acid composition among individual proteins (dots). Intracellular proteins (*n* = 100) are reported in green while extracellular proteins (*n* = 57) in orange. MDS enables a standardized unit-less representation of variation among data points in 2D space (along perpendicular axes MDS1 and MDS2): location of proteins within the plot is indicative of potential bias toward specific amino acids (blue characters in one letter code), proteins with more similar amino acid composition are displayed closer to each other. The dashed circles indicate 95% confidence intervals for each group. Stress for 2D representation: 0.21.](msy116f3){#msy116-F3}

Population Frequency of Divergent HLA Alleles
---------------------------------------------

The above-described results support historical pathogen-mediated selection through divergent allele advantage at the human MHC. However, we were also interested in exploring whether the divergent allele advantage was still maintaining diverse HLA allele pools in present day human populations. We hypothesized that, under the divergent allele advantage, alleles that on an average yield a more divergent genotype (conferring higher fitness) when paired with another allele in a heterozygote individual, would be selected for and thus exhibit higher frequencies in a given population. A similar observation had been made in the allele pool of a social marine mammal whose reproductive success is partly predicted by the divergence of its MHC genotype ([@msy116-B44]). In order to test this hypothesis across the five classical HLA loci, we calculated for each allele of a given locus the average pairwise amino acid sequence divergence to the most common alleles at this locus (≥5% allele frequency, representing alleles most likely to be forming a heterozygote with the allele in question). This average sequence divergence was then correlated with the allele's population frequency. For this analysis, we focused again on alleles defined as "*common*" in the CWD catalogue ([@msy116-B49]), assuming that very rare alleles are more susceptible to stochasticity and neutral demographic processes. In line with our expectation, a significant positive correlation between the average pairwise sequence divergence of an allele and its population frequency was observed for a number of HLA genes (*HLA-B*, *-C*, and *-DRB1*) across different European populations (USA European Caucasian, German, and Polish) ([fig. 4](#msy116-F4){ref-type="fig"} and [supplementary table S9](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

![Population frequency of divergent HLA alleles. Correlation between the average Grantham pairwise divergence to the most common alleles and the allele frequency in the USA NMDP European Caucasian (*N* = 1,242,890), German (*N*= 39,689), and Polish (*N*= 20,653) populations, for four classical HLA loci with available allele frequency data in AlleleFrequencies.net. Significant associations that persisted after Bonferroni correction across populations are reported with a solid line, while dashed lines indicate associations that are only nominally significant (*P* \< 0.05 before Bonferroni correction; for exact values see [supplementary table S9](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).](msy116f4){#msy116-F4}

Discussion
==========

Here, we used computational antigen-binding prediction on a large data set of potentially antigenic pathogen peptides to investigate whether pairs of highly diverged MHC alleles together bind more different antigens than more similar alleles. Such an association is predicted by the *divergent allele advantage* hypothesis ([@msy116-B83]) and had previously been investigated only for the *HLA-DRB1* gene in humans ([@msy116-B41]; [@msy116-B39]). The observed positive correlation between the genetic distance of two alleles and the combined number of peptides they bind together confirmed and extended the predictions of the divergent allele advantage for all five investigated human MHC genes. These results support the hypothesis that enhanced sequence diversity between alleles in a heterozygous MHC genotype increases the range of potential MHC-presented peptides, thus raising the chance to recognize specific antigens and consequently enhance immune surveillance.

Our observation that *HLA-B* and *HLA-C* (and potentially *HLA-DQB1*) exhibit the strongest associations when considering antigens originating from their most plausible targets is intriguing and lends further support to the biological relevance of this mechanism. It might indicate that the exceptional sequence divergence commonly observed among alleles of a given MHC locus has evolved specifically in response to selection by pathogens that are processed through the major protein degradation and antigen-presentation pathways that this given locus is associated with: alleles at *HLA-B* and *-C* loci have apparently evolved to bind specifically peptides derived from intracellular pathogens, while alleles at the *HLA-DQB1* locus may have evolved to bind a wider array of peptides from extracellular pathogens. The fact that we did not observe such a pathogen group-specific bias at the *HLA-A* and *HLA-DRB1* loci might indicate either that divergent allele advantage has not played a significant role in their evolution, or that they are less specific with regard to the pathogen origin of the peptides they present. It has indeed been shown that autophagy of intracellular components can promote the presentation of endogenous antigens by MHC class II molecules ([@msy116-B61]; [@msy116-B46]; [@msy116-B55]; [@msy116-B70]). Furthermore, several studies have reported a potential role for DRB1 molecules in viral infections ([@msy116-B52]). For instance, HLA-DR variants have been associated with spontaneous clearance of HBV and HCV infections ([@msy116-B80]; [@msy116-B54]), with protective effect against dengue shock syndrome (DSS) development ([@msy116-B58]) and with HIV suppression ([@msy116-B50]). In this light, the lack of a particular bias by DRB1 alleles toward either of the pathogen groups may indicate that this locus evolves under selection by both intra- and extracellular pathogens.

In addressing the frequency of divergent HLA alleles in different human populations, we observed in some human populations, and for specific HLA loci, significant correlations between an allele's population frequency and its average pairwise sequence divergence. These results might suggest still ongoing selection for divergent HLA genotypes, at least in some modern human populations, possibly depending on population-specific differences in historical pathogen communities. However, allelic age may also contribute to the observed pattern, as, in principle, older alleles are both more likely to have reached high frequencies, even under neutrality, and to have accumulated more point mutations (thus being more divergent). On the other hand, HLA genes are known to undergo frequent recombination and gene conversion events, yielding novel alleles with high divergence from their origin at the very start (i.e., at low frequency). It is thus unclear to what extent novel alleles contribute to the observed pattern, warranting further research to explore the effect of genetic drift on the frequency of divergent HLA alleles. Furthermore, HLA alleles that on an average form more divergent allele combination, and which have been maintained in the population because of their increased capacity in presenting pathogen-derived peptides, might also be advantageous in case of newly emerging and fast-evolving pathogens (i.e., HIV).

Humans share similar MHC allelic lineages with closely related species ([@msy116-B33]; [@msy116-B40]). This observation is a typical feature of MHC genes, compatible with the theory of transspecies evolution: ancestral lineages present in the common ancestor are inherited through successive speciation events, persisting over long periods of time ([@msy116-B35]). The ancestral and highly diverged MHC variants are assumed to be adaptive and selectively maintained as a polymorphism by balancing selection ([@msy116-B26]; [@msy116-B41]). Recently, the role of adaptive introgression has been proposed to contribute to the exceptional level of polymorphism at the MHC ([@msy116-B1]; [@msy116-B84]). So far, explanations for the maintenance of introgressed MHC alleles have largely relied on the idea that such alleles were somehow locally adapted and thus beneficial. For instance, it has been suggested that modern humans might have maintained introgressed archaic HLA variants because they conferred an advantage against local pathogens ([@msy116-B1]). However, another explanation appears also plausible: MHC alleles from another species are, on average, likely to have diverged significantly from the species' own allele pool. Thus, any allele that introgresses from another species is likely to lead to highly divergent MHC genotypes. Following the divergent allele advantage hypothesis (and our results), such introgressed alleles should then confer a significant advantage and should consequently be selected for in the new species. This scenario would easily explain the maintenance of introgressed MHC alleles, but further research is necessary to support this hypothesis.

While the present analysis focuses exclusively on the divergent allele advantage, in reality, selection at MHC genes is a dynamic process that involves additional mechanisms apart from the divergent allele advantage. Conceptually, the divergent allele advantage can be considered a quantitative mechanism, which works independent of specific pathogen species or strains. It can act over long evolutionary time scales, promoting the maintenance of ancient allelic lineages in natural populations ([@msy116-B41]) and facilitating immunity against the constant simultaneous barrage by many different pathogens. In contrast, negative frequency-dependent selection (NFDS) is a qualitative mechanism in which specific alleles can be selected by specific pathogens ([@msy116-B76]; [@msy116-B42]). This mechanism likely works on a shorter time scale, for instance affecting MHC evolution in humans in very recent history ([@msy116-B48]; [@msy116-B36]). Both mechanisms, the divergent allele advantage and NFDS might also act in parallel, but at different time scales, creating an intriguing combination of shared polymorphism but distinct allele pools among populations and possibly even species ([@msy116-B47]). Local adaptation plays another significant role in MHC evolution and might modulate the effect of the above mechanisms ([@msy116-B15]). The simultaneous action of these additional mechanisms might occasionally mask the effect of the divergent allele advantage and potentially explain the only sporadic evidence for this mechanism in the population frequency analysis reported here. Nevertheless, our results strongly support the divergent allele advantage as a meaningful quantitative mechanism through which pathogen-mediated selection contributes to the evolution of MHC diversity.

Materials and Methods
=====================

MHC Loci and Alleles Included in Analyses
-----------------------------------------

Five key classical human MHC genes (*HLA-A*, *-B*, *-C*, *-DRB1*, and *-DQB1*) were analyzed in this study. Alleles at each locus were defined at second field (four-digit) resolution and only alleles annotated as "*common*" in the CWD catalogue ([@msy116-B49]) were included in the analyses. The allele annotation "*common*" in the CWD catalogue does not specifically indicate a high population frequency but more the extent and quality of documentation available for the given allele. This category indicates that there is universal agreement about the identity of this allele because it has been observed in multiple populations and there is sufficient data for robust frequency estimation ([@msy116-B49]). These criteria resulted in the analysis of 63 alleles for *HLA-A*, 123 for *HLA-B*, 40 for *HLA-C*, 73 for *HLA-DRB1*, and 21 for *HLA-DQB1* ([supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

Pathogen Proteins
-----------------

Binding prediction analyses were performed on a data set of representative human pathogen proteins. Pathogens were selected from the Gideon database ([@msy116-B5]) based on the following criteria: a global distribution, a potential for high mortality and/or morbidity, and a significant impact over the course of human history ([@msy116-B85]). The rational for these criteria was that such pathogens are likely to have contributed significantly to human evolution in general and to the evolution of MHC genes in particular. [@msy116-B85] provided a comprehensive list of infectious diseases with the greatest evolutionary and historical significance. From that list, we have taken the majority of pathogens in our data set. However, to assess mortality and morbidity, epidemiological data were also collected from two published reports: the Annual report of the European Centre for Disease Prevention and Control ([@msy116-B16]) and the WHO Global Health Estimates ([@msy116-B86]). First, pathogens with the highest current mortality were included. However, not just mortality, but also nonfatal morbidity can be historically and evolutionarily significant. Indeed, morbid pathogens can reduce the fitness of their host in different ways (e.g., by increasing the sterility), thus pathogens considered morbid were also included. Finally, eradicated pathogens known to be important in human history were taken into account. Here, we used protein sequences of present day pathogens to explore signatures of historical selection, even though ancient pathogen strains might have differed slightly in their antigen repertoires. While we do not expect an effect on the general patterns observed here, it might be interesting to explore subtle differences in future work. We further aimed for a balanced representation of different groups of pathogens (i.e., viruses, bacteria, parasites). Based on these criteria, we identified 27 pathogens (10 viruses, 10 bacteria, 7 macroparasites) that were classified into three groups: extracellular, intracellular, and intra-extracellular, based on their primary environment in the human body ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Then, for the selected pathogens, amino acid sequences of 232 pathogen proteins (8.5 ± 5.8 per pathogen) known to be antigenic ([@msy116-B82]) and/or likely exposed to the host immune system (mostly secreted and surface proteins) ([@msy116-B66]) were obtained from GenBank (for accession numbers see [supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

Peptide Binding Prediction Algorithms
-------------------------------------

Computational antigen-binding prediction algorithms for MHC molecules were used to determine pathogen peptides potentially bound by the MHC alleles under investigation. Binding prediction was computed for all alleles at each of the five human MHC genes. Furthermore, as prediction analysis are likely to be more accurate for the core of the binding groove, which is known to be nine residues long and contributes the most to the recognition of the antigens, binding prediction was performed considering all possible 9mer pathogen-derived peptides. The data set of 232 representative human pathogen proteins described above resulted in a total of 118,097 unique pathogen-derived 9mer peptides that were analyzed using two different algorithms: NetMHCpan (v2.8) ([@msy116-B25]) for the alleles at class I loci (*HLA-A*, *-B*, *-C*) and NetMHCIIpan (v3.0) ([@msy116-B31]) for the alleles at class II loci (*HLA-DRB1*, *-DQB1*). For alleles at the two class II loci (*HLA-DRB1* and *HLA-DQB1*), we repeated the binding prediction analysis considering all possible 15mer pathogen-derived peptides. The predicted binding affinity between pathogen peptides and MHC molecule variants (defined in nanomolar IC50, i.e., half maximal inhibitory concentration) are ranked by the respective software, based on comparison with a large pool of naturally occurring peptides, and a rank percentage score (%rank) is assigned to each peptide. To define "bound" peptides, we used the default %rank threshold of 2, which includes weak and strong binders. All analysis were also repeated using another established binding threshold (%rank of 0.5) which includes only strong binders. The allele *HLA-A*\*30:04 was predicted to bind about four times as many peptides as the other 62 *HLA-A* alleles ([supplementary fig. S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online) and was thus excluded as an outlier from subsequent analysis in order to prevent distortion of results. The binding prediction analyses were performed first on the complete data set of pathogen proteins (*n* = 232), and then considering proteins within three groups separately: extracellular (*n* = 58), intracellular (*n* = 100), and intra-extracellular (*n* = 75).

Sequence Divergence
-------------------

Allele divergence was computed on the same set of alleles used in the binding prediction analysis reported in [supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online. Protein sequences of HLA alleles were obtained from IMGT/HLA database ([@msy116-B69]). Exons forming the variable region in the peptide binding groove (i.e., exon 2 and 3 for class I alleles and exon 2 for class II alleles) were selected following the annotation obtained from Ensemble database ([@msy116-B2]). Amino-acid sequence alignments were performed using MUSCLE ([@msy116-B14]), and sites containing alignment gaps at the beginning or the end of sequences were removed. Genetic distances between alleles for all possible allele pairs at each locus were determined removing missing sites in pairwise comparisons and using five different pairwise parameters of allele divergence: p-distance ([@msy116-B23]), DayHoff ([@msy116-B12]), JTT ([@msy116-B30]), Grantham ([@msy116-B22]), and Sandberg ([@msy116-B72]). Pairwise amino acid p-distance, DayHoff and JTT distances were calculated in MEGA 7 ([@msy116-B37]). Grantham and Sandberg sequence distances were calculated using a custom Perl script that required two input files: a FASTA file with aligned HLA alleles and a specific amino acid distance matrix. Grantham amino acid distance matrix was constructed from [@msy116-B22]. Sandberg amino acid distance matrix was calculated based on Euclidian distances between all 20 amino acids, using the Euclidian distance method in R version 3.4.1 ([@msy116-B65]) according to the five physicochemical z-descriptors described in [@msy116-B72]: z1 (hydorphobicity), z2 (steric bulk), z3 (polarity), z4, and z5 (electronic effects). Our perl script (together with the Grantham amino acid similarity matrix) is freely available for download from SourceForge (<https://granthamdist.sourceforge.io/>). It can be used for calculation of pairwise Grantham divergence for any set of aligned MHC alleles of any species.

Allele Frequencies
------------------

Information about HLA allele frequencies in different human populations where obtained from the Allele Frequency Net Database (AFND) ([@msy116-B21]). We considered only populations of European ancestry with large sample sizes and for which frequencies of alleles at second field resolution were available: USA NMDP European Caucasian (*N* = 1,242,890), German (*N*= 39,689), and Polish (*N*= 20,653) populations. Furthermore, as with the analyses above, we focused on alleles defined as "*common*" in the CWD catalogue, which led to exclusion of some alleles with a frequency \<1%. For each population, we first determined the most common alleles (allele frequency \>= 5%) and for all the alleles under investigation in a given population, we calculated the average Grantham pairwise divergence to the most common alleles, considering all possible heterozygote genotypes.

Statistical Analyses
--------------------

### Correlation Tests

The Shapiro--Francia test was performed for all the parameters under investigation (i.e., measures of genetic distance, combined number of bound peptides and average Grantham pairwise amino acid divergence to the most common alleles) to explore samples' distribution. As parameters were not normally distributed and tied ranks could be detected within our data, the nonparametric Kendall correlation was used to test for associations between parameters. When testing the association between sequence divergence and functional divergence, all *P* values were adjusted for multiple testing using a sequential Bonferroni correction across the number of alleles tested at each locus as well as across the number of different loci tested. When testing the association between the allele's average divergence and its population frequency, *P* values were corrected across the number of populations tested. Correlations were performed in R version 3.4.1 ([@msy116-B65]).

### Permutation Tests

To test for significant differences in the strength of correlation between allele divergence and the binding to pathogen group-specific peptides, we performed permutation tests. For this analysis, the set of 232 representative human pathogen proteins were randomly shuffled among the three groups of pathogens, maintaining the same number of proteins as observed in the original data (extracellular *n* = 57, intracellular *n* = 100 and intra-extracellular *n* = 75). For each group of pathogens, permuted proteins were used to perform binding prediction analyses and compute correlation values between genetic distances and combined number of bound peptides counted for all possible allele pairs for the five HLA genes (analogous to original analysis). Each permutation was run 1,000 times, and the difference between correlation coefficients for intracellular and extracellular proteins for the five HLA genes was recorded. If there was no significant bias for intracellular or extracellular pathogens, on average this difference should be zero. The distribution of permuted differences was then used to infer the significance of our initial observations using a one-tailed test with a 0.05 cut-off.

### Artificial Proteins

Four sets of artificial proteins were created and analyzed to test for potential differentiation of the amino acid composition (AAC) among the three groups of pathogens. The first set of artificial proteins was created by randomly shuffling amino acids within each pathogen protein by using the Shuffle Protein program ([@msy116-B78]), thus maintaining the AAC of each protein intact. Three more sets of artificial proteins were created in R version 3.4.1 ([@msy116-B65]) by assembling random amino acids while maintaining several features as they occurred within each of the three pathogen groups used in the initial test (i.e., the number of proteins, the average length of sequences, the SD of the length and the minimum and maximum length). The second set of artificial proteins was created from random amino acids but maintaining the AAC as it occurred within each group of pathogen proteins. The third set of artificial proteins was created from random amino acids, while maintaining amino acid frequencies as they occur in the whole data set of pathogen proteins. Finally, amino acid composition computed from UniProtKB/Swiss-Prot data bank ([@msy116-B19]; [@msy116-B8]) was used to create the fourth set of artificial proteins.

### Multivariate Analysis of Variance

Multidimensional scaling is a multivariate statistical technique that can be used to display and summarize a high-dimensional data set in 2D graphical form. The technique was here applied to explore associations between subsets of pathogen proteins and amino acids. A nonparametric, permutational multivariate analysis of variance (PerMANOVA) was used to test for differences in the amino acid composition between pathogen groups. The PerMANOVA, based on a Bray--Curtis dissimilarity distance matrix, was run with 999 permutations to tests for statistical significance. Both procedures are implemented in the vegan package ([@msy116-B59]) in R version 3.4.1 ([@msy116-B65]).

### Comparison of Average Amino Acid Compositions

Comparison of mean amino acid compositions between the two groups of pathogen proteins (extracellular and intracellular) were performed using one-way analysis of variance; all *P* values were adjusted for multiple testing using Bonferroni correction across the number of amino acids tested.

Supplementary Material
======================

###### 

Click here for additional data file.

We thank Martin Kalbe and Jamie Winternitz for their precious input and suggestions. We also thank Matteo Fumagalli and an anonymous reviewer for very constructive comments on an earlier version of this article. This work was supported by the Max Planck Society and an Emmy Noether grant from the German Research Foundation (DFG; LE 2593/3-1 to T.L.L.).

[^1]: [Note]{.smallcaps}.---MHC allele-specific peptide binding was predicted computationally for each locus. Total number of peptides: 118,097.

[^2]: Note - Average correlation values (Kendall's tau coefficient) across the five human MHC genes (*HLA-A*, *-B*, *-C*, *-DRB1* and *-DQB1*) and rank analysis across the five parameters of sequence divergence.

[^3]: [Note]{.smallcaps}.---Correlation values (Kendall's tau) between combined number of bound peptide and Grantham genetic distance between all possible allele pairs across the five key classical MHC genes. Binding prediction was performed on the complete data set of pathogen proteins (*n* = 232) as well as considering proteins separately within three groups of pathogens: extracellular (*n* = 57), intracellular (*n* = 100), and intra-extracellular (*n* = 75).

[^4]: Tau, Kendall's tau coefficient; *P*~adj~, *P* value after Bonferroni-correction across multiple alleles tested at each locus and number of loci.
