Introduction

The Major Histocompatibility Complex (MHC) is a large locus containing many immune genes that is shared among the jawed vertebrates (Radwan et al., 2020). In humans, the MHC is also known as the HLA (Human Leukocyte Antigen) region; it spans about 5 megabases (Mb) on chromosome 6 and contains 412 genes (Figure 1A) (MHC, 2022; O’Leary et al., 2016). Many of these are part of the MHC gene family, a large group of evolutionarily-related genes with varying functions. The “classical” MHC genes are responsible for presenting protein fragments for inspection by T cells. MHC peptide presentation allows T cells to monitor the body for the presence of foreign peptides, which might indicate infection or cancer; this is crucial for vertebrate immune surveillance (Neefjes et al., 2011). “Non-classical” MHC genes are essential to the innate immune system, where they perform a variety of niche roles.

Figure 1.

The MHC locus is extraordinarily polymorphic. Haplotypes can vary widely in gene content, and thousands of distinct alleles are observed at the classical genes in humans and other primates (Maccari et al., 2017, 2020; Robinson et al., 2019). Different alleles are functionally diverse, with distinct peptide-binding affinities and, consequently, allelic differences in pathogen detection (Neefjes et al., 2011; Adams and Luoma, 2013). Given this huge diversity of functionally distinct alleles, the MHC is by far the most important locus in the genome for inter-individual variation in both infectious and autoimmune disease risk, with thousands of GWAS hits (Buniello et al., 2019; Smith et al., 2024). In our companion paper (Fortier and Pritchard, 2024), we explored the relationships between the different classical and non-classical genes. Here, we look within 17 specific genes to characterize trans-species polymorphism, a phenomemon characteristic of long-term balancing selection.

Historically, the MHC provided some of the first clear examples of positive selection in early studies of molecular evolution. By the 1980s and 1990s, researchers had noted an excess of missense variants (i.e., dN/dS > 1) in the peptide-binding regions of classical MHC genes (Hughes and Nei, 1988, 1989), alleles shared across species (Arden and Klein, 1982; Mayer et al., 1988), and high nucleotide diversity across the region (Wakeland et al., 1987; Nei and Hughes, 1991) in rodents and primates. Indeed, modern data show that nucleotide diversity in the human MHC (HLA region) exceeds 70-times the genome-wide average near the classical genes, suggesting ancient balancing selection (Figure 1A-C). Meanwhile, the MHC also features prominently in genome-wide scans for short-term directional selection (Mathieson et al., 2015; Field et al., 2016; Allentoft et al., 2022; Cong et al., 2022; Okada et al., 2018; Yasumizu et al., 2020).

In the present paper, we explore a particularly striking feature of the selection signals at MHC, namely the evidence for extremely deep coalescence structure. Some alleles (haplotypes) are more closely related to corresponding alleles from another species than they are to distinct alleles from their own species. This phenomenon is referred to as trans-species polymorphism (TSP).

TSP is rare overall in humans. Across most of the genome, human alleles coalesce to a common ancestor well within the human lineage, typically around 2 million years ago (Mallick et al., 2016). Indeed, only ∼100 loci genome-wide show compelling evidence for sharing of ancestral alleles between humans and our closest relatives, chimpanzees (Leffler et al., 2013). TSP among humans and more distantly-related species is even rarer; besides the MHC, the only other clear example of deep TSP is at the ABO locus (which influences blood type) (Azevedo et al., 2015). At this locus, both the A and B alleles are shared by descent throughout the apes, implying that the A and B lineages date back to at least the divergence-point of humans and gibbons 23 million years ago (Ségurel et al., 2012; Kuderna et al., 2023). Such deep coalescence is extraordinarily unlikely under a neutral model, and instead points to some form of balancing selection.

Meanwhile, TSP is evident at multiple MHC genes, and in many different phylogenetic clades. TSP at this locus was first proposed in the 1980s on the basis of unusual sequence similarity between mice and rats (Klein, 1980; Arden and Klein, 1982; Klein, 1987; Figueroa et al., 1988; McConnell et al., 1988; Wakeland et al., 1987), and between humans and chimpanzees (Lawlor et al., 1988; Mayer et al., 1988). Later work has reported likely TSP between humans and apes (Slierendregt et al., 1995; Boyson et al., 1996; McKenzie et al., 1999; Wroblewski et al., 2017) and humans and old world monkeys (Mayer et al., 1992; Brändle et al., 1992; Kupfermann et al., 1992; Slierendregt et al., 1992; Geluk et al., 1993; Satta et al., 1996; Otting et al., 2000; Kriener et al., 2000; Gyllensten et al., 1990; Kriener et al., 2001; Otting et al., 1992; Otting and Bontrop, 1995; Otting et al., 1992, 2002); deep TSP is also consistent with the high levels of genetic diversity within the MHC. Such ancient TSP would make the MHC unique compared to any other locus in the genome. However, most previous work has not fully accounted for the inherent uncertainty in phylogenetic inference, especially given the potential for convergent evolution at functional sites. Although there is clear evidence for TSP, its exact age at each gene is still uncertain.

To address these questions, we used data from the IPD-MHC/HLA database—a large repository for MHC allele sequences from humans, non-human primates, and other vertebrates—along with supplementary sequences from NCBI RefSeq (Maccari et al., 2017, 2020; Robinson et al., 2019). This represents the most complete sampling of primate MHC genes to date, spanning the entire primate tree (Figure 1D). We account for the uncertainty in phylogenetic inference using a Bayesian MCMC approach (BEAST2) which is well-suited to handle highly variable and rapidly-evolving sequences. In our companion paper (Fortier and Pritchard, 2024), we compared genes across dozens of species and learned more about the orthologous relationships among them. That process helped inform this work, where we assess support for TSP within individual genes.

We find support for TSP among the African apes for genes MHC-C, -DPA1, and -DRB3, among the great apes for MHC-DPB1, and among all apes for MHC-B. We also find conclusive evidence for TSP at least back to the ancestor of humans and OWM in MHC-DQB1, implying—remarkably—that alleles have been maintained by balancing selection for at least 31 million years. Rapidly-evolving sites are mainly located in the critical peptide-binding regions of the classical genes, but are spread throughout the coding region of the non-classical genes. Moreover, the most rapidly-evolving sites are also frequently associated with immune phenotypes and diseases in the literature, connecting our evolutionary findings with their functional consequences. These results highlight the contrasting roles of ancient balancing selection and short-term directional selection within the peptide-binding regions of the classical genes, and motivate further evolutionary and functional studies to better understand this unique system.

Results

Data

We collected MHC nucleotide sequences for all genes from the IPD-MHC/HLA database, a large repository for MHC alleles from humans, non-human primates, and other vertebrates (Maccari et al., 2017, 2020; Robinson et al., 2024). Although extensive, this database includes few or no sequences from important primates such as the gibbon, tarsier, and lemur. Thus, we supplemented our set of alleles using sequences from NCBI RefSeq (O’Leary et al., 2016). Because the MHC genes make up an evolutionarily-related family, they can all be aligned (Kaufman, 2022; Adams and Luoma, 2013). In our companion paper (Fortier and Pritchard, 2024), we utilized these large multigene alignments for Class I, Class IIA, and Class IIB to compare genes. Here, we analyze subsets of those alignments, each focusing on a gene or group of closely-related genes.

