Common variants in the HLA-DRB1-HLA-DQA1 HLA class II region are associated with susceptibility to visceral leishmaniasis.

To identify susceptibility loci for visceral leishmaniasis, we undertook genome-wide association studies in two populations: 989 cases and 1,089 controls from India and 357 cases in 308 Brazilian families (1,970 individuals). The HLA-DRB1-HLA-DQA1 locus was the only region to show strong evidence of association in both populations. Replication at this region was undertaken in a second Indian population comprising 941 cases and 990 controls, and combined analysis across the three cohorts for rs9271858 at this locus showed P(combined) = 2.76 × 10(-17) and odds ratio (OR) = 1.41, 95% confidence interval (CI) = 1.30-1.52. A conditional analysis provided evidence for multiple associations within the HLA-DRB1-HLA-DQA1 region, and a model in which risk differed between three groups of haplotypes better explained the signal and was significant in the Indian discovery and replication cohorts. In conclusion, the HLA-DRB1-HLA-DQA1 HLA class II region contributes to visceral leishmaniasis susceptibility in India and Brazil, suggesting shared genetic risk factors for visceral leishmaniasis that cross the epidemiological divides of geography and parasite species.

Leishmania are protozoan parasites that live in macrophages. They are transmitted by sand flies and cause severe and debilitating cutaneous, as well as fatal visceral disease in subtropical/tropical regions of Old and New Worlds. Leishmaniasis is classified by WHO as one of the neglected tropical diseases. It affects 12 million people and there are an estimated 1.5 million new cases annually[1]. Of these, 500,000 are cases of potentially fatal visceral leishmaniasis caused by the Leishmania donovani complex, 90% of which occurs in three foci in India/Bangladesh/Nepal, Sudan, and Brazil. Skin-tests and lymphocyte responses indicate that only 1 in 5-10 infected individuals develop clinical disease[2-4]. The importance of host genetic factors is indicated by familial clustering[5] and high sibling risk ratios[6]. However, human genetic studies undertaken to date (reviewed[7-9]) provide inconsistent results. As part of the Wellcome Trust Case Control Consortium 2 (WTCCC2) study of 15 complex disorders and traits, we report the first genome-wide association study (GWAS) of visceral leishmaniasis across major foci of disease caused by L. donovani in India and L. infantum chagasi in Brazil.

Subjects for the India discovery GWAS (Supplementary Table 1, online methods) were recruited from Bihar state in northeast India. Patients and controls were matched for self-reported age, sex, religion, caste and geographic region of recruitment. The Brazilian family-based sample was collected as part of the Belém Family Study and from study sites near Natal[6,10,11] (Supplementary Table 1, online methods).

All individuals were genotyped at the Wellcome Trust Sanger Institute on the custom Illumina Human660W-Quad chip (online methods). After performing stringent quality control (QC) procedures[12] (online methods) the Indian discovery dataset comprised 2078 individuals (989 cases and 1089 controls) genotyped at 526,731 SNPs and the Brazilian discovery dataset comprised of 1970 individuals (357 cases; 308 families) genotyped at 553,323 SNPs.

As ancestry differences and close relationships between case and control individuals can confound association studies, we used a variance components method (sometimes referred to as a mixed model), which models the pair-wise relatedness between individuals to account for population structure in the data[13]. This approach accounts for relatedness on different scales, from close relatives to distant ancestral structure, and it can thus be applied both to the case-control samples in India and the family data from Brazil. Similar models have been recently used in several other association studies[14-18]. An advantage of our implementation of the mixed model is that we are able to estimate the effect sizes on the log-odds scale. Using this method, the genomic inflation factor (λ) was 1.03 for the India analysis and 1.02 for the Brazil analysis, showing that it successfully handled any population structure or relatedness in the data. The mixed model approach was used for all analyses presented here, for both discovery and replication data and in further dissection of association signals.

Previous studies attempting to identify the genetic components of visceral leishmaniasis susceptibility have typically been small, underpowered, family studies or candidate gene studies[19-23]. We present the results from our discovery GWAS and replication data (see below) at these previously identified loci in Table 1. None of the loci previously found to be associated with visceral leishmaniasis consistently show an association with visceral leishmaniasis in our three datasets.

