Immunogenomics: Molecular hide and seek

Abstract

Similar to other classical science disciplines, immunology has been embracing novel technologies and approaches giving rise to specialised sub-disciplines such as immunogenetics and, more recently, immunogenomics, which, in many ways, is the genome-wide application of immunogenetic approaches. Here, recent progress in the understanding of the immune sub-genome will be reviewed, and the ways in which immunogenomic datasets consisting of genetic and epigenetic variation, linkage disequilibrium and recombination can be harnessed for disease association and evolutionary studies will be discussed. The discussion will focus on data available for the major histocompatibility complex and the leukocyte receptor complex, the two most polymorphic regions of the human immune sub-genome.

Introduction

The ongoing 'hide and seek' between pathogens and their host immune systems has led to a molecular arms race which makes immunology one of the most interesting, but also one of the most complex, areas of research. By contrast with the harmless children's game, the outcome of this molecular race determines health versus disease and, in many cases, survival versus death. A prerequisite for immunogenomic research is the availability of comprehensive and fully informative sequence variation and gene maps. A high-quality sequence map of the human genome and successive generations of annotation have been available since 2001 [1, 2]. Using this sequence as a reference, global and population-specific variation maps of increasing resolution were subsequently generated by the International SNP Consortium and HapMap Projects [3, 4]. Because of the extreme polymorphism encountered in some immune regions, separate efforts were focused on regions such as the major histocompatibility complex (MHC; reviewed by Allcock et al. [5]), the leukocyte receptor complex (LRC; reviewed by Yawata et al. [6]) and others, to complement the global variation map. In addition, several specialised databases, including the Immuno Polymorphism Database [7] and the International Immunogenetics Information System [8], catalogue many of the classical immune genes and their allelic variants. Using immune ontology definitions [9], the first comprehensive gene map and database of the human immune sub-genome was recently reported to consist of 1,562 genes (about 7 per cent of human genes) that are distributed at varying densities across all chromosomes except the Y chromosome [10].

In the following review, some of the benefits and caveats of using a holistic approach to immunogenomic analysis will be discussed with respect to genetic and epigenetic variation and linkage disequilibrium (LD).

Genetic variation

Based on both experimental data [11–15] and inferences from population genetics data [16–18], there is evidence that recombination is not uniformly distributed across the human genome. In addition to the heterogeneity of recombination rates, it was suggested that the human genome could be subdivided into relatively short fragments with little or no evidence of historical recombination [19]. Under this premise, a genome-wide LD map is being constructed by the International HapMap Project, which aims to generate a resource of DNA variation in human populations with different ancestry [3]. In its first phase, which is close to completion, genotyping data at one single nucleotide polymorphism (SNP)/five kilobase (kb) density will provide enough resolution to portray the LD landscape of the genome, including regions of the immune sub-genome.

LD maps of complete chromosomes are already available; [20–22] however, detailed LD patterns at high SNP density covering large chromosomal regions encoding immunerelated gene clusters have long been awaited. LD data have so far been reported for two such clusters: the LRC and the MHC. Within the LRC, the haplotypic nature of the killer immunoglobulin-like receptor (KIR) gene cluster causes a DNA size variation of 120-200 kb as a consequence of the presence/absence of genes in some haplotypes [23–26]. This feature, in addition to intron similarity owing to recent duplications, restricts the availability of informative SNPs for LD studies [27] and poses particular restrictions on large-scale genotyping protocols. Consequently, a detailed LD map for the LRC is not available. For the MHC, a comprehensive LD map is available [16] as a direct consequence of existing SNP resources [28].

