A Possible Mechanism behind Autoimmune Disorders Discovered By Genome-Wide Linkage and Association Analysis in Celiac Disease

Celiac disease is a common autoimmune disorder characterized by an intestinal inflammation triggered by gluten, a storage protein found in wheat, rye and barley. Similar to other autoimmune diseases such as type 1 diabetes, psoriasis and rheumatoid arthritis, celiac disease is the result of an immune response to self-antigens leading to tissue destruction and production of autoantibodies. Common diseases like celiac disease have a complex pattern of inheritance with inputs from both environmental as well as additive and non-additive genetic factors. In the past few years, Genome Wide Association Studies (GWAS) have been successful in finding genetic risk variants behind many common diseases and traits. To complement and add to the previous findings, we performed a GWAS including 206 trios from 97 nuclear Swedish and Norwegian families affected with celiac disease. By stratifying for HLA-DQ, we identified a new genome-wide significant risk locus covering the DUSP10 gene. To further investigate the associations from the GWAS we performed pathway analyses and two-locus interaction analyses. These analyses showed an over-representation of genes involved in type 2 diabetes and identified a set of candidate mechanisms and genes of which some were selected for mRNA expression analysis using small intestinal biopsies from 98 patients. Several genes were expressed differently in the small intestinal mucosa from patients with celiac autoimmunity compared to intestinal mucosa from control patients. From top-scoring regions we identified susceptibility genes in several categories: 1) polarity and epithelial cell functionality; 2) intestinal smooth muscle; 3) growth and energy homeostasis, including proline and glutamine metabolism; and finally 4) innate and adaptive immune system. These genes and pathways, including specific functions of DUSP10, together reveal a new potential biological mechanism that could influence the genesis of celiac disease, and possibly also other chronic disorders with an inflammatory component.


Introduction
Celiac disease (CD) is a common chronic disease and even though most often diagnosed in early childhood, it can present itself at any age. Most of the individuals with CD remain undiagnosed and an estimated 2% of the Swedish population is affected without having been diagnosed [1]. Ongoing disease will increase the overall risk for developing other chronic inflammatory diseases, neurological manifestations and malnutrition disorders. CD is the only autoimmune disorder where the actual genes responsible for the association in HLA are known (HLA-DQA1 and HLA-DQB1) [2]. In the past few years Genome Wide Association Studies (GWAS) have had tremendous success in identifying new genes, or gene regions, that influence common diseases. These studies use several hundreds of thousands of genetic markers (single nucleotide polymorphisms, SNPs) across all human chromosomes in order to pin down the chromosomal locations of genes, which could influence the disease.
A large joint effort has been done, not the least in CD, and 40 new CD-associated genetic regions marked by SNPs have been discovered [3][4][5][6][7]. However, these genes cannot account for all CD heritability, and part of the genetic variance that influences disease development is still unknown [8].
Most GWAS so far have been performed on case control samples. A case control study design has some advantages compared to using a family study design. For example, in a case control design it is possible to select a perfectly matched set of controls to increase the chance of discovering susceptibility genes, and furthermore, cases and controls are usually easier to collect than individuals from the same family. However, using a family material can be a very good complement to a case control design. First of all, families with several affected members are likely to have a stronger genetic component compared to sporadic cases. Familial cases tend to be enriched for disease-predisposing alleles and there is an increased power especially for detecting rare genetic variants [9]. Another important fact is that statistical analyses based on family data are robust against population stratification. Already in their paper from 1996, Risch and Merikangas suggested that all sib-pair families collected for nonparametric linkage analysis in complex diseases, should be re-run ''Genome-Wide'' using SNP markers and the potentially more powerful Transmission Disequilibrium Test (TDT) [10]. The TDT test in sib-pairs is a test of linkage in the presence of association. Hereafter we refer to whole genome sibling TDT as ''Linkage GWAS''.
In this study, we aimed to uncover additional genetic factors in CD by performing a Linkage GWAS using 206 affected children (sib-pairs) within 97 nuclear families using the TDT test. In addition to the Linkage GWAS we explored gene-gene interactions and pathway analyses. We also performed a non-parametric linkage (NPL) analysis and compared the results with the published linkage analysis, with microsatellite markers, performed in the same set of families previously [11]. Furthermore, quantitative PCR was used to investigate levels of gene expression in small intestinal biopsies from additional patients with CD autoimmunity (CAI) and control patients. Finally, we stratified the TDT analysis on HLA genotype. It has been shown that carrying DQB1*02 on both chromosomes (i.e. being homozygous), confers higher risk of developing CD as compared to heterozygote individuals [12]. It is therefore conceivable that heterozygote individuals may require more additional risk factors outside HLA, in order to accumulate sufficient risk to develop CD, compared with homozygous individuals. Based on this assumption we stratified the patient in an HLA low-risk group and an HLA high-risk group. By stratifying the Linkage GWAS, we expected to uncover even more of the so-called ''missing heritability'' in CD. This strategy could identify different risk factors all together or perhaps a more likely scenario is that the same risk factors outside HLA would just be more common in the HLA low-risk group.