Table 1

Signals of association at previously reported non-HLA loci for visceral leishmaniasis

SNP	Chr	Locus	Position (bp)	RiskAllele	Popb	Ref	This Study	OddsRatio	95% CI	P-value
rs4674259	2q35	CXCR2	218,699,250	A	India	19	Indian Discovery	0.92	0.81-1.04	0.195
							IndianReplicationc	0.87	0.77-1.00	0.040
							Brazil Discovery	1.24	1.04-1.49	0.017
rs4672875a	2q35	CXCR1	218,729,630	A	India	19	Indian Discovery	1.05	0.89-1.23	0.582
							IndianReplication	NA	NA	NA
							Brazil Discovery	0.82	0.65-1.05	0.104
rs2276631	2q35	SLC11A1	218,957,257	A	Sudan	23	Indian Discovery	NA	NA	NA
							IndianReplication	0.96	0.80-1.15	0.633
							Brazil Discovery	NA	NA	NA
rs3731865	2q35	SLC11A1	218,958,247	G	Sudan	23	Indian Discovery	NA	NA	NA
							IndianReplication	0.99	0.82-1.19	0.933
							Brazil Discovery	NA	NA	NA
rs17860704	6q27	DLL1-FAM120B	170,444,354	G	Sudan	22	Indian Discovery	NA	NA	NA
							IndianReplication	1.14	0.97-1.34	0.117
							Brazil Discovery	NA	NA	NA
rs159273	17q12	CCL1	29,714,798	C	Brazil	21	Indian Discovery	1.05	0.95-1.15	0.362
							IndianReplication	NA	NA	NA
							Brazil Discovery	1.15	0.86-1.56	0.334
rs2063979	17q12	CCL16	31,327,679	A	Brazil	21	Indian Discovery	0.99	0.87-1.13	0.913
							IndianReplication	1.05	0.92-1.19	0.467
							Brazil Discovery	0.92	0.77-1.10	0.335
rs854680	17q12	CCL16	31,333,164	C	Brazil	21	Indian Discovery	1.00	0.88-1.13	0.938
							IndianReplication	1.07	0.94-1.22	0.277
							Brazil Discovery	1.16	0.97-1.39	0.103
rs228953	22q12	IL2RB	35,861,382	G	Sudan	20	Indian Discovery	NA	NA	NA
							IndianReplication	0.97	0.85-1.10	0.615
							Brazil Discovery	NA	NA	NA

Odd ratios are given for the previously reported risk allele. Positions are in NCBI human genome build 36 coordinates. ‘NA’ indicates that the SNP was not genotyped in this sample.

r2=1 (HapMap Utah residents with ancestry from northern and western Europe (CEU)) and 0.71 (HapMap Gujarati Indians in Houston, Texas (GIH)) with previously reported rs3138060.

Population in which the initial SNP association was reported.

The Indian replication sample reported in the table overlaps significantly with the Mehrotra et al (2011)[19] case-control population in which rs4674259 was genotyped as part of a separate Sequenom panel.

Because different pathogen species are responsible for disease in India and Brazil, we first examined each GWAS separately. Our primary focus will be on regions showing association in both GWAS studies. This would be appropriate if similar host genetic factors affected susceptibility to the two different parasite species. Observing the same association in different populations, with different study design (population-based and family-based association, respectively) also adds to the robustness of potential findings.

We combined the India and Brazil GWAS through a fixed effects meta-analysis (Fig. 1c). Two regions had a Pcombined <10−6, one of which was the HLA-DRB1-HLA-DQA1 class II region (Supplementary Table 4). Details of the regions showing association (P<10−5) in only one of the two studies are shown in (Fig 1, and Supplementary Tables 2, 3 and 4). Apart from the MHC, none of the association regions with P <10−5 in the separate India or Brazil analysis showed significant association in the combined analysis.

Figure 1