Historically, considerable attention has been paid to the MHC region, resulting in a wealth of valuable information in a variety of areas -- that is, origin and maintenance of nucleotide diversity in classical human leukocyte antigen (HLA) genes; gene conversion, average recombination rates, presence and distribution of the recombination hot spots; and extended LD in long-range haplotypes [11, 14, 15, 29–31]. The high-resolution LD map available for the MHC permits the detailed study of the allelic correlation, which can be particularly useful for the selection of tagSNP sets for disease association studies. Since the history of recombination between two SNPs can be estimated using D', DNA segments showing no evidence of historical recombination have been defined as LD-blocks [19]. Consequently, haplotypes within such LD-blocks are likely to share a common ancestral haplotype, and the genealogy of an LD-block might be different from the genealogy of neighbouring LD-blocks. Phylogenetic analysis of haplotypes observed in a ~110 kb LD-block (rs3115569-rs2022533) at the MHC class II-III boundary, the region of lowest recombination rate in the MHC (0.0093 centimorgan [cM]/ megabase [Mb]), surrounded by two recombination hot spots [16], showed that the pattern of nucleotide differences between haplotypes within this LD-block is difficult to account for by mutation only, but could readily be explained by recombination events (Miretti et al., unpublished). Following detailed inspection of the haplotype nucleotide sequence alignment at variable positions (Figure 1B), putative recombination sites can be identified. The subsequent splitting of the LD-block at a recombination site (Figure 1C) resulted in haplotypic networks presenting distinct clustering in both fragments (Figure 1D), which were also dissimilar from that of the complete LD-block. While it needs to be examined in a larger sample of LD-blocks from diverse populations, the presence of recombination within LD-blocks suggests that the occurrence of LD between two SNPs might not be sufficient to unambiguously detect historical recombination beyond the history of the sample population where LD is being measured. Recombination hot spots within LD-blocks -- that is, hot spots that have not left an imprint on LD -- have recently been postulated to be either old hot spots en route to extinction or too young to leave a mark on haplotype diversity in Europeans [18]. Additional experimental data based on sperm cell recombination analysis would help to confirm the signature of recombination identified here and to determine if it is an example of a potentially evolutionarily very young hot spot [18].

Figure 1

Potential recombination sites within a linkage disequilibrium (LD) block. Patterns of LD for the complete major histocompatibility complex were obtained for a panel of 180 Centre d'Etude du Polymorphisme Humain founder chromosomes [16] and LD-blocks defined according to Gabriel et al. criteria [19] employing Haploview [32]. (A) The region selected for the analysis constitutes an LD-block extending over a ~110 kilobase (kb) DNA segment (National Center for Biotechnology Information built 35 coordinates: 32321766-32431723) containing most of the C6ORF10 gene, with an average space between single nucleotide polymorphisms (SNPs) of 1.6 kb. This region also represents a recombination cold spot -- with recombination rates = 0.00931 centimorgans/megabase, ten times lower relative to that on the neighbouring area -- surrounded by two recombination hot spots identified in NOTCH4 and C6ORF10 genes [16]. (B) Ten haplotypes (rows) observed in this sample, which were recognised based on alleles at 68 variable positions (SNPs, in columns) within the selected LD-block. Haplotype names and frequencies (per cent) are given on the left. (C) A break in the haplotype alignment suggesting the presence of ancestral recombination within the LD-block. Haplotypes were then split at a potential recombination site into two sub-regions (with 41 and 27 SNPs, respectively) aiming to construct a phylogenetic network in order to check for consistency. (D) The distribution of variation among haplotypes as represented by median-joining trees obtained for both sub-regions employing the NETWORK program [33]. Haplotypes F and D are grouped closely related within each network but belong to different clusters when comparing networks from the two sub-regions. The number of substitutions involved in these branching differences could preferably be explained by recombination events rather than owing to mutations occurring between haplotypes. The distribution pattern of these differences adds consistency to this view and strongly suggests the presence of recombination. The circle size is proportional to the haplotype frequency, and the number of substitutions between nodes is represented by the number of lines between them. The site where the haplotypes were split can be changed to few positions on either side of the current site without significantly modifying the topology of the network.

