Introduction

Hundreds of common variants have been established for immune-mediated diseases, mainly from genome-wide association studies (GWAS) in populations of European ancestry.1 A large proportion of the European association signals replicate within non-European populations. Yet, additional risk variants can be found using multi-ethnic cohorts.2, 3, 4, 5

Coeliac disease (CeD) is a common, autoimmune disorder, with a similar prevalence among Europeans and north Indians (1–3%).6 The largest genetic risk to CeD is conferred by the HLA haplotypes, which account for approximately 35–40% of the risk.7 Recent progress in understanding the genetic architecture of CeD has identified 39 non-HLA loci.8, 9, 10, 11 Although these analyses provided much insight into genetic risk factors for CeD, they were largely based on European populations, so that the risk of these associations in other non-European populations needed to be assessed.

We explored the contribution of reported risk variants and searched for new associated variants using a cohort of 1233 north Indian individuals. We acknowledge the fact that about half of this cohort contributed to the previously published Immunochip study. However, owing to the large size of the European cohorts, it is likely that any effect derived from the north Indian population would have been overwhelmed by the European effects (620 north Indian samples, which is only 2.5% of the total 24 269 samples). In this study, we have doubled the sample size of the north Indian cohort and analysed the transferability of known European variants to this population. We also took advantage of the higher degree of consanguinity in the north Indian population12 and assessed the excess homozygosity to explore the possibility of recessively acting risk variants.

Materials and methods

Ethical statement

This study was approved by the respective institutional and university ethical committees. Informed written consent was received from all participants.

Study populations

We used two ethnically distinct cohorts from northern India and the Netherlands. The Indian cases (n=497) and controls (n=736) were recruited from Dayanand Medical College and Hospital, Ludhiana, in Punjab (northern India). The Indian controls included blood donors who tested negative for CeD serology. All the Indian subjects had self-reported northern Indian ethnicity. The Dutch cases (n=1150) and controls (n=1173) were recruited from three university medical centres (Utrecht, Leiden and VUmc Amsterdam) in the Netherlands. A small proportion of cases were recruited via the Dutch CeD patients’ society. In both cohorts, Indian and Dutch CeD patients were diagnosed according to standard clinical, serological and histopathological criteria, following the ESPGHAN criteria.13 Most of these samples have been described elsewhere.11

DNA extraction and genotyping

Most of the DNA samples were derived from blood, or from saliva for a small proportion of the Dutch cases and controls. Samples were hybridized on the Immunochip platform, a custom-made chip with 196 524 markers.11 Genotyping was carried out according to Illumina’s protocol at the genotyping facility of the University Medical Centre Groningen. The variant calls were exported from Genome Studio with human genome GRCh36/hg18 coordinates and further converted into GRCh37/hg19 coordinates using the UCSC liftOver tool.14

Genotype data quality control

We manually inspected cluster plots and adjusted them if required for about 150 000 markers. For about 20 000 SNPs, the assay performed poorly, and those variants were marked as uncalled and excluded from the final analysis. This left 178 111 high-quality SNPs that were exported from Illumina’s Genome Studio. We required samples to have ≥98% call rate. Individuals showing a high degree of genetic relatedness (estimated kinship coefficient >0.2) or with a discordant sex were removed. Population outliers were identified by principal component analysis (PCA) implemented in EIGENSOFT 5.0.1.15, 16 We excluded markers with a call rate <99% or with significant deviation from Hardy–Weinberg equilibrium (P>1 × 10−3) and were left with 176 799 variants. A large number of Immunochip variants was of low frequency, but our study was underpowered to assess their association in the north Indian cohort; we therefore filtered out SNPs with a minor allele frequency (MAF) <10%. This excluded 88 083 and 87 442 variants in the north Indian and Dutch cohorts, respectively. We further only used the set of variants shared across the two cohorts, resulting in 88 716 SNPs in the final analysis. Genomic inflation for each cohort was estimated using 3016 ‘neutral’ variants present on the Immunochip array for the replication of the GWAS for reading and writing skills and therefore unlikely to be confounded by the immune signals.

Statistical analysis

Chip-wide analysis