Plot of genome-wide association results in the separate (a,b) and combined (c) discovery GWAS using a variance components method. SNPs in red show regions with replicated association to visceral leishmaniasis susceptibility. a) Indian discovery data at 526,731 SNPs b) Brazilian discovery data at 553,323 SNPs c) Plot of the meta-analysis genome-wide association results.

A second Indian cohort was used in a replication analysis of the findings from the two discovery GWAS. These samples were genotyped on the Immunochip [24], a purpose-designed Illumina chip, some content of which was chosen by WTCCC2 for deep replication of its GWAS findings. After QC (online methods), the Indian replication dataset comprised 1931 individuals (941 visceral leishmaniasis cases and 990 controls). Outside of the MHC, none of the regions that were associated at P<10−5 in either the separate Indian or Brazilian GWAS or the combined discovery meta-analysis replicated at P<0.1 in the Indian replication data (Supplementary Tables 2, 3 and 4). However, only four of the regions associated in the Brazil GWAS were available on the Immunochip.

The MHC was the only region which replicated in all three cohorts (rs9271858, Pcombined=2.76×10−17 and OR(95%CI)=1.41(1.30-1.52)). Interestingly while there is shared signal at the MHC in both discovery GWAS, which clearly replicates, there are also differences between the India and Brazil results (Figure 2, Supplementary Table 4).

Figure 2

Regional association plots of the signal at the MHC in the a) India discovery b) Brazil discovery c) Meta analysis of the discovery cohorts. The −log10 P values are shown on the upper part of each plot. SNPs are colored based on their r2 with the labeled hit SNP, calculated in the a) Indian controls b) Brazilian founders c) no r2 colouring. The SNP coloured in green is the hit SNP in the a) Brazil cohort b) India cohort c) India and Brazil cohorts. The bottom section of each plot shows the fine scale recombination rates estimated from individuals in the CEU HapMap population, and genes are marked as horizontal blue lines. Genes flanking the hit SNP are coloured red and are labeled.

In the India discovery GWAS, the two SNPs most associated with visceral leishmaniasis in the MHC were in strong LD (rs9271842 and rs9271858, r2=1 in Indian controls). Only one of these SNPs was present in the replication data (rs9271858), therefore for convenience we focus below on this SNP, and refer to it as the top SNP. This SNP is also significant in the Brazil data (PBrazil=2.04×10−4), but conversely the top MHC signal in the Brazil data (rs9268878) was not significant in the India data (PIndia=0.71, PReplication=0.58 at rs4428528, r2=1). The top MHC SNP signal in Brazil is 163.6 Kb upstream of the India signal, and is in almost complete linkage equilibrium with the top India SNP (r2= 0.06 in the Brazil founders). Further replication in a Brazilian dataset is required to assess this additional MHC signal: it could be a population or parasite species-specific association, or alternatively a false-positive finding.

In the Indian discovery cohort, when conditioning on the top SNP from the Indian discovery GWAS (rs9271858), there was still residual association signal at the MHC; two highly correlated SNPs both had P<10−3 (rs9271252 and rs9271255, r2=1 in Indian controls). For convenience, because it is typed in the discovery and replication samples, we focus in what follows on one of these, rs9271255, with a P value of 8.79×10−4 in the conditional analysis. This SNP is not highly correlated with rs9271858 (r2=0.21 in Indian controls). When conditioning on both rs9271858 and rs9271255, no association signal remains with P <10−3. In the India replication data, when conditioning on rs9271858, rs9271255 is still significant (P=5.72×10−3). This conditional analysis was repeated in the Brazil discovery data; when conditioning on rs9271858, rs9271255 was significant with a P value of 1.91×10−5.

We turn now to a more detailed examination of this replicated MHC signal (at rs9271858 and rs9271255). We focus on the Indian data, for which the Indian replication cohort is from a comparable population and can therefore be used to replicate cohort specific association signals. To further explore the Indian MHC signal, we estimated phase for the SNP genotypes for all individuals in the Indian discovery data. We subsequently fitted a model to the phased haplotypes with different risk parameters for the four haplotypes defined by the phased genotypes at the two aforementioned SNPs, rs9271858 and rs9271255. After additional exploration of the model (Supplementary Results), it was apparent that combining two of the four haplotypes, to give a model with three risk classes, fitted the data as well as the model with four risk classes (Supplementary Figure 1), and also gave a better fit than the two-SNP model (Supplementary Results).