We considered 16 gene groups and two or three different genic regions for each group: exon 2 alone, exon 3 alone, and/or exon 4 alone. Exons 2 and 3 encode the peptide-binding region (PBR) for the Class I proteins, and exon 2 alone encodes the PBR for the Class II proteins. For the Class I genes, we also considered exon 4 alone because it is comparable in size to exons 2 and 3 and provides a good contrast to the PBR-encoding exons. See the Methods for more detail on how gene groups were defined. Because few intron sequences were available for non-human species, we did not include them in our analyses.

Trans-Species Polymorphism is Widespread

For each gene group and genic region, we used the Bayesian phylogenetics software BEAST2 (Bouckaert et al., 2014, 2019) with package SubstBMA (Wu et al., 2013) to infer phylogenies. One major advantage of BEAST2 over less tuneable methods is that it can allow evolutionary rates to vary across sites, which is important for genes such as these which experience rapid evolution in functional regions (Wu et al., 2013). We also considered each exon separately to minimize the impact of recombination as well as to compare and contrast the binding-site-encoding exons with non-binding-site-encoding exons.

We can visualize each set of phylogenies as a single summary tree, which maximizes the sum of posterior clade probabilities (BEA, 2024). Three of these summary trees are shown in Figure 2, constructed from the second exons of classical Class I gene MHC-C, classical Class II gene MHC-DQB, and non-classical Class II gene MHC-DOA, respectively (see Figure 2—figure Supplement 1-Figure 2—figure Supplement 10 for the other exons and genes).

Figure 2.

In Figure 2, each tip represents a sequence, labeled with its allele name (see Appendix 1 for more details on nomenclature) or its RefSeq accession number. Groups of tips containing similar alleles are collapsed for clarity. Preceding the allele label is a colored rectangle with a four-letter abbreviation indicating the species; the species key is in the middle of the plot and also displays the known species tree (Kuderna et al., 2023).

MHC-C is a classical Class I gene which duplicated from MHC-B in the ancestor of the great apes (Piontkivska and Nei, 2003; Fukami-Kobayashi et al., 2005; Abi-Rached et al., 2010; Adams and Parham, 2001; Lugo and Cadavid, 2015). Its protein product participates in classical antigen presentation and also serves as the dominant Class I molecule for interacting with killer cell immunoglobulin-like receptors (KIRs) in innate immunity (Adams and Parham, 2001; Guethlein et al., 2015; Vollmers et al., 2021). MHC-DQB is a classical Class II gene which pairs with MHC-DQA. Apes have two MHC-DQ copies, MHC-DQA1/MHC-DQB1 and MHC-DQA2/MHC-DQB2, while the second copy was deleted in OWM. NWM can have two or three sets of MHC-DQ genes, depending on species, but it has been unclear whether any of them are 1:1 orthologous with the ape or OWM genes. Lastly, MHC-DOA is a non-classical Class II gene whose protein product modulates MHC-DM activity, indirectly affecting Class II peptide presentation (Heijmans et al., 2020; Neefjes et al., 2011). The genes’ differing roles result in different patterns in the phylogenetic trees.

Critically, we observe that, at classical genes MHC-C (Figure 2A) and MHC-DQB (Figure 2B), the alleles fail to cluster together according to species, as indicated by the mixed-color clades throughout the trees. In MHC-C, human (HLA; red), chimpanzee (Patr; dark pink), bonobo (Papa; light pink), and even gorilla (Gogo; orange) alleles can be found throughout the tree, indicating that variation in this gene is almost as old as the gene itself.

Figure 2B displays the BEAST2 tree consisting of MHC-DQB1, -DQB2, and outgroup -DQB alleles all together. It shows many mixed-color clades throughout, consisting of ape (red/orange/yellow) and OWM (green) alleles grouping together. Alleles often group by first-field name instead of by species, indicating that some allelic lineages have been maintained since before the split of humans and OWM—at least 31 million years. An example of this is shown in the inset to the left of this tree, “Example of Human-OWM TSP.” Here, human alleles coalesce with OWM alleles before they coalesce with each other. Near the bottom of the tree is a clade consisting of ape and NWM MHC-DQB2 sequences, suggesting that they are orthologous; however, NWM species have expanded their MHC-DQ regions so these genes may not be 1:1 orthologous (see our companion paper, Fortier and Pritchard (2024)). Additionally, Strepsirrhini sequences do not group with either the MHC-DQB1 or -DQB2 clade, showing that duplications of these genes must have happened in or after the Simiiformes ancestor.

Figure 2C shows the BEAST2 summary tree for non-classical MHC-DOA. In this tree, alleles group exclusively by species (clades are collapsed for clarity) and the branching order of the species deviates only slightly from the species tree. This shows that not all MHC genes are affected by long-term balancing selection, despite the complicated linkage disequilibrium across classical and non-classical genes in the region (Smith et al., 2024). This also suggests that antigen presentation specifically, as opposed to a general role in immune function, is the driving force behind this long-term balancing selection.

While the BEAST2 summary trees in Figure 2 are suggestive of deep TSP, they do not directly quantify the statistical confidence in the TSP model. Moreover, standard approaches to quantifying uncertainty in trees, such as bootstrap support or posterior probabilities for specific clades, do not relate directly to hypothesis testing for TSP. We therefore implemented an alternative approach using BEAST2 output, as follows (see the Methods for details).

We performed formal model testing for TSP within quartets of alleles, where two alleles are taken from a species (or taxon) A, and two alleles are taken from a different species (or taxon) B. If the alleles from A group together (and the alleles from B group together) in the unrooted tree, this quartet supports monophyly of A (and of B). But if alleles from A group more closely with B in the unrooted tree, then this comparison supports TSP. Since BEAST2 samples from the posterior distribution of trees, we counted the number of trees that support TSP versus the number that support monophyly as an estimate of the posterior support for each model. We then summarized the relative support for each model by converting these to Bayes factors (see Methods for more detail). The precise interpretation of Bayes factors depends on one’s prior expectation; however, following standard guidelines (Jeffreys, 1998) we suggest that Bayes factors > 100 should be considered as strong support in favor of TSP. Bayes factors <1 are evidence against TSP. For each comparison of two taxa, we report the maximum Bayes factor across the possible quartets, as we are interested in whether any quartet shows compelling evidence for TSP.

Gene conversion, the unidirectional transfer of short tracts of DNA from a donor to an acceptor sequence, can affect the inferred trees. In particular, acceptor sequences may group more strongly with donor sequences than with sequences that share DNA by descent. This can make it difficult to distinguish trees influenced by trans-species polymorphism from those influenced by gene conversion. Thus, we inferred gene conversion tracts using GENECONV (Sawyer, 1999) and excluded significant gene-converted acceptor alleles from the Bayes factor calculations. While GENECONV cannot possibly infer all past events, this procedure should ameliorate any biasing effects. Additionally, we consider each exon separately; analyzing short tracts reduces the effect of recombination on the tree (see our companion paper for more specifics (Fortier and Pritchard, 2024)).