Overall recombination rates within the MHC are known to be lower than across chromosome 6 and the genome average [34, 35], and regions of relatively high recombination rates (hot spots) and low recombination rates (cold spots) have been described in males based on microsatellite typing [15]. The fine-scale experimental recombination analysis in a 200 kb segment within the MHC class II region revealed that a few hot spots shape the distribution of LD throughout the segment, that hot spots tend to be delimited to short DNA segments (usually < 5 kb) and that a portion of the recombination is contributed by gene conversion [11, 31, 36, 37]. The question therefore arises as to whether the recombination pattern observed in this segment can be extended to the complete MHC and, ultimately, if it represents a general characteristic of the genome. The heterogeneous model of recombination seems to be a general feature of the genome [17], where the landscape is moulded by the local distribution and intensity of recombination. According to the recombination rates inferred from population genetic data [16], the experimentally verified recombination hot spots in the 200 kb segment of the MHC class II region [11] might not represent a paradigm for the entire MHC, as their intensities are ~ 2-100 times higher than hot spots in the rest of the MHC. Figure 2 shows the uneven distribution of hot spots across the complete MHC -- integrated with gene annotation, LD-blocks and tagSNP distribution, as represented by the GLOVAR genome browser http://www.glovar.org -- and also that the hot spot intensity is far from being uniform. Additional large-scale experimental evidence of hot spot intensity and distribution is required to more accurately assess the concordance with predictions based on population genetic data. Importantly, hot spot features will significantly influence the selection of tagSNPs for disease association studies involving any particular region of the immune sub-genome.

Figure 2

Data integration in the major histocompatibility complex (MHC) region as represented by the GLOVAR genome browser showing the distribution and intensity of recombination hot spots relative to gene annotation, linkage disequilibrium blocks and tag single nucleotide polymorphism distributionwww.glovar.org. Note that recombination hot spots within the MHC class II region are comparatively more intense than in the rest of the MHC. Most gene names were excluded owing to space constraints.

From the detailed distribution of recombination and LD across the MHC, it would be possible to exploit the common variation observed in approximately 80 per cent of the MHC sequence, which is in high LD, to conduct association studies. Even with the current SNP density (1 SNP/1.9 kb), however, between 10-27 per cent of the DNA sequence in the MHC sub-regions is contained in regions of low LD, where tagSNP efficiency -- the number of genotyped markers divided by the number of tagging SNPs [38] -- would not entail any benefit, meaning that almost every SNP would need to be typed. Figure 2 illustrates the dependence of the tagSNP distribution on the SNP density and LD pattern, which is particularly variable across the MHC. Long regions of strong LD are common in the extended class I sub-region, where 75 tagSNPs are sufficient to recover haplotypic information provided by 408 SNPs. Conversely, 204 tagSNPs are necessary to recover the equivalent information from 466 SNPs in the class II sub-region, where the LD pattern is more interrupted and the decay of LD is more pronounced [16]. Also, the tagSNP sets are derived from studies which mostly dismiss loci of minor allele frequency (MAF; < 5 per cent), which seem to be present at increased frequency in low LD regions [39], thus excluding rare variants from the analyses. This raises the question of whether failure to detect disease variants simply results from the exclusion of SNPs with rare alleles from the analysis (discussed below). Furthermore, tagging effectiveness -- the proportion of 'hidden' SNPs being detected in LD with SNPs from the tagging set -- can be substantially increased by incorporating SNPs with rare alleles (MAF < 5 per cent); this effect was even more prominent in low LD regions [39].

Epigenetic variation

Much less is known about epigenetic variation, both in terms of causality and function. The most frequent and stable form of epigenetic variation is differential DNA methylation. DNA methylation, which was first discovered in 1948, occurs naturally at the carbon-5 position of cytosine (5-methyl cytosine) at CpG dinucleotides [40]. In subsequent years, it was proposed that DNA methylation played an important role in the regulation of gene expression [41, 42] and disease aetiology, particularly cancer [43]. The discovery of CpG islands (sequences enriched in cytosine-guanosine dinucleotides) suggested candidate regions in the genome for epigenetic modulation [44]. DNA methylation occurs predominantly at CpG sites [45] and, in addition to gene regulation, is also involved in phenomena such as X-chromosome inactivation in female mammals, parent-of-origin-specific, mono-allelic gene expression (imprinting) and epigenetic reprogramming in mammalian development [46]. The key technology for detecting 5-methyl cytosine by DNA sequencing was developed in 1992, and is known as bisulphite sequencing [47].

The need to consider epigenetic variation alongside genetic variation in the context of disease has been highlighted by the finding of high discordance rates in monozygotic twin studies and numerous other studies, confirming that epigenetic factors play a decisive role in the aetiology of virtually all human pathologies (reviewed by Robertson and Wolffe) [48]. Therefore, when the SNP Consortium announced in 1999 that it would generate a first-generation genetic variation map of the human genome [49], the time and opportunity seemed right to generate a first-generation epigenetic (methylation) variation map alongside -- which is what the Human Epigenome Consortium announced that it would do [50]. One aspect of the Human Epigenome Project (HEP) that is of particular interest in the context of this review, is the aim and ability to identify methylation variable positions (MVPs), which, together with SNPs, promise to significantly advance our ability to understand and diagnose human disease. MVPs are defined as differentially methylated CpG sites that have the statistical power to discriminate, for example, between biological states such as active versus inactive or healthy versus diseased.