At rs9271858 and rs9271255 respectively, the haplotype with highest risk is defined by the alleles GA and the haplotype with intermediate risk by AA, with the haplotypes AG and GG being protective relative to the other two haplotypes, and having statistically indistinguishable risk. This model with three risk classes had P=1.65×10−11 in the Indian discovery data, and where the protective haplotypes are set as the baseline with OR=1, the intermediate haplotype OR(95%CI)=1.28(1.09-1.50) and the highest risk haplotype OR(95%CI)=1.72(1.48-2.00). This analysis replicated in the Immunochip data: the model with three risk classes was significant at P= 4.78×10−6, where the intermediate haplotype OR(95%CI)=1.19(1.00-1.40) and the highest risk haplotype OR(95%CI)=1.50(1.27-1.76).

Our results establish that common polymorphisms in the HLA-DRB1-HLA-DQA1 region of the MHC are genetic risk factors for visceral leishmaniasis, and importantly appear to cross the epidemiological divides of geography and parasite species. For GWAS associations, the top SNPs are not necessarily the causal variants, but are expected to be correlated with the functional variants. There are SNPs in the MHC which are in moderate LD with rs9271858 (r2 ~ 0.2) and are located up to 500kb away (e.g. rs408359, r2=0.19 calculated in the Indian controls). This extensive linkage disequilibrium is common in the MHC, and makes it particularly difficult to assess possible functional candidates within the region of association[25].

One possibility is that the association signal could be driven by functional variation in the adjacent highly polymorphic classical HLA class II genes, DRB1 and DQB1, each of which is a natural immunological candidate. The classical alleles at these loci were genotyped in a subset of individuals in both of the discovery cohorts, and phased (statistically) onto the haplotypes estimated from all the SNP data (see Supplementary Results). In this phased data set, chromosomes carrying the protective haplotype at rs9271858 and rs9271255 always carried one of the three classical DRB1 alleles 15, 16 and 01, and vice versa (Figure 3, Supplementary Table 5). We note that the DRB*01 and *16 alleles are very rare in the Indian data (Supplementary Table 5), so we have limited power to assess their possible role. This correlation was also perfect in the 142 phased Brazil discovery sample chromosomes (Online Methods) where all of the chromosomes carrying the protective haplotype defined by rs9271858 and rs9271255 also carried one of DRB1 *15, *16 or *01 (Supplementary Table 5, Supplementary Figure 2). This lends support to the possibility that the association signal could be driven by classical HLA alleles, and the biological importance of endogenous processing and presentation of antigen from infected macrophages and dendritic cells to CD4 T cells that drive the immune response in visceral leishmaniasis[26,27]. Nevertheless, additional research will be needed to establish whether or not these DRB1 alleles are themselves functionally involved in protection against visceral leishmaniasis, or merely correlated with an as-yet uncharacterised functional variant. Further details of the HLA analysis are presented in the Supplementary.

Figure 3

Schematic representation of the HLA and SNP phased Indian discovery haplotypes in the MHC region. Each row represents an individual haplotype, the three middle columns show different risk models as described in the text, with the shading denoting which class in the model that haplotype is in: column 2 shows the classification of haplotypes based solely on the allele (risk or protective) present at the top SNP rs9271858; column 3 shows the classification of haplotypes based solely on the allele (risk or protective) present at the top SNP in the conditional analysis, rs9271255; column 4 shows the three risk groups as described in the text, which can also be defined as the allele combinations (risk-risk, protective-risk, protective-protective) at SNPs rs9271858 and rs9271255 that most parsimoniously capture their phase-known haplotype associations with disease (see text and Supplementary Figure 1). In India (N=3732), the frequencies of the 2-SNP haplotypes are A_A 0.301, A_G 0.310, G_A 0.372, G_G 0.017. In Brazil (N=1026), the frequencies of the 2-SNP haplotypes are A_A 0.406, A_G 0.097, G_A 0.390, G_G 0.105. The two flanking columns show the HLA alleles at DRB1 and DQB1 and are colour-coded as shown in the legends. The protective risk class correlates perfectly with haplotypes carrying DRB1 *15, *16 and *01. This analysis used the subset of individuals in the Indian discovery data for whom classical HLA alleles were available, and the plot above shows only the 112 haplotypes for which the posterior probability of the estimated phase at both DRB1 and DQB1 was >0.9. Note the haplotype carrying DRB1*0101 is rare and not shown in the plot as it did not occur on a haplotype with a successfully phased DQB1 allele.