Bayes factors are shown in Figure 3 and Figure 4. For each genic region (x-axis), we tested for TSP among humans, non-human primates, and non-primates with each of the comparisons listed on the y-axis. Each table entry is colored and labeled with the maximum Bayes factor among all tested quartets of alleles belonging to that category. See Figure 2—figure Supplement 1-Figure 2—figure Supplement 16, Figure 3—figure Supplement 1-Figure 3—figure Supplement 12, and Figure 4—figure Supplement 1-Figure 4—figure Supplement 10 for examples of high-Bayes-factor quartets for each comparison. In Figure 3 and Figure 4, orange indicates support for TSP among the species in that category, while teal indicates evidence against TSP.

Figure 3.

Figure 4.

At the Class I genes (Figure 3), MHC-C shows strong support for TSP within the African apes: human, chimpanzee, and gorilla. Having arisen fairly recently in the ancestor of human and orangutan, MHC-C has thus maintained some allelic lineages for most of its history. TSP has not previously been reported for this gene.

For MHC-A, Bayes factors vary considerably depending on exon and species pair. Many sequences had to be excluded from MHC-A comparisons because they were identified as gene-converted in the GENECONV analysis or were previously identified as recombinants (Hans et al., 2017; Gleimer et al., 2011; Adams and Parham, 2001). Importantly, for MHC-A we do not see concordance in Bayes factors across the different exons, whereas we do for the other gene groups. Similar Bayes factors across all exons for a given comparison is thus evidence in favor of TSP being the primary driver of the observed deep coalescence structure (rather than recombination or gene conversion).

The non-classical Class I genes MHC-E, -F, and -G (bottom row of Figure 3) are interspersed with the classical Class I genes in the MHC region, but their products have niche functions in innate immunity. Their indirect involvement in adaptive immunity means they experience different selective pressures. They exhibit lower polymorphism and dN/dS < 1, reflecting the fact that they have not been subject to the same pathogen-mediated balancing selection. The Bayes factors for all three of these genes show strong evidence against TSP, as expected. Interestingly, despite its non-classical role, MHC-E has a known balanced polymorphism in humans; the two main alleles are at similar frequencies worldwide but may have different expression levels and peptide preferences (Paganini et al., 2019; Grant et al., 2020). Our approach—meant to detect ancient TSP—does not reveal balancing selection in MHC-E, showing that this balanced polymorphism is young. For MHC-G, there were not enough sequences available to perform many of the tests (at least two from each species group are required). However, we do not expect to see evidence of TSP in this gene.

Each Class II MHC molecule has an α and β component which are encoded by an A and B gene, respectively. Bayes factors for the Class IIA genes are shown in the top row of Figure 4, while those for their Class IIB partners are shown in the bottom row. The non-classical MHC-DM and -DO molecules assist the classical Class II genes with peptide loading and are not believed to be shaped by balancing selection (Heijmans et al., 2020; Neefjes et al., 2011). As expected, we see strong evidence against TSP between humans and all other primate species for these genes (first two columns of Figure 4).

In contrast, we find evidence for deep TSP within the classical Class II genes. MHC-DPA1 shows TSP between human, chimpanzee, and gorilla in exon 2, but not in exon 3. We find that this TSP in MHC-DPA1 is less deep than has been previously suggested with non-Bayesian methods (Otting and Bontrop, 1995), underscoring the importance of this methodology for handling the MHC. Meanwhile, its partner MHC-DPB1 shows strong evidence for TSP between human and orangutan in exon 3 and suggestive evidence in exon 2; our work provides the first evidence of TSP between humans and other apes for this gene (Slierendregt et al., 1995).

The MHC-DR genes behave somewhat differently than the other classical Class II molecules. While the α and β components of all the other molecules engage in exclusive binding, there are many different MHC-DRβ molecules which all bind to the same MHC-DRα. The MHC-DRA gene is conserved across species with little polymorphism, while the MHC-DRB region is highly variable both in gene content and allelic diversity. Consistent with this, Bayes factors for the MHC-DRA gene reveal strong evidence against TSP for all species pairs, while MHC-DRB1 shows strong evidence in favor of TSP between human, chimpanzee, gorilla, OWM, and even NWM in exon 2. However, because the Bayes factors only support TSP between humans and OWM/NWM in MHC-DRB1 in exon 2, but not in exon 3, this could mean alleles are not actually that ancient. We show in our companion paper that individual MHC-DRB genes are short-lived, and only three are truly orthologous between apes and OWM (Fortier and Pritchard, 2024). These pieces of evidence suggest previous work may have overstated the extent of TSP at this locus.

Remarkably, the MHC-DQB1 gene shows definitive evidence for TSP back to at least the ancestor of humans and OWM. While this result has been presented previously, we confirm it with decisive evidence (Bayes factor >100) for allelic lineages being maintained for over 31 million years (Otting et al., 1992, 2002; Loisel et al., 2006; Simons et al., 2017).

Many of our BEAST trees also showed intermingling of sequences from different species within the OWM or NWM (e.g. Figure 2—figure Supplement 4), regardless of whether they formed transspecies clades with ape sequences. This could indicate trans-species polymorphism within the OWM or NWM that is still ancient, but not old enough to be shared with the apes. Therefore, we also calculated Bayes factors across different clades of OWM and NWM (Figure 3—figure Supplement 13, Figure 3—figure Supplement 14, Figure 4—figure Supplement 11, and Figure 4—figure Supplement 12). We see very strong evidence for TSP between all groups of OWM for MHC-DPA1, -DPB1, -DQA1, and -DQB1, indicating that allelic lineages at these genes have been maintained within the OWM for at least 19 million years (MY). Unexpectedly, we also see evidence for TSP in some non-classical genes. MHC-E, a gene which is non-classical in humans and is presumed non-classical in OWM (yet is duplicated in some species), shows evidence for 15-MY-old TSP within the OWM. In non-classical MHC-DMB, we also observe TSP within the OWM as old as 11 MY. Due to the uncertainty of locus assignments for alleles of the OWM MHC-A, -B, and -DRB genes and of the NWM genes, we cannot make definitive conclusions about TSP within these groups for these other genes.

In summary, the phylogenetic analyses point to ancient TSP in classical genes MHC-B, -C, -DPA1, -DPB1, -DQB1, -DRB1, and -DRB3. Bayes factors for the non-classical genes MHC-E, -F, -G, -DMA, -DMB, -DOA, and -DOB suggest there is no TSP involving apes at these loci—as expected for non-classical genes—confirming the usefulness of our approach. However, there may still be TSP at these genes within other groups, such as the OWM, hinting at possible functional differences. Overall, TSP is more ancient among the Class II genes than the Class I genes, consistent with their older age.