As a result of the HEP and other studies, a detailed genomic methylation map is already available for the MHC [51] and models for the epigenetic control of the KIR expression repertoire have been proposed (reviewed by Uhrberg [52]). For the MHC, methylation profiles have been generated essentially for all expressed genes, demonstrating tissue-specificity (eg for the C2 locus) and inter-individual heterogeneity (eg for the TNF locus). Previously, it has been shown that transcription profiles can be associated with specific haplotypes [53], with epigenetic states [54] or even with specific epi-alleles in the absence of DNA variation [55]. Overall, the methylation profile of the MHC appears to be strongly bimodal, with good correlation between hypo- and hypermethylation within the upstream regions of genes and their transcriptional activity (Figure 3). For the LRC, fewer data are available and are mostly restricted to the KIR genes. KIR genes exhibit a high degree of polymorphism and are expressed in a clonally restricted fashion [56]. Each cell expresses a mono-allelic repertoire which is highly variable with respect to number and combination of receptors. KIR genes are regulated by three types of promoters, corresponding to the three different modes of expression: one for KIR2DL4, which is constitutively expressed; one for KIR3DL3, which is expressed at a very low level; and one which is common to all clonally distributed KIR genes [57]. Given the absence of promoter variability among the clonally distributed KIR genes, it is, in fact, likely that the variegated expression of these genes is mainly regulated epigenetically, and a model has been proposed to this effect [52]. The model proposes four stages of sequential DNA and histone modifications leading to the observed mosaicism of clonally restricted KIR expression but requires some still unproven processes, including locus- and allele-specific DNA demethylation.

Figure 3

Correlation between DNA methylation and gene expression. DNA methylation data derived from the Human Epigenome Pilot Project were grouped according to 5'-UTR being methylated or unmethylated and the body of the gene being methylated or unmethylated. Each group was then plotted against the corresponding Genomics Institute of the Novartis Research Foundation (GNF) http://symatlas.gnf.org/SymAtlas expression level, which is set to anything below 200 as not being expressed and anything above 200 as being expressed. This cluster analysis shows that genes with methylated 5'-UTRs are silent and genes with unmethylated 5'-UTRs are expressed. In this study, no correlation was found between expression and the methylation level within the body of the gene. Reproduced from Rakyan et al. [51]

Compared with genetic data, the ability to analyse and interpret epigenetic data is still fairly simplistic, although enormous progress has been made, particularly in the past few years. What is becoming increasingly clear, however, is that the two are inextricably intertwined, and efforts towards an integrated (epi)genetic approach to common disease are well underway [58, 59].

Application of immunogenomic data for disease association studies

Major progress in genotyping technology has improved our ability to harvest genotypic information from a comprehensive number of loci and has laid down the foundation for the genome-wide common variation survey [3], enhancing our understanding of human genetic variation. Current genotyping costs, however, make genome-wide association studies still challenging for the large sample sizes that will be required for adequately powered studies. The interdependencies of marker numbers, allele frequency, LD, cohort size and power in disease association studies have recently been reviewed [60–62] and will therefore not be discussed in depth here. Essentially, genome-wide studies are still restricted by the allelic spectra underlying complex diseases. The frequency distribution of the disease variants and the proportion of the trait variance for which they are responsible determine the potential power of genetic association studies and, therefore, the feasibility of the study. Disease susceptibility variants associated with Mendelian disorders tend to be modelled by purifying selection and present low population frequencies ( < 1 per cent) [63]. Thus, rare alleles -- also including lowfrequency, mildly deleterious variants together with those alleles in high LD -- might be under-represented in genomic LD surveys and tagSNP sets excluding loci with MAF < 5 per cent. It has been proposed, however, that genetic risk for common diseases and most of the clinically important traits are determined by the joint contribution of diseasepredisposing high-frequency alleles that are shared between unrelated affected individuals in the population: the common-disease/common-variant (CD/CV) model [60, 64, 65]. A less extreme model based on the divergence of the allelic spectra of disease susceptibility variants relative to that of all variants has been proposed for evaluating the impact of allele architecture on common diseases [61]. Ultimately, the applicability of the CD/CV model in disease association studies will eventually rely on the selected markers and on the relationship between the four parameters that affect the apparent size at the marker locus, namely the odds ratio of the disease allele, the disease allele frequency, the marker allele frequency and the LD between the marker and the disease locus (reviewed by Zondervan and Cardon [62]).