Genotyping and Imputation
We included single nucleotide polymorphism (SNP) markers that had a call rate above 97%, which led to the exclusion of 1.3% of the Omni Express and 0.6% of the 660W-Quad SNP markers. Out of the 127,535,126 imputed genotypes, 88.3% had a posterior probability of over 0.95. Approximately 90% of the 944,512 SNP markers had a minor allele frequency of at least 0.01 after imputation.

Transmission Disequilibrium Test (TDT)
All markers from the TDT analysis are shown in Figure 1a. As expected, the region around the CD associated HLA genes on chromosome 6 showed the strongest association with the most significant p-value reaching 4.9610 221 at marker rs424232. In Table 1, we present the 35 most significant associations found outside of HLA (HLA defined as SNP markers located within 27-34 Mb on chromosome 6). The most significant finding outside of the HLA region was the marker rs12734338 on chromosome 1, including the PPP1R12B gene.

HLA Stratified Transmission Disequilibrium Test (TDT)
In Figure 1b and Table 2, we present results from the TDT analysis stratified on the HLA-DQ risk factor. For this analysis 115 affected offspring trios were included in the ''low-risk'' group and 88 trios were put in the ''high-risk'' group. A region including the DUSP10 gene (also known as MKP5) reached genome-wide significance (p-value = 3.8610 28 ) in the low-risk group. Figure 1c presents this region including the most associated SNPs plotted on the x-axis using SNAP.