Previous studies of HLA associations in human leishmaniasis have provided confusing and often contradictory results generally compromised by small sample size and lack of power (reviewed[9]). One small study did find protection associated with DRB1*15/*16 for visceral leishmaniasis caused by Old World L. infantum[28], the same parasite species that causes visceral leishmaniasis in Brazil[29]. Protection has also been associated with DR2 (=DRB1*15/*16) in small studies of cutaneous and mucocutaneous forms of leishmaniasis caused by L. mexicana[30] and L. braziliensis[31,32], respectively.

There have been relatively few successful GWAS studies of infectious disease susceptibility[33-36]. Our GWAS has identified HLA class II region polymorphisms as significant genetic risk factors for visceral leishmaniasis, with the risk associated with the high-risk haplotype being high by GWAS standards: OR(95%CI) estimated from replication data of 1.50(1.27-1.76). Importantly, genetic variation in the MHC has not been found to be associated with all infectious disease. Our observations focus attention on determining the precise nature of HLA class II variation in determining visceral leishmaniasis susceptibility and offer the potential to contribute to development of novel strategies for disease control.

ONLINE METHODS

Samples

Supplementary Table 1 provides characteristics of samples. In India, discovery and replication sets comprised cases and controls from villages located within a radius of ~100km from the city of Muzaffarpur covering the districts of Muzaffarpur, Vaishali, Samastipur, Saran, Sheohar, East Champaran and Sitamarhi in Bihar State in northeast India. Further epidemiological and demographic details relating to the study site are described elsewhere[37]. An annual incidence rate of 2.49 clinical visceral leishmaniasis cases/1,000 persons has been reported in the region[37]. L. donovani sensu strictu (zymodeme MON-2) was confirmed as the causative agent of visceral leishmaniasis in the study region, in accordance with other reports on clinical isolates from kala-azar patients in the state of Bihar[38-41]. The controls had no history of visceral leishmaniasis or a family history of visceral leishmaniasis among first-, second-, or third-degree relatives. Informed written consent in Hindi was obtained from all participating individuals and from parents of children under 18 years old. Approval for the study was provided by the Ethical Committee of the Institute of Medical Sciences, Banaras Hindu University, Varanasi, India.

From northeast Brazil, 64 families were collected during 1991-1994 as part of the Belém Family Study[6,42]. These families were ascertained from medical records of the Fundacão National de Saude in the States of Para, Maranhão, Piaui based on data from the 1983-85 and 1993-94 epidemics. A further 244 families were collected from study sites near Natal where transmission of the parasite is focal and transient, allowing identification of neighbourhoods with ongoing or recent transmission. These families were ascertained from medical records of the Fundacão National de Saude in Rio Grande do Norte. All families were of equivalent socioeconomic status. The causative agent of visceral leishmaniasis was confirmed as L. infantum chagasi in a subset (~10% over 4 States) of visceral leishmaniasis patients[6]. Further clinical, epidemiological, family structure and demographic information relating to the families are described elsewhere[6,10,21,42,43].

Ethical approval for the Belem Family Study was obtained originally from the local ethics committee at the Instituto Evandro Chagas, Belém, Para, Brazil[6,42]. Approval for continued use of the Belem Family Study samples, and for collection and use of the samples from Natal, has been granted from the local Institutional Review Board at the Universidade Federal do Rio Grande do Norte (CEP-UFRN 94-2004), nationally from the Comissão Nacional de Ética em Pesquisa (CONEP: 11019), and from the Ministerios Cencia e Tecnologia for approval to ship samples out of Brazil (portaria 617; 28 September 2005). Informed written consent for sample collection was obtained from adults, and from parents of children <18 years old. The Brazilian populations studied are long-term (>200 years) admixtures of Caucasian, Negroid and Native Indian ethnic backgrounds, as confirmed in recent analysis of the Natal families[11].

Individuals from both India and Brazil were classified as affected when they had been diagnosed with clinical visceral leishmaniasis that responded positively to specific anti-leishmanial treatment. Data on sub-clinical disease or asymptomatic infections were not included. Diagnosis of clinical visceral leishmaniasis was made on the basis of clinical, parasitological (splenic aspirates in India; bone marrow aspirates in Brazil) and serological criteria as described elsewhere[6,10,21,42,44-46].

DNA sample preparation

For India non-invasive buccal swabs were collected and genomic DNA prepared in the Cambridge laboratory. For Brazil, blood was collected by venepuncture from all available family members. Genomic DNA was prepared directly from blood in Natal[10], or from Epstein-Barr virus-transformed B cells cultured from fresh or cryopreserved peripheral blood mononuclear cells in Cambridge[21]. Genomic DNA for all cases was shipped to the Sanger Institute, Cambridge. Quality was validated using the Sequenom iPLEX assay designed to genotype 4 gender SNPs and 26 SNPs present on the Illumina Beadchips. DNA concentrations were quantified using a PicoGreen assay (Invitrogen) and an aliquot assayed by agarose gel electrophoresis. A DNA sample was considered to pass quality control if the DNA concentration was greater than or equal to 50 ng/μl, the DNA was not degraded, the gender assignment from the iPLEX assay matched that provided in the patient data manifest and genotypes were obtained for at least two thirds of the SNPs on the iPLEX.

Genotyping

Indian discovery cases and controls and the Brazil families were genotyped at the Sanger Institute on the Illumina Infinium platform using the Human660-Quad, a custom chip designed by WTCCC2 and comprising Human550 supplemented with 60,000 additional probes that were intended to allow the genotyping of common CNVs from the Structural Variation Consortium[47]. Replication genotyping for India was carried out at the Sanger Institute using the Illumina Immunochip, a custom chip designed by WTCCC2 and the Immunochip Consortium and comprising 196,524 SNPs. For both chips, bead intensity data was processed and normalized for each sample in BeadStudio; data for successfully genotyped samples was extracted and genotypes called using Illuminus.

Quality controls

SNPs. SNPs were excluded if the Fisher information for the allele frequency was not close to unity (information <0.98) or if the minor allele frequency (MAF) was very low (defined as <0.01%), or if the missingness >0.02, or for extreme departures from Hardy-Weinberg equilibrium (HWE P-value <1×10−20). After applying the above filters in the Indian discovery 526,731 autosomal SNPs remained for further analysis. For the Brazilian families, the above filters left 553,323 autosomal SNPs. Of the two discovery cohorts, 521,134 overlapping SNPs passed QC in both cohorts. Cluster plots of the hit SNPs described in this study are shown in Supplementary Figure 3.

Samples

For quality control of samples typed on the Indian and Brazilian Illumina Human660 quad data, a Bayesian clustering method was used to infer and exclude outlying individuals on the basis of call rate, heterozygosity and signal intensity[12,48]. To remove signal intensity outliers observed for raw intensity data, the difference between the A channel intensity and the B channel intensity was averaged over all SNPs on autosomes for each sample. A similar approach was used taking intensity measures from the A channel on the non-pseudo autosomal X chromosomes to identify outliers and infer gender for all Indian and Brazilian samples. Samples were removed if their inferred gender was discordant with the recorded gender after cross-checking with original database entries, or if less than 90% of the SNPs typed by Sequenom on entry to sample handing (see above) agreed with the genome-wide data. Duplicated individuals were also excluded. For the Indian discovery sample a total of 121/2199 individuals were removed following these quality control checks. For the Indian replication sample 111/2042 individuals were removed. For the Brazilian family study, additional quality control checks were carried out in order to verify (and if necessary adjust) the expected pedigree relationships between the individuals. Genome-wide identity by descent (IBD) sharing probabilities for all pairs of individuals were estimated using a subset of 11,177 autosomal markers using the software PLINK[49]. These estimates were used to assist in pedigree checking and resulted in 49 additional post quality control individuals being removed due to unresolved pedigree relationship discrepancies. In total, for the Brazilian family study, 189/2159 individuals were removed.

Statistical analyses

We carried out association analyses using a novel[13] variance components method (a linear mixed model, similar to Kang et al.[17]). The linear mixed model explicitly accounts for correlations in individuals’ phenotypes due to their relatedness through specification of the phenotypic variance/covariance matrix in terms of parameters representing underlying genetic and environmental components of variance (estimated during the model fitting procedure) and pairwise relatedness measures (kinship coefficients), estimated prior to model fitting on the basis of the genome-wide genotype data. In this analysis we used a standard linear mixed model where =(y1,…,y) is a vector of responses (coded 1/0 for case/control) on n subjects, =(x) is the n × K matrix of predictor values for any variables to be modelled as fixed effects (including variables representing covariates and genotypes at any SNPs currently under test), =(β1, … β) are regression coefficients (to be estimated) representing the linear effects of the predictors on the response, and Q and ε are random effects assigned the distributions and ε~N(0, σI), where and are parameters (to be estimated) representing the genetic and environmental components of variance respectively, I is the n × n identity matrix and Φ is the n × n matrix of pairwise kinship coefficients.

Although use of this linear mixed model was originally proposed for pedigrees with known relationships[50-53], recently this approach has gained popularity for use with samples of unknown or uncertain relationship[13,17,18,54], including apparently unrelated samples who may nevertheless display distant levels of common ancestry. For this purpose, the pairwise kinship coefficients modelling either close or distant relatedness are estimated (prior to fitting the linear mixed model) on the basis of genome-wide genotype data, rather than being fixed at known theoretical values.

Computational considerations have led to the development of several faster approximations for constructing tests of the fixed effects of interest in the above linear mixed model[15-17,53,55,56]. These approximate tests have been implemented in various software packages including EMMAX, TASSEL, FaST-LMM and GenABEL. For the analyses presented here, we used our own implementation of the linear mixed model developed for a previous study[13] that has the advantage, in common with the recently-developed GEMMA package[57] of fitting the full (rather than an approximate) model, which in principal can lead to a small increase in power, depending on the underlying true model. Our implementation also allows a transformation between the parameters of the linear model shown above and the logistic model (which is the usual model for case/control data), allowing the convenient estimation of effects on the log odds scale and the generation of resulting odds ratio estimates. We confirmed that the Indian discovery results did not change significantly after the inclusion of the first 10 PCs in the linear mixed model (data not shown).

Combined P-values across Indian (discovery and replication) and Brazilian datasets were calculated in R using an inverse variance fixed effects meta-analysis.

In addition, for the Brazilian family study we verified our results using two complementary alternative methods implemented in the software packages FBAT[58] and ROADTRIPS[59]. Results were comparable across the different methods (data not shown).

HLA typing

The HLA DRB1 and DQB1 typing of 100 selected Indian discovery individuals was performed using a sequencing-based method[60]. To select these individuals, SNPs in the class II region DRA-QB2 (32,514,320 to 32,840,188 NCBI Build 36) were first pruned in PLINK[49] to generate a set of 78 SNPs in approximate linkage equilibrium (r2<0.8) and with MAF>0.1. LD blocks across these 78 SNPs were then visualised using Haploview[61] and the top 59 haplotype tagging SNPs were selected and used to generate haplotypes for all 1866 unrelated individuals in fastPHASE[62]. The subset of 100 individuals selected for sequence-based HLA typing represented the 135 most common 59-SNP haplotypes in the region. Automated sequencing was carried out on ABI Prism 3730 or 3730xl Genetic Analysers (Applied Biosystems; Foster City, CA, USA) and HLA-DRB1 analysis was carried out using ASSIGNV4.0.1.36 (Conexio Genomics; Applecross, WA, Australia). Where possible, individuals were assigned DRB1 or DQB1 specificities to the 4-digit level[63], i.e. to a specific HLA protein encoded on each chromosome (for example, DRB1*1501, where *1501 is a specific functional DRB1 protein). Where this was not possible, individuals were assigned to a 2-digit level allele group, e.g. HLA-DRB1*15 indicates an allele group for which all specific functional DRB1 proteins share a common ancestry.

HLA-DRB1 and DQB1 typing was further performed in a subset of 71 unrelated individuals from the Brazilian families using the same sequence-based method as described above for the Indian samples.

PHASING the SNP and HLA data

SNPs in the MHC region and HLA-DRB1 and -DQB1 alleles were phased using IMPUTE2[64] and PHASE[65]. The Indian discovery and Brazilian replication cohorts were treated separately throughout the phasing pipeline.

For the Indian discovery data, IMPUTE2[64] was used to phase all the SNPs on chromosome 6 for 1866 individuals, using HapMap3 populations as the haploid reference panel. For the Brazilian data 5,667 SNPs in the MHC (20,733,613- 41,984,313bp) were phased. The phasing of the Brazilian data was performed in two stages, firstly with an unrelated set of 498 parents and secondly an unrelated set of 268 offspring, totalling 766 individuals. The trio relationships were then used to check the accuracy of both SNP and HLA allele phasing.

Treating the IMPUTE2 SNP phasing as fixed, PHASE[65] was then used to phase the HLA alleles typed at HLA-DRB1 and HLA-DQB1 onto HLA region SNP haplotypes, selecting only haplotypes phased at probability >0.9. PHASE[65] was used to phase the multi-allelic HLA-DRB1 and HLA-DQB1 alleles onto the known SNP haplotypes, specifying known phase at all of the SNPs and the PHASE options –MR and -d1. The HLA alleles were given positions: - DRB1 was set at 32,597,662bp, between rs28756238 and rs41546317 and HLA-DQB1 was set at 32,735,635bp between rs28724231 and rs1063355. Both HLA-DRB1 and HLA-DQB1 alleles were phased in the Indian data, and HLA-DRB1 was phased in the Brazilian data.

After QC of the PHASE data, in the Indian discovery set there was a total of 3,732 SNP haplotypes, 154 with HLA-DRB1 alleles and 153 with HLA-DQB1 alleles. Of these, 112 were phased to 4-digit level for both DRB1 and DQB1. In the Brazilian dataset there was a total of 1,481 SNP haplotypes, 142 of which were phased with 2 digit HLA-DRB1 alleles.

Haplotype analyses

Haplotype associations were investigated by applying two methods; a conditional analysis (see Results) and a Bayesian approach. Both methods reached the same result, each fitting a three-haplotype model of risk.

GENECLUSTER, which adopts a Bayesian approach[66] was used to look for primary and secondary association signals at known and putative SNPs. This method analyses the genealogy of the case-control sample to find evidence for causal mutation(s). It uses the genealogy of a reference haplotype panel (in this case the HapMap3 GIH haplotypes) to approximate the genealogy of the case-control sample at positions across the region by clustering the case-control haplotypes under the leaves of the reference genealogy. Placing a disease mutation on a branch of the reference genealogy segregates the haplotypes at the leaves into two different risk groups: those case-control haplotypes that fall under the disease mutation (and therefore carries the mutation) and those that do not. For a given mutation in the tree, GENECLUSTER carries out a Bayesian test of association. The model can also be extended to two disease mutations in the genealogy. For this data, the genealogical tree was estimated on chromosome 6 at 32,678,817 bp. The log10(Bayes Factor) for the 2-mutation (three haplotype) model was 7.36, compared to log10(Bayes Factor) 5.42 for a single mutation (two haplotype) model.

Following the linear mixed model conditional analysis, we phased the SNP data for all individuals in the India discovery data and found that the two SNPs associated in the conditional analysis (rs9271858 and rs9271255) tag the haplotypes from the GENECLUSTER analysis. The first mutation perfectly correlated with the protective haplotype and the second mutation separated the risk halotype into two, which are well tagged by the two SNP intermediate and risk haplotypes (r2 protective=1, intermediate=0.49, risk=0.60) (overall correlation 0.86).