Diseases related to defective immune response, such as autoimmune disorders, might comply with the CD/CV model of allelic frequency distribution for susceptibility loci, showing increased allele frequency after being under positive selection for infectious disease resistance or heterozygote advantage [61, 66]. Direct extrapolations to the immune sub-genome and autoimmune diseases could be elusive, however; that is, it would be difficult to prove that alleles having an effect on the function of an immune gene -- and also contribute to the risk of autoimmune diseases -- are beneficial on a population level because they increase the diversity in the immune response [67]. In fact, the overall level of polymorphism in immune genes is, in general, similar to that of non-immune genes [10], and only a few loci among the genes constituting the immune sub-genome largely contribute to the diversity of the immune response -- disregarding loci with somatic recombination. In addition, modes of selection maintaining extreme allelic diversity levels -- heterozygote advantage, balanced selection -- have been proposed for these few loci, including the classical HLA class I and II and KIR genes [29, 30, 66, 68]. Consequently, it is probably premature to generalise any disease model for immune genes. A recent study shows, for instance, that multiple disease loci with predominant disease variants are not necessarily required to develop a complex immune disorder [69]. It has been shown that sarcoidosis -- a multi-systemic immune disorder initially associated with several classical HLA markers -- is associated with a single disease locus independently responsible for most of the predisposing influence [69]. The authors of this study demonstrated that a single transition in exon 5 of the BTNL2 gene, which alters a splice site and causes a premature stop codon, has a profound effect on the mature protein structure and function.

A number of autoimmune disease-susceptibility loci have been mapped to the MHC and LRC regions [70–72]. Autoimmunity and cancer are associated with the progressive decline in immune functions, such as decline in T cell function, dysregulation of T cell apoptosis and immune senescence during ageing [73, 74]. Tolerance induction of T cells is mediated by cytotoxic T lymphocyte antigen 4. Polymorphism in this inhibitory receptor has also been associated with autoimmune diseases [75, 76].

Outlook

Since the term haplotype was coined in 1967 [77], a vast amount of literature has been published employing haplotypes to construct genetic and LD maps, which have enabled the emerging field of immunogenetics to make major contributions to immunology, population genetics and medicine [78]. For the MHC, a haplotype hierarchy has been defined based on pedigree and population genetic data, where physically linked classical HLA genes constitute 'blocks' of relatively short DNA segments containing alleles in LD at different loci, namely HLA-B/C, TNF, complotype and HLA-DQB/DRB. Combinations of these four basic blocks are frequently found in LD, constituting 'conserved extended haplotypes' (CEH) [79, 80] or 'ancestral haplotypes' [81]. More recently, local fine-scale LD blocks based on high-density SNP typing have been described and meiotic recombination hot spots inferred from population genetic data. The integration of long-range LD in CEH and local fine-scale LD data will provide new and exciting opportunities to explore the evolution and disease associations of the immune sub-genome. While the jury is still out on the validity of the CD/CV model, the recently announced Wellcome Trust Case Control Consortium (WTCCC) is likely to provide the answer to this ongoing debate. The WTCCC proposes to analyse 19,000 samples for common genetic variants in eight common diseases using data generated by the HapMap Consortium. Integration of the genetic and epigenetic datasets discussed here is a necessary next step towards a more holistic approach to immunological research, which is one of the aims of the recently formed International Immunomics Society, IMMIS http://research.i2r.a-star.edu.sg/IIMMS/. The IMMIS was formed with the main objective of promoting the science of immunomics, which is the interdisciplinary field spanning immunology, immunoinformatics, genomics, proteomics, bioinformatics and related scientific fields.