Interaction Analyses
Since some markers just below genome-wide significance are still expected to be true findings, we wanted to try and separate these from the, in fact, true negative findings (those that show linkage and association close to genome-wide significance just by chance). In total, 603 SNP markers from 383 independent regions and their surrounding genes were identified by three inclusion criteria ( Fig. 2 and Table S1). These genes were subsequently used for pathway and two-locus interaction analyses.
Two-locus interaction analysis. Two-locus interaction analysis, identified 582 SNP pairs with a p-value of less than 1.0610 24 for the test comparing the model M 0 of no association and the general two-locus model M G . Out of these, 101 pairs from 87 regions deviated significantly (p,0.05) from a purely multiplicative model (M M ), which is the best fitting model when at least one of the SNP markers is false. Under the null hypothesis we expect to find 29 such pairs. The 101 pairs showed either epistasis (individuals carry both risk alleles) or evidence of heterogeneity (individuals carry either the one or the other risk allele from the two loci).
The results with a p-value ,1.0610 24 for epistasis and those with high p-value (.0.05), which represent pairs that did not show convincing deviation from the heterogeneity model are listed in Table 3 and 4. Several loci were in an epistatic relationship with HLA; rs4899272 (ACTN1), rs1073933 (COX7C), rs10482751 (TGFB2), rs571879 (APPL1) and rs7590305 (FABP1). Also, previously identified susceptibility loci for CD were involved in Figure 1. Manhattanplot of the TDT p-values. a) The location of all genotyped SNPs on chromosomes 1-22 and X plotted on the x-axis. -log10(p-value) result for each SNP and all transmissions on the y-axis. b) The location of all genotyped SNPs on chromosomes 1-22 and X plotted on the x-axis. -log10(p-value) result for each SNP and all transmissions, to children in the low risk group, on the y-axis. c) Regional plot of association results and recombination rates, within the region surrounding DUSP10, generated by SNAP (http://www.broadinstitute.org/mpg/snap/ldplot.php). The x-axis show 500 kb around the most associated SNP. Genomic locations of genes within the region of interest (NCBI Build 36 human assembly) were annotated from the UCSC Genome Browser (arrows  [13] are shown in Table 5, 6 and 7. Several clusters were significant after correction for multiple comparisons. The most significant network implicated by IPA included DUSP10 ( Fig. 3 and Table 8). The second top network included the MHC complex (HLA) and the third top network included LPP, which is located within the most significantly non-HLA associated region identified in CD so far [3].

Gene Expression
Out of the 34 selected target genes, three were from the top associated SNPs (DUSP10, SVIL and PPP1R12B) and the remaining were genes identified from the two-locus and pathway analysis. Eight genes showed significant up-or down-regulation after correction for multiple testing using Bonferroni correction (Fig. 4). For the top associated genes, several transcript variants were tested (Table 9). For the PPP1R12B gene, Isoform c and d (transcript variants NM032103.2 and NM032104.2) also known as the small subunit (sm-M20) of myosin light chain phosphatase, show significant up-regulation in patients with CD autoimmunity compared to control patients. An additional ten genes showed nominally significant differences in expression (Table 9).

Non-parametric Linkage (NPL)
The strongest linkage outside of HLA was detected in chromosome regions 5q23.2-q33.1, and 1q32.1. In total, thirteen regions with an NPL point wise p-value below 0.01 were detected ( Fig. 5 and Table 10). In our previous linkage-scan, using almost the same set of families, we detected only one region (11q23-25) with a point wise p-value below 0.01 [14]. The reason for the improved results is mainly the almost perfect information content achieved by a dense set of highly successful SNP markers compared to a relatively sparse set of less successful microsatellite markers. Also in the NPL analysis, the PPP1R12B gene was located in one of the top regions (1q32.1).

Discussion
This study confirmed some previous GWAS findings and in addition, it established a new genome-wide significant region containing the DUSP10 gene. The top markers, rs12144971 and rs4240931 showed a substantial effect size in the HLA low-risk group with a transmitted versus non-transmitted allele ratio of 3.11 (Table 2).

DUSP10, TNF-a and Tissue Transglutaminase (TGM2)
The protein product of DUSP10 preferentially binds to the stress-activated p38 MAPK (mitogen-activated protein kinase) and plays an important role in regulating chemokine induction after infection by various pathogens [15], and in coordinating MAPK activity in response to oxidative stress [16]. In previous studies, both p38 MAPK and DUSP10 have been shown to activate TNFa [17,18], of which one also demonstrates that TNF-a upregulates TGM2 (the gene encoding the main autoantigen in CD [19]) in intestinal mucosa from untreated CD patients [17]. Whether this up-regulation of TGM2 is of importance for the immune response leading to formation of IgA-tTG and IgG-tTG autoantibodies, the serological markers for CD is still unresolved.

Pathway Analyses
In order to discover possible functional connections between DUSP10 and other genes, we analyzed genes surrounding the top 603 markers. A total of 845 genes were used in the analysis. Ingenuity pathway analysis (IPA) included DUSP10 within the most significant network. Also part of this network were GLS and RGS1, two genes previously identified within significant GWAS loci [3], as well as the insulin (INS) gene, and the immune regulatory nuclear factor kappa B (NF-Kb) complex ( Fig. 3 and Table 8). The second top network included the MHC complex (HLA) and also several genes within already identified GWAS loci: ACTN1, CD247, CCR5, ICOS and STAT1 [3]. In addition, both IPA and GeneTrail [13] identified T2D genes as the most significantly overrepresented gene cluster after correction for multiple testing (Table 5 and 6). Among this set of genes surrounding the 603 markers, many genes belonged to growth and nutrient signaling pathways, for example, INS, INSR, EGF, POMC, TIPRL and PRR5L. There were also related genes directly involved in energy metabolism; PDK1, COX7C, COQ3 and GLS.

Overlapping Results with Other GWAS Findings
Surprisingly, four out of six top loci identified by a GWAS for anorexia nervosa [20] and two out of three loci involved in plasma glucose levels in type 1 diabetic patients [21] were among our 603 and 35 best SNP markers respectively. One of the genes in anorexia, namely AKAP6, is also associated to fasting insulin-related traits as well as the autoimmune disease Ankylosing spondylitis [22]. Of the 40 identified regions in CD, seven regions overlap with our 603 SNP list (LPP, STAT4/GLS, RGS1, CCR1/CCR3, PUS10, ICOS/CTLA4 and CD247). Out of the 69 regions reported in the GWAS catalog for type 1 diabetes, eight overlap with the regions reported in this study and out of those eight, CTLA4/ICOS also overlap with the previously reported CD associations.
We compared minor allele frequencies between the previous CD GWAS by Dubois et al. and our GWAS. In their top 42 associations, there was no SNP below a minor allele frequency of 0.08. In our top 42 associations, we identified five SNPs with a minor allele frequency below 0.06. This observation could just be a chance finding or perhaps an indication that rare variants are easier to discover using families. We also identified a relatively rare variant in the LPP gene region (rs17283813), with a minor allele frequency of 0.075. This SNP was not at all significant in the GWAS by Dubois et al. (Table S1).
Neither was there an association with the DUSP10 region in the GWAS by Dubois and co-workers. The associated markers in the DUSP10 region in our GWAS have a minor allele frequency around 0.5 and are hence very common in the population. It is difficult to say if this is a population specific effect or if DUSP10 could be detected in an HLA stratified population from another ethnicity. Interestingly, the DUSP10 region has also been identified as a risk factor for colon cancer by a meta-analysis of three GWAS from the UK. This is an indication that colon cancer and CD could share genetic risk factors.

Key Metabolic Regulators as well as the Top Associated gene PPP1R12B were Differently Expressed in CD Cases Compared to Controls
Another important finding was the difference between cases and controls and their gene expression patterns in the small intestine. Eight of the 34 candidate genes selected for quantitative measurements of gene expression, including PPP1R12B, PDK1, GLS, PRR5L and the INSR, showed significant up or down regulation of mRNA levels in cases compared to controls (Fig. 4). Table 2. Cont.   This could very well be a consequence of an ongoing inflammation or possibly also indicate an underlying metabolic difference. Glutamine is converted to glutamate by the enzyme glutaminase (GLS). In turn, glutamate can be converted to proline and subsequently catabolized by the enzyme proline dehydrogenase (PRODH) resulting in the production of reactive oxygen species and apoptosis [23]. In the present study, we show that the expression of GLS is down-regulated and PDK1 is up-regulated in cases. Interestingly, a previous study has shown that cell lines with a known familiar mutation for amyotrophic lateral sclerosis (ALS) have the same expression pattern, with up-regulated PDK1 and down-regulated GLS, as compared to the wild-type cell line [24]. PRR5L (also called Protor-2) belong to the TOR signaling pathway. Our results show an up-regulated PRR5L expression in cases (Fig. 4). Like DUSP10, the protein product from PRR5L has been shown to stimulate an increased TNF-a expression [25]. Another gene, connected to the MAPK pathway and which was identified both by our two-locus interaction analysis and in significant biological functions implied by IPA, was the APPL1 gene. APPL1 is a binding partner of the protein kinase Akt2 and a key regulator of insulin signaling [26]. It takes part in adiponectin signaling to stimulate activity of p38 MAPK in muscle cells [27] and is a critical regulator of the crosstalk between adiponectin signaling and insulin signaling pathways [28]. We could detect expression of both APPL1 and APPL2 in small intestinal biopsies and a significantly lower expression of APPL2 was detected in the CD autoimmunity cases as compared to controls (Fig. 4). Lower expression of APPL2 levels lead to enhanced adiponectin stimulated glucose uptake and fatty acid oxidation [29]. A SNP (rs10861406) included in the top 603 list was located upstream of the APPL2 gene, however the promotor of this gene was on the opposite side of a recombination hotspot and therefore not included in the gene list for pathway analyses.
The most significant finding from our non-stratified linkage GWAS analysis was the association with the PPP1R12B gene region. PPP1R12B is involved in smooth muscle contractibility and mediates binding to myosin [30]. Myosin light chain phosphatase from smooth muscle consists of a catalytic subunit (PP1c) and two non-catalytic subunits, M130 and M20. The two non-catalytic subunits are both encoded by the PPP1R12B gene. The M130 transcript was not differentially expressed between CD autoimmunity and control patients while the small subunit ''M20'' showed a significantly higher expression in patients with CD autoimmunity. (PPP1R12B_22 in Fig. 4) Several other genes located close to top markers such as the PPP3CA, ACTN1, MYO1B, MYO5A, MAPK1, PRKCH, PRKCQ, PRKACB, PRR5L and NTS genes, are connected to smooth muscle when examining their function by using KEGG [31] and Gene Ontology [32]. The second most significant region in the HLA-stratified analysis after DUSP10 contains the SVIL gene. The product of this gene has been suggested to bind LPP [33]. In our two-locus interaction analysis, the LPP locus and a locus containing KIF13A was one of the 101 interaction pairs. KIF13A is a motor protein, which shuttles vesicles containing AP-1 and the mannnose-6phosphate receptor [34]. KIF13A was significantly down-regulated in intestinal biopsies from CD patients in our gene expression analysis (Fig. 4). SVIL is associated with cell-focal adhesions (substrate contacts), which are important for rapidly moving cells such as for example immune cells but also for motility and polarity of intestinal epithelial cells. SVIL mRNA was down-regulated in our gene expression analysis, however, not significant after correction for multiple testing.

Proline and Glutamine Metabolism -Part of a ''Danger Signal''
Amoebiasis was one of the nominally significant pathways in the GeneTrail analysis of genes surrounding the two-locus interaction SNPs ( Table 7). Several of these genes were also present together with DUSP10 and the MHC class II genes in the two most significant IPA generated networks (marked in bold text in    Table 8). Another gene present in these networks was the gene encoding for the immune molecule CD40 (associated SNP rs6065961, Table S1). CD40 has been shown to regulate immune responses to another parasite, Leishmania Major, by shared signaling through p38 MAPK and ERK1/2 [35]. CD40 also regulates DUSP expression and activity, which in turn contribute to anti-leishmanial functions [35]. It has been suggested that Leishmania Major inhibits CD40-triggered p38 MAPK signaling as part of an immune evasion strategy [36]. Another overrepresented category from GeneTrail was the extracellular matrix (ECM) ( Table 6). Also, in the two most significant Ingenuity networks from the 603 marker analyses, ECM molecules and matrix metalloproteinases (MMPs) were included ( Table 8). The ECM represents a major barrier to parasites like amoebas and leishmania. Parasites produce a wide variety of proteases to break down the ECM in order to access essential nutrients and invade host tissue [37]. A different situation when the ECM is degraded is during nutrient deprivation. In this way the ECM can provide energy for starving host cells. Just like gluten, the ECM has an unusually high proline content. MMPs are enzymes, which break down ECM making proline readily available as a nutritional source. Pandhare and co-workers have shown that energy or nutrient stress activates MMPs as well as the degradation of proline and furthermore demonstrated that, as the levels of glucose decreased to 1 mM and lower in the medium, intracellular proline increased almost 2-fold [38]. If gluten lingering in the intestine conveys a signal of ECM degradation (due to increased proline levels), several other mechanisms will most likely signal that there is food available at the same time (salivary secretion as one example is shown in Table 6). In this case, the immune system will rule out starvation as a possibility and the only other sensible option would be to search for an invasive intruder breaking down the ECM. The autoantigen in CD, TGM2, counteracts proteolysis and degradation of ECM by crosslinking ECM proteins [39]. If DUSP10 and PRRL5 upregulate TNF-a and subsequently TGM2 [17,18,25], in CD, the purpose may very well be for TGM2 to help prevent an apparent or illusory pathogenic invasion. It has also been shown that downregulation of SVIL protects against ECM invasion by pathogens [40]. In our gene expression analysis SVIL was nominally significantly down-regulated in cases (Table 9).
When the body ''senses'' a pathogen disturbing energy balance or breaking down ECM, but there are no pathogenic antigens present, maybe there could be a risk that ''self'' antigens become our immune systems futile attempt to rid the perceived pathogen. In HLA-DQA1*02/05 and HLA-DQB1*02 carriers, peptides derived from TGM2 could constitute such ''self'' antigens. It is possible that individuals carrying other HLA molecules still respond to this ''phantom pathogen'' and that under these circumstances, various other antigens present in the intestine at the time could become triggers of other autoimmune diseases. If the expression or presence of an autoantigen, like TGM2, was stimulated by the disturbed proline/glutamine homeostasis, it can explain why symptoms in CD also disappear by withdrawal of gluten.

Conclusion
At least four major functional components together with gluten, all seem to play a role in forming an individual's risk for CD: 1) polarity and epithelial cell functionality, e.g. nutrient/vesicle transport, proliferation and apoptosis, important for cell migration from the crypt to the shedding (apoptosis) at the apical villi. 2) intestinal smooth muscle, which is important for the movement of the bowel as well as the villi. 3) growth and energy homeostasis, which includes proline and glutamine metabolism, and finally 4) the innate and adaptive immune system.  A slight dysfunction combining these categories together with gluten consumption would result in a metabolic imbalance which in turn could convey enough stress or ''danger signal'', to trigger the immunological process and tissue destruction. A schematic illustration showing a rough outline of a possible disease model is presented in Figure 6.
In this study, we identified DUSP10 to be significantly associated with celiac disease. We also identified mechanisms, which we believe influence the risk of developing disease. Our data points towards genes that are involved in cancer as well as metabolic and cardiovascular diseases. Besides understanding how they work in celiac disease, our findings could also have consequences for these other common diseases.
Whole genome analysis allows for discovering completely unknown mechanisms behind disease. Even if the discovered genes and gene variants won't be able to predict who will develop disease in the future, they can be used to identify the underlying molecular pathways that influence disease. These molecular pathways would then be valuable targets for drug intervention. Our data provides new insights and hypotheses to the research field of CD and autoimmunity. However, the functional variants behind associations as well as mechanisms causing differences in gene expression and if and how these are relevant for disease, remains to be identified.

Ethics Statement
The regional ethics board in Gothenburg approved this study and participants in the study gave written informed consent after being fully informed about the aim of the study. For all children in the study, parental written consent was obtained.

Study Population
A total of 106 families with multiple affected individuals, mostly nuclear families with an affected sib pair (ASP) were collected from Sweden and Norway. There were 403 subjects and 97 families with DNA to complete the analysis. A total of 226 of the family members had CD, including 20 parents. The makeup and selection process regarding the families has been described previously in detail [41].
Small-intestinal biopsies, for the gene expression analysis, were collected at four pediatric clinics in Sweden: Skåne University Hospital in Malmö, Sach's Childrens' Hospital and Karolinska   [42,43]. Patients being both tTG IgA and IgG positive (.4 U/ml) were included as cases (we use the term ''CD autoimmunity'' for these cases [44]) and all other individuals were included as disease controls. The LPP gene was run in a first set of 42 cases (30 females, 12 males) and 38 control patients (21 females, 17 males). Cases and controls had a mean age of 7.3 and 11.1 years respectively. The remaining genes were run in a second set of 46 cases (25 females, 21 males) and 52 control patients (27 females, 25 males) with a mean age of 7.4 years for cases and 11.8 years for control patients at the time the biopsy was taken.

Gene Expression Analysis
We performed quantitative gene expression analysis using duodenal biopsies from CD autoimmunity patients and control patients. Biopsies were immediately put in RNAlater solution (Life Technologies, CA, USA). Total RNA was extracted using the miRNeasy Mini Kit (QIAGEN, Germany). RNA was converted to cDNA and quantitative PCR was run using TaqMan chemistry and the ABI7900 SDS instrument (Life Technologies, CA, USA). Seven control genes were evaluated using GeNorm [45] (ACTB (Hs00357333_g1), B2M (Hs99999907_m1), EPCAM (Hs00158980_m1), GUSB (Hs99999908_m1), HPRT1 (Hs99999909_m1), MUC1 (Hs00159357_m1), PGK1 (Hs99999906_m1)) and the geometrical mean of ACTB, EPCAM, and PGK1 were selected as reference for the relative quantification analysis (Delta-Delta Ct method). A total of 34 expressed genes located close to some of the most significantly associated SNPs were evaluated ( Table 9). The top associated genes from our Linkage GWAS (DUSP10, SVIL, PPP1R12B) were selected as well as several genes from the two-locus interaction analysis and pathway analyses including LPP which is the top associated from the GWAS by Dubois et al. Also the RGS genes (RGS1, 2 and 5) and GLS show genome-wide association in the study by Dubois et al and are also present in our two-locus and pathway analyses.

Genotyping and Imputation
Samples were genotyped using two different SNP arrays, 211 samples with Human Omni Express and 192 samples with Human 660W-Quad (Illumina Inc, CA, USA). A total of 308,246 markers were available on both arrays and were therefore genotyped in the entire material. For the remaining 682,470 and for sporadic missing values we performed genotype imputation using the Impute 2 software [46], with the Hapmap 2 (rel. 24 Build 36) as a reference.
All individuals in the same family were located on the same plate. Quality control was first performed separately for the two arrays. SNP markers with less than 97% call rate in either of the two arrays were excluded.
Mendelian errors were detected by PLINK [47], 125,874 family-wise mendelian inconsistencies were set to missing (in each family the genotypes were set to missing for all subjects if there were any mendelian inconsistencies for a specific SNP).

Statistical Analysis
Linkage. For the linkage analysis only markers from both platforms were considered. From this set of 271,078 common SNP markers a LD pruned set of 105 539 SNPs were selected using PLINK. Parameters were a window size of 50 and R 2 ,0.5. The Decode genetic map as supplied by Illumina was used to run nonparametric linkage using Merlin version 1.1.2 [48] with the NPL Expression (e.g. mRNA levels) of these genes was either up-or down-regulated in small intestinal biopsies from CD cases compared with control patients. Effect direction is presented for cases with control group as a reference. The selection column indicates if the gene was selected due to its presence in two-locus or pathway analyses. ''Top'' indicates top SNP in the present GWAS and ''top previous'' indicates that it was present in the GWAS by Dubois et al.
a Gene not included in the gene list used for pathway analyses due to recombination between associated SNP and gene promotor: rs10861406.
(APPL2), rs882820 (ADCY9). However, possible regulatory site could be close to associated SNP and influence gene expression. all method [49]. Marker allele frequencies were estimated from the founders.
The imputation analysis provides us with (posterior) probabilities for each of the possible genotypes at each locus and to utilize all posterior probabilities, we performed an analysis where we use the expected values of the transmission counts. The test statistic will then have the following form, T imp has approximately the same distribution as the test statistic T in [50].
Stratified TDT. We implement a stratified TDT analysis where trios are split into a low-risk and a high-risk group based on the HLA genotype of the affected offspring. Children carrying the HLA-DQA1*02/05 risk allele and homozygous for the HLA-DQB1*02 risk allele (i.e. individuals carrying the DR3/DR3 or the DR3/DR7 haplotypes) were put in the ''high-risk'' group and the remaining children were put in the ''low-risk'' group. The rationale behind this is explained in the introduction and further information about this stratification can be found in our previous Linkage study [14]. A standard TDT analysis, with the 0.95 cut-off for imputation probabilities, was applied to each of these groups using PLINK [47].
Test for two-locus interaction. To examine possible interactions between marker variants, we used a pairwise test based on the one introduced by Kotti [51]. Consider two biallelic markers without linkage disequilibrium between their alleles. In the general model M G the penetrance matrix has 9 parameters,

5,
Let n be the 363 matrix of genotype counts among the cases for the two markers, and let m be the corresponding matrix for the non-transmitted allele combinations. The likelihood for the models is For this analysis we use one affected subject from each family and markers were chosen based on the expected counts TDT (equation 1) and three different inclusion criteria: 1. P-value less than 3.0610 24 . 2. P-value less than 0.01 in our analysis and with a p-value less than 0.05 in the GWAS by Dubois et al. [3] and if the product of these p-values were less than 5.0610 25 and the association were in the same allelic direction. 3. An allele transmission ratio of ,0.2 or .5 combined with a pvalue less than 2610 23 .
We defined 383 regions using the inclusion criteria above (Fig. 2 and Table S1) a region consisted of a set of markers where the distance between adjacent markers was less than 100 kb. With these regions defined we analyzed all pairwise interactions using a Likelihood Ratio (LR) tests comparing the following four models: N M 0 : None of the two loci is associated with CD, N M R : Heterogeneity model [52], with penetrance where a i and b j are the penetrance factors for the genotypes A i and B j [53] respectively. The restricted model used in [51] is the multiplicative model. We use the M G -versus-M M test to filter out false positives, based on that if one or both of the SNPs were marginally significant by chance, then the joint distribution (penetrance) of these markers should follow a multiplicative model.
We have the likelihood ratio statistic T jk~{ 2 log max L j max L j T jk will follow a x 2 distribution under the restricted model if M j is nested in M k . The maximum likelihood estimates of the penetrance parameters and allele frequencies do not have a simple explicit expression, so to maximize the likelihoods we use the function optim in the statistical software R.

Gene Selection
Out of the 603 SNPs selected from the three inclusion criteria ( Fig. 2 and Table S1), we were able to identify genes surrounding Figure 6. Proposed disease model. Illustrating a possible scenario for disease development. Genetic variation contributing susceptibility to disease can be found in at least four, somewhat overlapping, biological functions. The result is an ''overload'' or imbalance of proline vs glutamine. Due to the abundance of proline within the extracellular matrix (ECM) as well as in gluten, the proline from gluten is interpreted as degradation of ECM. When the body is not starving, the ECM is normally not degraded, unless there is a pathogen attempting to break through this barrier. The immune system mounts an attack against an invasive ''phantom pathogen'' which is believed to degrade the ECM. When proline is catabolized, reactive oxygen species (ROS) are released. In order to start re-building and crosslinking ECM molecules, Tissue transglutaminase (TGM2) expression is up-regulated by TNFa which in turn is stimulated by DUSP10 and Protor-2 (PRR5L). This rebuilding of the ECM counteracts the degradation by the imagined pathogen. However, the phantom pathogen remains and the adaptive immunity is brought in. Searching for antigens, it finds an abundance of TGM2 beside the ECM and forms antibodies against its own soldier. Some susceptibility genes can be found in the center of this model and some can be found within the spiral. Genes like HLADQ and other genes from the adaptive immunity are likely to be found in the spiral. doi:10.1371/journal.pone.0070174.g006 444 SNPs using GRAIL [54]. Grail uses known recombination hotspots in order to limit the region of interest surrounding each SNP marker. Genes around the remaining SNPs were identified with the Genome Browser (http://genome.ucsc.edu) and the 5 closest genes within 250 kb from the associated SNPs were included. In cases where there were no genes within this distance we included the closest gene.

Pathway Analysis
We analyzed connections between genes in different regions, using GeneTrail [13] and the Ingenuity Pathway Analysis (IPA) software (Ingenuity Inc. CA, USA). Within each associated region, all but one gene from the same gene family were removed. This was done in order not to amplify the significance of homologous gene clusters, (i.e. chemokine receptor-, interferon-and histonegene clusters).

Supporting Information
Table S1 A selection of 603 top associated SNPs and the two top HLA SNPs. Based on three inclusion criteria, 603 SNP markers and 383 regions were identified. Our TDT results, where T and U are the expected transmission counts (based on all the posterior imputation probabilities) and corresponding results from Dubois et al. [3]. (DOCX)