From Evolution to Function

Alongside the evidence for ancient TSP, the MHC region is also notable for its high rate of missense substitutions (dN/dS > 1) (Hughes and Nei, 1988, 1989) and its large number of GWAS hits for autoimmune and infectious diseases (Buniello et al., 2019; Kennedy et al., 2017). We next aimed to understand how these observations relate to signals of TSP and known features of the MHC proteins.

To explore these questions, we first estimated the per-site evolutionary rates within each gene. As in our TSP analysis, we used the BEAST2 package SubstBMA, which estimates evolutionary rates at every site concurrently with a tree. We averaged these rates over all states in the chain to get per-site evolutionary rates for each site, then calculated their fold-change relative to the average rate among mostly-gap sites in the alignment (“baseline”; see Methods).

Figure 5A shows the substitution rate fold-change for each nucleotide along the concatenated coding sequence of Class I genes MHC-B, -C, and -E. Exon boundaries are indicated by orange vertical dashed lines. In classical genes MHC-B and -C, nearly all the rapidly-evolving sites lie within exons 2 and 3, which encode the protein’s peptide-binding domain. While exons 2 and 3 make up only ∼50% of the gene’s length, they contain 94% and 90% of the sites evolving at more than 4 times the baseline evolutionary rate for the classical genes MHC-B and MHC-C, respectively. In MHC-B, exons 2 and 3 each show significantly higher proportions of rapidly-evolving sites compared to the “other” exons (exons in the gene excluding 2, 3, or 4), while the difference is not significant for MHC-C (Figure 5—figure Supplement 2). This result could reflect the relatively young age of MHC-C or its additional role as the dominant Class I molecule for interacting with KIRs (Adams and Parham, 2001; Guethlein et al., 2015; Vollmers et al., 2021; Piontkivska and Nei, 2003).

Figure 5.

In contrast to these classical genes, non-classical MHC-E primarily presents self-peptides for recognition by NK cell receptors, and its peptide-binding groove is tailored to accommodate a very specific set of self-peptides—leader peptides cleaved from other Class I MHC proteins during processing (Miller et al., 2003). As shown in Figure 5A, this gene shows fewer rapidly-evolving sites than the classical genes. These sites are also relatively evenly distributed across the gene, with exons 2 and 3 (which cover ∼50% of the gene’s length) containing 45% of the sites evolving at over 4 times the baseline evolutionary rate. Interestingly, exon 4—a non-peptide-binding exon of equal size—displays a significantly lower proportion of rapidly-evolving sites compared with the “other” exons (Figure 5—figure Supplement 2). These results support that MHC-E has been remarkably conserved across the primates and that its evolution may not be driven by differential peptide binding (Heijmans et al., 2020).

We then examined where the rapidly-evolving sites lie within the physical protein structures. To do this, we averaged the per-site rates within each codon to get per-amino-acid rates, then mapped these onto the known human protein structures. Unfortunately, there are few non-human primate protein structures in the Protein Data Bank (Berman et al., 2000), but the macaque structures we found were nearly identical to those of human. Figure 5B shows structures for human HLA-B, -C, and -E; this view features the peptide (black) sitting in the peptide-binding groove (flanked on top and bottom by helices). In MHC-B and -C, rapidly-evolving amino acids (orange) tend to be located within the peptide-binding groove. To quantify this, we measured the minimum distance between each amino acid and the bound peptide. We averaged these distances over several structures, which are listed in the Methods. For all three proteins, amino acids closer to the peptide have significantly higher evolutionary rates than amino acids further from the peptide, as shown in Figure 5C (see also Figure 5—figure Supplement 4 and Figure 5—figure Supplement 5). The effect is much less pronounced in MHC-E, where even the amino acids closest to the peptide do not exhibit high evolutionary rates. These results are consistent with the expectation that rapid evolution and diversity at the classical MHC genes would be mediated by selective pressures for changes in peptide binding.

The rapidly-evolving sites for the Class II genes are shown in Figure 6. Panel A shows the substitution rate fold-change for each nucleotide along the concatenated coding sequence of the Class II MHC-DRA, -DQA, -DRB, and -DQB gene groups. Both MHC-DQA and -DQB are extraordinarily polymorphic, but MHC-DRA is conserved compared to its multiple, highly-variable MHC-DRB partners. Indeed, Figure 6A shows that rapidly-evolving sites are concentrated in binding-site-encoding exon 2 for MHC-DRB, -DQB, and to a lesser extent MHC-DQA. Exon 2, which makes up ∼30% of the coding region, contains 32% of sites evolving at more than twice the baseline rate in MHC-DRA, but 57% of such sites in MHC-DQA, 61% in MHC-DQB, and 73% in MHC-DRB. Comparing across exons, exon 2 contains a significantly higher proportion of rapidly-evolving sites compared to the “other” exons in classical MHC-DQA, -DQB, -DRB, and -DPB, but also—curiously—in non-classical MHC-DMA and -DMB. It is interesting that MHC-DM appears to be evolving rapidly in its binding-site-encoding exons, despite the fact that it is not thought to bind peptides. Instead, it is responsible for assisting with peptide loading onto the classical genes. Co-evolution with MHC-DR in particular seems possible; the interaction between the MHC-DM and -DR molecules depends on the affinity of the peptide trying to bind with MHC-DR. MHC-DM thus shapes the repertoire of peptides presented by MHCDR, favoring high-affinity peptides (Dijkstra and Yamaguchi, 2019; Schulze and Wucherpfennig, 2012). It is plausible that the host-pathogen evolution shaping the MHC-DRB genes has resulted in co-evolution of MHC-DMA and -DMB to maintain this regulatory interaction.

Figure 6.

We again mapped the evolutionary rates onto human protein structures, shown in Figure 6B. In each molecule, the α chain is positioned at the top and encompasses the upper helix forming the binding site, while the β chain is oriented toward the bottom and encompasses the lower helix. The peptide is shown in black. Rapidly-evolving sites are concentrated in each protein’s binding site, although in MHC-DR this is more prominent in the bottom helix (MHC-DRB).

We then measured the distance between each amino acid and the bound peptide, shown in Figure 6C. MHC-DRA did not show a significant relationship between evolutionary rate and distance, as expected by its relatively uniform distribution of evolutionary rates. For the other three proteins, amino acids closer to the peptide had significantly higher evolutionary rates. This held true for the classical MHC-DPA and -DPB genes as well (Figure 5—figure Supplement 4). Again, this is consistent with differential peptide binding and TCR responsiveness driving the diversity and longterm balancing selection at the classical genes. We could not measure distances for non-classical MHC-DOA, -DOB, -DMA, and -DMB genes because they do not engage in peptide presentation.