For both populations, we applied additive and recessive models using logistic regression, including the first four principal components and gender as covariates. We used an Immunochip-wide P-value cutoff of P=3.5 × 10−6 calculated as 0.05 divided by 14 136 linkage disequilibrium (LD)-independent SNPs (r2<0.1 computed across 100 SNPs, with a 10-SNP sliding window) with MAF>10% and based on the north Indian cohort to determine a significant association. We used P<1 × 10−4 to report a suggestive association. To reassess the significance of the associations detected in this test, we also performed a minimum of 100 000 permutations to compute empirical P-values.

Replication

We used two approaches to test for replication of the 39 non-HLA CeD loci that have been reported in Europeans11 in the north Indian cohort: the exact SNP replication, and the locus transferability.

  • Exact SNP replication: the reported index SNPs (representing primary as well as secondary, independent signals) were tested for association in the north Indian cohort. We considered replication to be significant at P<0.05. During the quality control, we removed 13 of the total of 57 independent SNPs associated to the 39 loci because of MAF<10%. The exact SNP replication was applied to 44 reported SNPs, representing 38 distinct loci.

  • Locus transferability test: we performed a transferability test to account for LD differences across different ethnic populations. If the causal SNP is poorly tagged by the associated SNP reported in the discovery population (European), the differences in LD between the north Indian and Dutch populations could result in either a lack of replication or a different localization of the association signal. Published LD boundaries were used to define each of the 39 non-HLA loci.11 First, we identified all the LD-correlated markers for each locus, measured by r2≥0.05 between the index SNP and all the variants present in CeD loci using the Dutch controls. We then checked the association status of these variants in the north Indian Immunochip data set for transferability. SNP rs13132308 from the IL2/IL21 locus was not present in the north Indian data set and was replaced by its best proxy SNP rs13151961 (r2>0.9 based on the Dutch controls). The ELMO1 locus and SCHIP1/IL12A locus were excluded because the index variants for these loci, or their proxies failed the 10% MAF filter. We required at least one significant variant per locus. The significance threshold was assessed in two ways. First, across all tested loci (37) we identified 142 LD-independent SNPs (pairwise r2<0.1) and used P=0.05/142=3.52 × 10−4 as the significance threshold. Second, we used a permutation-based test as an alternative way to account for LD structure. We performed 1000 permutations for which we randomized phenotype labels. We identified a nominal P-value for each of the 37 loci in each of the permutation results. We then defined the 5th percentile of the nominal P-values as a significance threshold.

Runs of homozygosity (ROH)

The north Indian population has a reported higher degree of consanguinity,12, 17, 18, 19 therefore we sought to identify if there were significantly associated runs of homozygous variants. We first pruned both data sets for LD using an r2≥0.8 as a threshold to avoid detecting false-positive regions because of homozygous stretches of LD-correlated SNPs.20, 21, 22 This is particularly likely to occur for the Immunochip, as it includes regions with a high density of markers that are often in strong LD with each other. Then, we identified regions with ROH, which were defined as at least 50 contiguous homozygous SNPs with no interruption (without heterozygous positions in between) and with less than 500-kb distance between neighboring SNPs. Association of ROH regions was performed using PLINK v1.07.23 We considered a locus significant after Bonferroni correction using the total number of ROH found per data set (8) in our analysis with a P<0.006. We assessed the false discovery rate (FDR) by running 10 000 permutations in which we randomly shuffled the phenotypes and reassessed the significance of association.

Data availability

Genotype data for the north Indian cohort are available through the European Genome-phenome Archive (https://www.ebi.ac.uk/ega/) under the accession number EGAS00001000849. The summary statistics can be obtained by contacting the authors directly.

Results

We obtained high-quality genotype data for 88 716 common SNPs that were shared between the north Indian and Dutch cohorts. PCA revealed different clustering patterns between the two cohorts, consistent with their different ethnicities. Furthermore, the case–control samples were homogeneous within each of the data sets (Supplementary Figure 1). We used the first four principal components to control for possible population substructure and did not observe significant genomic inflation in any of the data sets (λIndia=1.036, λNetherlands=1.039).

The strongest association in the north Indian cohort was present at the HLA locus. Comparable to the Dutch sample, the SNP rs2854275 was the top signal (pIndia=8.17 × 10−49, ORIndia=6.97 (5.38–9.04)); pDutch=4.413 × 10−121, ORDutch=10.08 (8.3–12.23)). This replicated the known, very strong effect of HLA in CeD. No other locus reached genome-wide significance (Figure 1 and Table 3). The rs4948256 at 10q21.2 in the ANK3 gene reached our Immunochip-wide significance threshold (P=9.28 × 10−7, ppermuted=1 × 10−6). The rs4758538 at 11p15.4 in the OSBPL5 gene (P=8.57 × 10−5, ppermuted=7.93 × 10−5) and rs17080877 at 18q22.2 near the DOK6 gene (P=2.68 × 10−5, ppermuted=1.7 × 10−5) both were suggestively associated. None of these associations could be replicated in the Dutch cohort.