Lastly, since the rapidly-evolving sites are likely involved in peptide binding, they also influence the response to pathogens and self-antigens, presumably affecting risk for infectious and autoimmune diseases. To bridge the gap between evolution and complex traits, we collected HLA fine-mapping studies for infectious, autoimmune, and other diseases, as well as for biomarkers and TCR phenotypes. These studies report associations between a disease or trait and classical HLA alleles, SNPs, and amino acid variants, often with multiple independent hits per gene. They demonstrate that HLA variation affects disease in complex ways—sometimes, a single variant is strongly associated with a condition, while other times, a combination of amino acids or even entire alleles are the strongest indicators of disease susceptibility or protection.

We were interested in whether our rapidly-evolving amino acids from the BEAST2 analysis corresponded with disease-associated amino acids from the literature. Table 1 lists disease, trait, and TCR-phenotype associations for the most rapidly-evolving amino acids (fold-change g 1) of the MHC-B group (see Table 1—source data 1 for the other genes). The majority of rapidly-evolving positions in MHC-B have at least one association.

Table 1.

Thus, in summary, we find that rapid evolution has primarily targeted amino acids within the peptide-binding region of each gene, and that these specific positions are likely the primary drivers of phenotypic associations at the MHC locus.

Discussion

The MHC region contains the clearest signals of balancing and directional selection in mammalian genomes, including extreme diversity, ancient trans-species polymorphism, and high rates of non-synonymous evolution between allelic lineages. In humans, MHC/HLA variation is associated with risk for infectious and autoimmune diseases and many other traits, and HLA matching is critical for successful tissue transplantation (Kennedy et al., 2017; Smith et al., 2024; Lee et al., 2007).

Despite its evolutionary and clinical importance, the extreme diversity of the MHC makes it challenging to study, and basic questions about its evolutionary history remain unresolved. Although past work has suggested ultra-deep TSP at this locus, in this study we re-examined the region with modern, comprehensive data and a unified analysis framework. Using Bayesian evolutionary analysis, we report conclusive evidence for long-term TSP in seven classical genes, including between humans and OWM at MHC-DQB1. Thus, remarkably, allelic lineages at this gene have been maintained for at least 31 million years.

Our evidence for TSP at MHC-DQB1 spanning at least 31 million years places it among the most ancient examples of balancing selection known in any species, and almost certainly the oldest in primates. Aside from MHC, the deepest example within primates is at the ABO locus controlling blood type; it exhibits trans-species polymorphism between humans and gibbons, an age of 23 million years (Ségurel et al., 2012). In various chimpanzee species, OAS1, which helps inhibit viral replication, contains alleles up to 13MY old (Ferguson et al., 2012). TSP between chimpanzee and human includes LAD1, a protein that maintains cell cohesion (6MY) (Teixeira et al., 2015), retroviral transcription factor TRIM5α (4-7MY in apes and > 8MY in OWM) (Cagliani et al., 2010; Newman et al., 2006), and ZC3HAV1, an antiviral protein leading to viral RNA degradation (6MY) (Cagliani et al., 2012), among others (Leffler et al., 2013).

Looking more broadly across the tree of life, ancient trans-species polymorphism occurs widely, albeit rarely. Several of the oldest examples are found in the MHC locus: MHC polymorphisms have been maintained for 35MY in cetaceans (Xu et al., 2009), 40MY in herons (Li et al., 2011), 48MY in mole rats (Kundu and Faulkes, 2007), 70MY in tree frogs (Zhao et al., 2013), and over 105MY in salmonid fishes (Kiryu et al., 2005; Grimholt et al., 2015). There are examples in non-MHC loci as well; in cyanobacteria, polymorphism at the HEP island controlling heterocyst function has been maintained for 74 million years (Sano et al., 2018), in plants, S-genes determining self-incompatibility exhibit TSP spanning 36 million years (Ioerger et al., 1990; Igic et al., 2006), and in Formica ants, alleles at a supergene underlying colony queen number have been maintained for over 30 million years (Purcell et al., 2021).

Paradoxically, given the extremely long-lived balancing selection acting in these lineages, many authors have also reported strong directional selection (Brandt et al., 2018; Nunes et al., 2021; Bhatia et al., 2011). Indeed, within the phylogeny we find that the most rapidly-evolving codons are substituted at around 2–4-fold the baseline rate. For the classical genes (except MHC-DRA), these rapidly-evolving sites lie within the peptide binding regions of the corresponding proteins, usually very close to the peptide-contact surfaces. This does not hold true for the non-classical genes, supporting the fact that selection at the MHC is mainly driven by the peptide presentation pathway. Non-classical MHC-DMA and -DMB are a surprising exception, showing significantly elevated proportions of rapidly-evolving sites in exon 2 similarly to the classical genes (even though the MHC-DM molecule does not bind peptides). This pattern may be caused by its co-evolution with the classical Class II genes, and more research is needed to address this.

The primary role of classical MHC proteins is to present peptides for T cell recognition; we found that the same rapidly-evolving amino acids are consistently associated with shaping T cell receptor (TCR) repertoires. Moreover, these amino acids are frequently associated with autoimmune and infectious diseases in HLA fine-mapping studies.

Taken together, we begin to see a comprehensive picture of the nature of primate MHC evolution. In response to rapidly-changing pathogen pressures, the PBRs of classical MHC proteins evolve to bind changing pathogen antigens and present them to TCRs. Broad lineages of MHC alleles are maintained over tens of millions of years by strong balancing selection, providing defense against a wide variety of different pathogens. Yet within these lineages, alleles turn over quickly in response to new specific threats. This reconciles evidence for TSP, the presence of thousands of alleles, and the existence of rapidly-evolving sites.

MHC molecules must be trained to detect pathogens with both specificity and sensitivity, and distinguishing self from non-self peptides is a challenging task. As MHC proteins evolve, there is an unavoidable flux between infection defense and autoimmune susceptibility. Additionally, many MHC proteins have roles in both innate and adaptive immunity. As a result, rapidly-evolving amino acids are associated with both infections and autoimmune conditions. In the future, disease studies within other primate species could provide insight into the trajectory of MHC evolution and might reveal evolutionary trade-offs. Perhaps balancing selection has kept the same amino acids disease-relevant across the entire primate evolutionary tree, or maybe the rapid turnover of MHC variation means different primate clades will have different disease associations.

Although the primate MHC has been of interest to evolutionary biologists for more than 30 years, there is still much to be done to more fully document the evolution of the MHC genes within and between species. Moreover, we still have a limited understanding of how sequence changes map to functional differences among alleles, and how these relate to allele-specific profiles of pathogen protection (and autoimmunity risk). However, functional and computational advances will provide key opportunities for progress on these problems (Radwan et al., 2020; Vizcaíno et al., 2020).

Methods and Materials

Data

We downloaded MHC allele nucleotide sequences for all human and nonhuman genes from the IPD Database (updated January 2023) (Barker et al., 2023; Maccari et al., 2017, 2020). To supplement the alleles available in the database, we also collected nucleotide sequences from NCBI using the Entrez E-utilities with query “histocompatibility AND txidX AND alive[prop]”, where X is a taxon of interest.

We wanted to provide “zoomed-in” versions of various subtrees within the multi-gene trees presented in our companion paper (Fortier and Pritchard, 2024). Thus, we included more species and more alleles per species than in the original trees. In each tree, we also included a “backbone” of sequences from the overall multi-gene tree to provide context for each expanded clade (lists of alleles provided as Supplementary Files).