Figure 1
figure 1

Additive and recessive association results in the north Indian population. Manhattan plot of association under (a) an additive and (b) a recessive model for the north Indian (blue dots) and Dutch populations (black dots). The Immunochip-wide significance cutoff (P=3.5 × 10−6) is depicted as a green dotted line and suggestive significance threshold as a red dotted line (P=1.0 × 10−4). We excluded the HLA region.

Full size image

Replication of European associations in the north Indian cohort using exact and transferability tests

The exact association in the north Indian cohort was observed at five previously reported loci, corresponding to FASLG/TNFRSF18 (rs12142280, P=0.002, OR=0.69 (0.54–0.87)), ICOSLG (rs58911644, P=0.004, OR=0.68 (0.51–0.88)), PFKFB3/PRKCQ (rs2387397, P=0.02, OR=0.77 (0.61–0.96)), SCHIP1/IL12A (rs2561288, P=0.03, OR=1.2 (1.01–1.44)) and ZMIZ1 (rs1250552, P=0.03, OR=1.2 (1.01–1.44)) (Table 1). At all these variants, the directions of association were the same as those previously reported (Supplementary Table 1). Three of these five loci, FASLG/TNFRSF18, ZMIZ1 and ICOSLG, previously reached P<0.05 in the north Indian samples used in the initial study by Trynka et al11 (Supplementary Table 1).

Table 1 Five loci that were significantly (P<0.05) replicated in the north Indian cohort in the exact SNP test
Full size table

Using the transferability analysis of 39 loci, we observed two associations, which passed two independent significance thresholds (Table 2). One was at the PFKFB3/PRKCQ locus represented by the SNP rs744254 (P=2.8 × 10−4, OR=0.70 (0.57–0.84)); the other was at the PTPRK/THEMIS locus, rs4142030 (P=3.4 × 10−4, OR=1.39 (1.16–1.66)) (Supplementary Table 2).

Table 2 Two loci significantly replicated in the north Indian cohort in the transferability test
Full size table

Recessive analysis in the north Indian population

The ROH analysis identified 16 genomic regions with >50 consecutive homozygous SNPs each; two regions showed suggestive association with CeD in the north Indian cohort. They were located on chromosomes 2q11.2 (P=4.1 × 10−3, FDR=0.11) and 6p21.33 (P=1.4 × 10−4, FDR=0.0001) (Supplementary Table 3). None of those regions replicated in the Dutch cohort.

Complementary to the ROH analysis, we also applied a recessive model to test for association across all 88 716 SNPs. Apart from the SNP rs9271850 (P=3.67 × 10−23, OR=0.12 (0.08–0.18)) located in the HLA locus, we observed four variants (MAF>0.25) at four loci that reached the threshold of suggestive association (P<1 × 10−4; Figure 1 and Table 3). The rs744253 at 10p15.1 (P=9.27 × 10−5, OR=0.3859; ppermuted= 6.7 × 10−5) corresponds to PFKFB3/PRKCQ locus that was replicated in the direct association and transferability tests described above. The remaining three loci were new suggestive associations, with the strongest observed at the intergenic SNP rs963872 (P=1.2 × 10−5, OR=2.01 (1.47–2.75); ppermuted=1.2 × 10−5) at 5p15.33, near IRX2 gene, followed by rs10994257 (P=1.8 × 10−5, OR=2.5 (1.67–3.99); ppermuted=5 × 10−6) at 10q21.2 mapping to the ANK3 gene and rs6010998 (P=6.17 × 10−5, OR=2.8 (1.67–4.66); ppermuted= 2.5 × 10−5) at 20q13.33 in the RTEL1 gene.