For Class I, we expanded the following clades: 1) MHC-A-related genes (MHC-A group), 2) MHC-B-related genes (MHC-B group), 3) MHC-C-related genes (MHC-C group), 4) MHC-E-related genes (MHC-E group), 5) MHC-F-related genes (MHC-F group), and 6) MHC-G-related genes (MHC-G group). For Class IIA, we expanded: 1) MHC-DMA-related genes (MHC-DMA group), 2) MHC-DOA-related genes (MHC-DOA group), 3) MHC-DRA-related genes (MHC-DRA group), 4) MHC-DPA-related genes (MHC-DPA group), and 5) MHC-DQA-related genes (MHC-DQA group). For Class IIB, we expanded: 1) MHC-DMB-related genes (MHC-DMB group), 2) MHC-DOB-related genes (MHC-DOB group), 3) MHC-DRB-related genes (MHC-DRB group), 4) MHC-DPB-related genes (MHC-DPB group), and 5) MHC-DQB-related genes (MHC-DQB group). These sets were inclusive of all orthologs and paralogs of a given human gene across all species we included (see our companion paper for more information (Fortier and Pritchard, 2024)). For example, the MHC-A group includes human HLA-A and its 1:1 orthologs in the apes, the expanded MHC-A and -AG paralogs of the OWM, chimpanzee specific Patr-AL, gorilla-specific Gogo/Gobe-OKO, orangutan-specific Poab/Popy-Ap, and pseudo-genes MHC-H and -Y. See Supplementary Files for all alleles included within each group.

We aligned each group separately using MUSCLE (Edgar, 2004) with default settings and manually adjusted, following the alignments we already produced for the multi-gene trees in our companion paper (Fortier and Pritchard, 2024).

Nucleotide Diversity

The classical MHC region is defined as chr6:28,510,120-33,480,577 (GRCh38) (MHC, 2022). Nucleotide diversity (π) was calculated on the modern human data using VCFtools (0.1.15) (Danecek et al., 2011). For the entire MHC region (Figure 1A), π was calculated in 5000bp sliding windows with a step size of 1000bp. For each gene separately (Figure 1B-C and Figure 1—figure Supplement 1), π was calculated in 50bp sliding windows with a step size of 10bp.

Bayesian Phylogenetic Analysis

We constructed phylogenetic trees using BEAST2 (Bouckaert et al., 2014, 2019) with package SubstBMA (Wu et al., 2013). SubstBMA implements a spike-and-slab mixture model that simultaneously estimates the phylogenetic tree, the number of site partitions, the assignment of sites to partitions, the nucleotide substitution model, and a rate multiplier for each partition. Since we were chiefly interested in the partitions and their rate multipliers, we used the RDPM model as described by Wu et al. (2013). In the RDPM model, the number of nucleotide substitution model categories is fixed to 1, so that all sites, regardless of rate partition, share the same estimated nucleotide substitution model. This reduces the number of parameters to be estimated and ensures that only evolutionary rates vary across site partitions, reducing overall model complexity. We used an uncorrelated lognormal relaxed molecular clock because in reality, evolutionary rates vary among branches (Bergeron et al., 2023).

Priors

For the Dirichlet process priors, we used the informative priors constructed by Wu et al. (2013) for their mammal dataset. This is appropriate because they include several of the same species and their mammals span approximately the same evolutionary time that we consider in our study. We also use their same priors on tree height, base rate distribution, and a Yule process coalescent prior. We did not specify a calibration point—a time-based prior on a node—because we did not expect our sequences to group according to the species tree.

Running BEAST2

We aligned sequences across genes and species and ran BEAST2 on various subsets of the alignment. For the Class I gene groups (MHC-A group, MHC-B group, MHC-C group, MHC-E group, MHC-F group, and MHC-G group), we repeated the analysis for 1) exon 2 only (PBR), 2) exon 3 only (PBR), 3) exon 4 only (non-PBR), and 4) exons 1, 5, 6, 7, and 8 together (non-PBR; “other” exons). For the Class IIA gene groups (MHC-DMA group, MHC-DOA group, MHC-DRA group, MHC-DPA group, and MHC-DQA group), we used 1) exon 2 only (PBR), 2) exon 3 only (non-PBR), and 3) exons 1, 3, 4, and 5 together (non-PBR; “other” exons). For Class IIB gene groups (MHC-DMB group, MHC-DOB group, MHC-DRB group, MHC-DPB group, and MHC-DQB group), we analyzed 1) exon 2 only (PBR), 2) exon 3 only (non-PBR), and 3) exons 1, 3, 4, and 5 together (non-PBR; “other” exons). In the following, each “analysis” is a collection of BEAST2 runs using one of these sets of exons of a particular gene group.

The xml files we used to run BEAST2 were based closely on those used for the mammal dataset with the RDPM model and uncorrelated relaxed clock in Wu et al. (2013)(https://github.com/jessiewu/substBMA/blob/master/examples/mammal/mammal_rdpm_uc.xml). Running a model with per-site evolutionary rate categories and a relaxed clock means there are many parameters to estimate. Along with the large number of parameters, the highly-polymorphic and often highly-diverged sequences in our alignments makes it difficult for BEAST2 to explore the state space. Thus, we under-took considerable effort to ensure good mixing and convergence of the chains. First, we employed coupled MCMC for all analyses. Coupled MCMC is essentially the same as the regular MCMC used in BEAST2, except that it uses additional “heated” chains with increased acceptance probabilities that can traverse unfavorable intermediate states and allow the main chain to move away from an inferior local optimum (Müller and Bouckaert, 2020). Using coupled MCMC both speeds up BEAST2 runs and improves mixing and convergence. We used four heated chains for each run with a delta temperature of 0.025. Second, we ran each BEAST2 run for 40,000,000 states, discarding the first 4,000,000 states as burn-in and sampling every 10,000 states. Third, we ran at least 8 independent replicates of each analysis. The replicates use the exact same alignment, but explore state space independently and thus are useful for improving the effective sample size of tricky parameters. As recommended by BEAST2, we examined all replicates in Tracer version 1.7.2 (Rambaut et al., 2018) to ensure that they were sampling from the same parameter distributions and had reached con-vergence. We excluded replicates for which this was not true, as these chains were probably stuck in suboptimal state space. Additionally, Tracer provides estimates of the effective sample size (ESS) for the combined set of states from all chosen replicates, and we required that the ESS be larger than 100 for all parameters. If there were fewer than 4 acceptable replicates or if the ESS was below 100 for any parameter, we re-ran more independent replicates of the analysis until these requirements were satisfied. We obtained between 5 and 18 acceptable replicates per analysis (median 8).

For some analyses, computational limitations prevented BEAST2 from being able to reach 40,000,000 states. In these situations, more replicates (of fewer states) were usually required to achieve good mixing and convergence. The first 4,000,000 states from each run were still discarded as burn-in even though this represented more than 10% of states in these cases.

This stringent procedure ensured that all of the replicates were exploring the same parameter space and were converging upon the same global optimum, allowing the g 4 independent runs to be justifiably combined. We combined the acceptable replicates using LogCombiner version 2.6.7 (Drummond and Rambaut, 2007), which aggregates the results across all states. We then used the combined results to perform downstream analyses.

Phylogenetic Trees

After combining acceptable replicates, we obtained 12,382 - 64,818 phylogenies per group/gene region (mean 34,499). We used TreeAnnotator version 2.6.3 (Drummond and Rambaut, 2007) to summarize each set of possible trees as a maximum clade credibility tree, which is the tree that maximizes the sum of posterior clade probabilities. Since BEAST2 samples trees from the posterior, one could in principle perform model testing directly from the posterior samples; the complete set of trees can typically be reduced to a smaller 95% credible set of trees representing the “true” tree (BEA, 2024). However, given the high complexity of the model space, all our posterior trees were unique, meaning this was not possible in practice. (Since the prior over tree topologies is unstructured, this effectively puts minuscule prior weight on trees with monophyly. Thus, sampling directly from the posterior provides an unacceptably high-variance estimator.)

Gene Conversion

We calculated gene conversion fragments using GENECONV version 1.81a (Sawyer, 1999) on each alignment. It is generally advisable to use only synonymous sites when running the program on a protein-coding alignment, since silent sites within the same codon position are likely to be correlated. However, the extreme polymorphism in these MHC genes meant there were too few silent sites to use in the analysis. Thus, we considered all sites but caution that this could slightly overestimate the lengths of our inferred conversion tracts. For each alignment, we ran GENECONV with options ListPairs, Allouter, Numsims = 10000, and Startseed = 310. We collected all inferred “Global Inner” (GI) fragments with sim_pval < 0.05 (this is pre-corrected for multiple comparisons by the program). GI fragments represent a possible gene conversion event between two sequences in the alignment.

For each GI fragment, we made an educated guess on which sequence was the donor sequence and which was the acceptor sequence by comparing how sequences clustered using sites within the fragment bounds to how sequences clustered using sites outside of the fragment bounds (but in the same exon). Sequences that were confidently determined to be acceptor sequences were excluded from the Bayes factor analyses for the relevant exon because their non-tree-like behavior has the potential to bias results. For sequence pairs where the direction could not be determined, both sequences were excluded from subsequent analyses.

Bayes Factors

Because we could not perform model testing directly on the full phylogenies, we used an alternative approach—computing Bayes factors for TSP within manageable subsets of the data, i.e., quartets of alleles. Let D be a sample of phylogenies from BEAST2, sampled from the posterior with uniform prior. For a chosen species, we have a null hypothesis H, that human alleles form a monophyletic group, and an alternative hypothesis, H^c, that is also the complement of H—that the human alleles do not form a monophyletic group. The Bayes factor, K, is a ratio quantifying support for the alternative hypothesis:

where the first term on the right hand side is the posterior odds in favor of the alternative hypothesis and the second term is the prior odds in favor of the null hypothesis. Bayes factors above 100 are considered decisive support for the alternative hypothesis (Jeffreys, 1998).

Because it is difficult to evaluate monophyly using a large number of alleles, we evaluate Bayes factors considering 4 alleles at a time: 2 alleles of a single species and 2 alleles of different species. For example, to assess support for TSP between humans and chimpanzees, we could use 2 human alleles and 2 bonobo alleles. Or, to assess support for TSP between humans and OWM, we could use 2 human alleles, one baboon, and one macaque allele. Because there are many possible sets of 4 alleles for each comparison, we tested a large number of quartets. We reported the maximum Bayes factor among all tested allele sets to represent evidence for TSP for that species comparison, because our aim was to find any evidence of TSP among any set of 4 alleles.

Next, we calculated the prior odds of the null hypothesis (that the chosen species, i.e. humans, form a monophyletic group). The prior odds , because if the trees were assembled at random, there is 1 possible unrooted tree where the 2 human alleles would form a monophyletic group and 2 possible unrooted trees where the 2 human alleles would not form a monophyletic group, as shown in Figure 7.

Figure 7.

The data, D, is the set of BEAST2 trees, so the posterior odds is the fraction of BEAST2 trees where the 2 human alleles do not form a monophyletic group divided by the fraction BEAST2 of trees where the 2 human alleles do form a monophyletic group. If either fraction is 0, we set its probability to , where n is the number of BEAST2 trees for that gene/sequence subset, and set the complement’s probability to 1 − p. This is the reason that some labels in Figure 3-Figure 4 contain a > sign (e.g. if no trees in a set of 14, 000 were monophyletic, then the Bayes factor must be at minimum 7, 000).

Bayes factors K were then computed as follows and interpreted according to the scale given by Jeffreys (1998).

For each gene and genic region, we tested for TSP between human and chimpanzee, gorilla, orangutan, gibbon, OWM, and NWM. For the Class II genes, which have orthologs beyond the primates, we also tested for TSP between human and tarsier, Strepsirrhini, Flying Lemur, Treeshrew, Glires, Laurasiatheria, and Atlantogenata. For the Class I genes, we considered outgroup species to be the tarsier and Strepsirrhini.

Rapidly-Evolving Sites

BEAST2 places sites into partitions and estimates evolutionary rates for each partition. We averaged these rates over all sampled states, resulting in an overall evolutionary rate for each nucleotide position. In designing the gene groups, we included common “backbone sequences” in every set, although we expanded one particular clade for each focused tree. The inclusion of genes spanning the whole MHC family caused the alignments to contain many gaps. To normalize the rates, we took advantage of the fact that every alignment had many mostly-gap sites. We defined “gappy” sites as those in which the alignment had a gap in all of the following human back-bone alleles: HLA-A*01:01:01:01, HLA-B*07:02:01:01, HLA-C*01:02:01:01, HLA-E*01:01:01:01, HLA-F*01:01:01:01, HLA-G*01:01:01:01, and HLA-J*01:01:01:01 for Class I, HLA-DRA*01:01:01:01, HLA-DQA1*01:01:01:01, HLA-DPA1*01:03:01:01, HLA-DMA*01:01:01:01, and HLA-DOA*01:01:01:01 for Class IIA, and HLA-DRB1*01:01:01:01, HLA-DQB1*02:01:01:01, HLA-DPB1*01:01:01:01, HLA-DMB*01:01:01:0 and HLA-DOB*01:01:01:01 for Class IIB. Because mostly-gap sites are not expected to affect the BEAST2 run very much and were common to all group alignments, we considered these sites’ BEAST2 evolutionary rates as a baseline. For each group and gene region, we calculated the mean rate among all of these baseline mostly-gap sites. Then, we expressed normalized per-site rates as fold-changes by taking the base-2 logarithm of the ratio of each site’s evolutionary rate to the baseline mean.

Protein Structures

To map the rapidly-evolving sites onto protein structures, we first translated our nucleotide alignments into protein sequences using Expasy translate (Gasteiger et al., 2003). We then aligned our translated sequences with amino acid sequences from selected Protein Data Bank (PDB) (Berman et al., 2000)(https://www.rcsb.org/) structures (listed below) using MUSCLE (Edgar, 2004) with default settings.

The translated alignments and the PDB structures’ sequences allowed us to determine which nucleotides in the alignment corresponded to which amino acids in each structure. We then calculated per-amino-acid evolutionary rates by averaging the per-site evolutionary rates among the sites composing each codon. We caution that the vast majority of PDB structures we considered were from human (and were functional), and thus they cannot capture all of the indels and structural changes that define different primate proteins. Additionally, some of our alignments included pseudogenes and diverse sets of related genes, for which it may seem reductionist to map to a single human protein structure. However, we note that despite significant sequence variation, the structures of different MHC proteins (from different genes within the same class) are extremely similar, and thus we do not expect such complications to alter the results.

We used PyMOL version 2.4.2 (Sch, 2021) to visualize the per-amino acid evolutionary rates on each gene’s protein structure. We used model 4BCE (Teze et al., 2014) for HLA-B, 4NT6 (Choo et al., 2014) for HLA-C, 7P4B (Walters et al., 2022) for HLA-E, 5JLZ (Gerstner et al., 2016) for HLA-DR and 2NNA (Henderson et al., 2007) for HLA-DQ from Protein Data Bank (Berman et al., 2000) to prepare the main figures (https://www.rcsb.org/).

We calculated the distances between all atoms of all amino acids of the HLA molecule and all atoms of all amino acids of the peptide in PyMOL, then took the minimum distance to represent each amino acid’s distance to the peptide. Where possible, we averaged the minimum distances over multiple alternative structures for each protein, to avoid relying too heavily on a particular structure. The structures chosen represented different alleles, different contexts (e.g. bound to a receptor or not, bound to a self or non-self peptide), and in a few cases different species.

We used models 1ZVS (Chu et al., 2007), 3JTT (Dai et al., 2010a), 3OX8 (Liu et al., 2011), 3OXR (Liu et al., 2011), 3OXS (Liu et al., 2011), 3RL2 (Zhang et al., 2011), 4HX1 (Niu et al., 2013), 6J1V (Zhu et al., 2019), 6J1W (Zhu et al., 2019), 6MPP (Flores-Solis et al., 2019), 6PBH (van de Sandt et al., 2019), 7SR0 (Finton et al., 2023), 7SRK (Finton et al., 2023), 7WT5 (Asa et al., 2022), and 8I5C (Lu et al., 2023) from Protein Data Bank (Berman et al., 2000) for MHC-A; 1JGD (Hillig et al., 2004), 3BVN (Kumar et al., 2009), 3KPL (Macdonald et al., 2009), 3KPN (Macdonald et al., 2009), 3LN4 (Bade-Döding et al., 2011), 3LN5 (Bade-Döding et al., 2011), 3RWJ (Wu et al., 2011), 3W39 (Yagita et al., 2013), 3X13 (Saunders et al., 2015), 4JQV (Rist et al., 2013), 4JRY (Liu et al., 2013), 4MJI (Motozono et al., 2014), 4O2E (Sun et al., 2014), 4PRA (Liu et al., 2014), 4PRB (Liu et al., 2014), 5EO0 (Du et al., 2016), 5IEK (Alpizar et al., 2016), 5VUD (Illing et al., 2018), 5VVP (Illing et al., 2018), 5VWF (Illing et al., 2018), 6IWG (Yamamoto et al., 2019), 6MTM (Grant et al., 2018), 6PYL (Sian et al., 2019), 6PYV (Sian et al., 2019), 6UZP (Schutte et al., 2020), 6VIU (Schutte et al., 2020), 6Y27 (Loll et al., 2020), 7R7V (Li et al., 2023), 7T0L (Vivian and Rossjohn, 2022), 7TUC (Jiang et al., 2022b), 7X1C (Huan et al., 2023), and 7YG3 (Jiang et al., 2022a) for MHC-B; 4NT6 (Choo et al., 2014), 5VGD (Kaur et al., 2017), 5VGE (Kaur et al., 2017), 5W67 (Mobbs et al., 2017), 6PAG (Moradi et al., 2021), and 7WJ3 (Asa et al., 2022) for HLA-C; 2ESV (Hoare et al., 2006), 3CDG (Petrie et al., 2008), 5W1V (Sullivan et al., 2017), 7P49 (Walters et al., 2022), and 7P4B (Walters et al., 2022) for HLA-E; 5IUE (Dulberger et al., 2017) for HLA-F; 1YDP (Clements et al., 2005), 2DYP (Shiroishi et al., 2006), and 3KYN (Walpole et al., 2010) for HLA-G; 1HDM (Mosyak et al., 1998), 2BC4 (Nicholson et al., 2006), 4FQX (Pos et al., 2012), 4GBX (Pos et al., 2012), and 4I0P (Guce et al., 2013) for HLA-DM; 4I0P (Guce et al., 2013) for HLA-DO; 3LQZ (Dai et al., 2010b), 3WEX (Kusano et al., 2014), 7T2A (Ciacchi et al., 2023), 7T6I (Sebastian Klobuch et al., 2022), and 7ZAK (Racle et al., 2023) for HLA-DP; 2NNA (Henderson et al., 2007), 4D8P (Tollefsen et al., 2012), 5KSA (Petersen et al., 2016), 5KSU (Nguyen et al., 2017), 6DIG (Jiang et al., 2019), and 6PX6 (Ting et al., 2020) for HLA-DQ; and 1BX2 (Smith et al., 1998), 1FV1 (Li et al., 2000), 1H15 (Lang et al., 2002), 1T5X (Zavala-Ruiz et al., 2004), 2Q6W (Parry et al., 2007), 3C5J (Dai et al., 2008), 4FQX (Pos et al., 2012), 4H1L (Yin et al., 2012), 5JLZ (Gerstner et al., 2016), 5V4M (Ooi et al., 2017), 6ATF (Scally et al., 2017), and 8EUQ (Kassardjian et al., 2023) for HLA-DR.

Disease Literature

We conducted a literature search for papers that used HLA fine-mapping to discover disease associations, limiting our selection to those including at least 1,000 cases and which identified putatively independent signals via conditional analysis. We included all independent amino acid signals identified as significant by the original authors. If there was more than one study for the same disease, but in different populations, we included all unique independent hits. We also collected associations between amino acids and TCR phenotypes.

Polymorphic Amino Acid Positions

Amino acid allele frequencies were obtained from Luo et al. (2021), who imputed HLA variation in 1000G using SNP2HLA (https://github.com/immunogenomics/HLA-TAPAS/blob/master/resources/1000G.bglv4.FRQ.frq).

Acknowledgements

We acknowledge support from NIH grants R01 HG011432 and R01 HG008140. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1656518. We appreciate helpful comments from Jeffrey Spence, the Pritchard lab, and the reviewers of the previous version of this work.