Table 3 Association results for Immunochip-wide additive and recessive models in the north Indian population compared with the same model applied to the Dutch cohort
Full size table

Discussion

Recent trans-ethnic studies have indicated that a large proportion of common disease- or trait-associated variants established in populations of European descent also contribute to the risk in other ethnic groups.4, 24, 25 Yet this observation cannot be generalized because only a few of the GWAS were conducted in non-European populations.26 Novel loci have been reported in these limited studies.27, 28 Thus, more multi-ethnic studies may help identify not only shared disease determinants, but also additional risk variants, which may provide leads toward a better understanding of the disease biology.29

To our knowledge, there have been very few genetic association studies for CeD performed in non-European populations, and those that have been published have focused only on the HLA effects based on very small sample size.30, 31, 32, 33, 34, 35, 36 In this study, we attempted to validate known genetic determinants of CeD and identify potentially novel loci in the north Indian population. We used the Immunochip genotyping platform for this study because it offers the highest coverage of regions associated to immune-related disorders, including those conferring risk to CeD.

Given our modest sample size and the limited power of our study, it is likely that the real replication rate of the European associations in the north Indian population is much greater than we can report here. In this context, it should be mentioned that although half of the samples of the north Indian cohort were previously analysed in the study by Trynka et al,11 the study focused strongly on the effects in samples of European ancestry (97% of the 24 269 samples were of European descent). In this study, we wanted to focus specifically on the role of the European variants in the north Indian population and to explore potential new associations that could have been diluted in the previous study because of the dominance of the European cohorts. Our study suffers from a small sample size and lack of replication, for these reasons we acknowledge that the results presented here need to be validated further in future replication studies.

We successfully replicated HLA and five of the 38 tested European SNPs (13%); these included the loci harbouring the FASLG/TNFSF18, SCHIP1/IL12A, PFKFB/PRKCQ, ZMIZ1 and ICOSLG genes. The transferability test further confirmed association to the PFKFB/PRKCQ locus and identified an additional association to a variant at the PTPRK/THEMIS locus with suggestive evidence of association. The top association at this locus, rs4142030, has a very low correlation with the European index variants (r2=0.083, D′=0.53) and is located in intron 2 of THEMIS, in contrast to the European-associated SNP, which is located in intron 30 of PTPRK. Given that our study was underpowered, we acknowledge that these suggestive associations must be further replicated in north Indian populations. However, it is tempting to speculate that this localization of the association signal could result from the independent effects of each of the genes in the different ethnicities. Both genes are of functional relevance to CeD aetiology and a recent study showed that the mRNA levels between these two genes have higher expression levels in CeD patients but not in controls.37

The analysis of the recessive effects implicated multiple genomic regions at a suggestive significance. With the run of homozygosity test, we found a suggestive association to a region at 2q11.2 that spans 1.14 Mb and includes the LYG1, TXNDC9, EIF5B, REV1, AFF3 and LONRF2 genes. This region has been associated with multiple immune-related phenotypes, including type 1 diabetes (T1D) and rheumatoid arthritis (RA).38, 39, 40 The second region maps to the HLA locus at 6p21 and contains the MICA gene, which is also a locus with a pleiotropic effect in phenotypes such as Behcet’s disease, RA or HCV-induced hepatocellular carcinoma.41, 42, 43, 44, 45, 46, 47 The Immunochip-wide analysis under a recessive model suggested four loci that have also been linked to other phenotypes, for example, ANK3, which has been associated with psychiatric disorders,48, 49 and with other loci that include PFKFB3/PRKCQ and RTEL1/TNFRSF6B reported to be associated with immune-related diseases such as T1D, RA and irritable bowel disease.39, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63 Apart from HLA, the ANK3 locus was the only one that passed Immunochip-wide significance threshold when we considered the additive model. None of these loci were replicated in the Dutch cohort, which could reflect a population-specific effect, but they require further replication in the north Indian population to robustly establish their association with a genome-wide significance. Although the genes mapping to the associated loci that we report here appear to be relevant candidates for CeD pathogenesis, our study was biased toward identifying biologically relevant genes because the Immunochip was designed to specifically target regions previously reported as associated to immune-related diseases.