Genome-wide association study identifies five novel susceptibility loci for Crohn's disease and implicates a role for autophagy in disease pathogenesis.

Two independent replication studies confirm novel IBD loci

While certainly encouraging, even an overall high-quality data set and robust analyses cannot be assumed to perfectly filter out all artifactual sources of false positives. We therefore tested the most significantly associated SNP in the 23 independently associated regions with 5 p<5×10⁻⁵ in an independent set of 650 individuals with ileal CD and their parents (Replication Cohort #1), to attempt replication of these findings in a robust family-based setting and with an alternative genotyping technology (iPlex; Sequenom). Although we did not succeed in typing three of these SNPs (and the replication of CARD15 and IL23R had already been established and are not considered here) the remaining 18 SNPs had high quality genotype data that passed our quality thresholds. After removal of DNAs with low genotype call rates and families with an excess of Mendel errors, the final analyses included 530 trios. Association testing of this final dataset from Replication Cohort #1 identified that six out of the 18 SNPs had significant evidence of replication (Table 1; see Supplementary Table 1 online): rs2241880 (p=0.00089); rs16853571 (p=0.00037); rs224136 (p=0.00021); rs4821544 (p=0.019); rs8050910 (p=0.024); rs11617463 (0.037). To provide further evidence for the association of these replicated SNPs, we typed these six in another independent cohort of 353 patients with ileal CD and 207 controls (Replication Cohort #2) on yet another independent genotyping platform (Taqman; Applied Biosystems). One of these failed to type but the other five were successfully typed and provided additional evidence of replication (Table 1; see Supplementary Table 1 online). Considering the two replication studies together (using Mantel-Haenszel statistics), two of these SNPs were unequivocally confirmed as novel CD risk factors (rs2241880, p<5×10⁻⁸; rs224136, p<5×10⁻⁷). In combination with the genome-wide scan data, these two SNPs are considerably in excess of the most conservative genome-wide thresholds (rs2241880, p<10⁻¹³; rs224136, p<10⁻¹⁰). Three additional loci (rs16853571, rs4821544, and rs8050910) showed very strong evidence of replication (p<0.01). Even when setting aside the two most significantly associated novel SNPs (rs2241880 and rs224136), observing replication of 3/16 SNPs with p<0.01 is highly unlikely to occur by chance (p=0.0005). Since there is only a 1% chance of observing even two such results by chance, it is highly likely that these three SNPs represent true CD susceptibility variants, but further confirmation is required to achieve the same level of certainty as attained for the IL23R, rs2241880 and rs224136 SNPs identified in this study.

Association mapping implicates novel IBD genes

Given the significant evidence of replication for these five novel loci, we examined the linkage disequilibrium structure surrounding the replicated SNPs using the genotype data from the GWA screen, and examined the association mapping data from the GWA, replication studies as well as some additional higher density typing in order to map these association signals to specific genomic regions and genes. As shown in Supplementary Figure 1 online, in four of these instances, the association mapping implicates specific genes, whereas in the fifth it implicates an intergenic region.

Specifically, the LD structure and association mapping around the most associated SNP (GWA, p=6.4E-08; Replication, p=4.1E-08), rs2241880, implicates a region on chromosome 6 2q37.1 containing a single gene known as ATG16 autophagy related 16-like 1 (S. cerevisiae) or ATG16L1 gene. In fact, this SNP encodes a non-synonymous amino acid change – an alanine to threonine substitution in exon 8 – Ala197Thr. Logistic regression analyses conditional on the Ala197Thr in the family-based samples indicates that this coding variant can fully explain the association signal to this locus – therefore we consider this to be the causal risk variant. In all of the datasets tested in this study, the threonine allele is the minor allele and has a protective effect; in the GWA dataset the estimated odds ratio is 0.692 (95% CI of 0.608-0.788). Interestingly, the threonine is the ancestral allele and is the major allele in the YRI, JPT, and HCB samples typed in the International HapMap Project (see Supplementary Fig. 2 online). ATG16L1 is part of a family of genes involved in autophagy – a biological process involved in protein degradation, antigen processing, regulation of cell signaling and many other pathways essential to the initiation and regulation of the inflammatory response – and therefore this association suggests that this biological process is likely to have an important role in disease pathogenesis.

The next SNP with incontrovertible association to ileal CD is rs224136 (GWA, p=7.9E-06; Replication, p=2.9E-07). The association mapping in this case, however, implicates a genomic region on chromosome 10q21.1 of approximately 70 kb that does not contain any known protein coding genes. This associated intergenic region is flanked on the centromeric side by the zinc finger gene known as ZNF365 and on the telomeric side by a predicted protein (c10orf22) and by the known early growth response gene EGR2, also a zinc finger protein. Although there is little known about the ZNF365 and c10orf22 genes, mutations in EGR2 are known to lead to the development of congenital hypomyelinating neuropathy (MIM#605253) or of Charcot-Marie-Tooth disease (MIM#607678). Potentially more relevant, however, is the observation that EGR2 appears to be a negative regulator of T cell activation in mice¹², and although speculative, genetic variation within this intergenic region could play a role in the control of expression of EGR2 and influence biological function.

The replicated SNP rs16853571 (GWA, p=7.7E-07; Replication, p=8.4E−03) is located in the promoter region of the Paired-Like Homeobox 2B (PHOX2B) gene on chromosome 4p13 (see Supplementary Fig. 1 online) - 2143 bp upstream of the first known exon and 10 bp upstream of a putative alternate first exon. PHOX2B is a homeodomain transcription factor, predominantly expressed in differentiating neurons. Human mutations in the PHOX2B gene have been implicated in disorders of neural crest development, congenital central hypoventilation syndrome and Hirschsprung disease^13,14. Inactivation of PHOX2B in mice disrupts noradrenergic differentiation throughout the nervous system. Autonomic ganglia fail to develop properly and these mice have striking defects in enteric neurons^14-16.

The genomic region identified by rs4821544 (GWA, p=2.9E-05; Replication, p=9.0E-03) also contains a single gene known as NCF4 (rs4821544 itself is in intron 1). NCF4 encodes the 7 p40phox protein that has previously been demonstrated to play an important role in NADPH oxidase activity and the generation of reactive oxygen species (ROS) production upon phagocytosis^17,18; both important for mounting an effective anti-microbial response.

Finally, rs8050910 (GWA, p=3.3E-05; Replication, p=8.5E-03) is located within the 4^th intron of the predicted gene FAM92B located on chromosome 16q24.1. This gene prediction is supported by multiple mRNAs and spliced ESTs, as well as significant conservation across multiple species (data not shown). This gene encodes a protein with no known function or recognizable motifs.

It should be noted that when using a stratified analysis approach, none of these three loci appeared to have any statistically significant epistatic interactions with the CARD15, IBD5 or IL23R genes (data not shown). Though some association with UC is seen for rs224136, little or no evidence is seen at the other newly identified SNPs (see Supplementary Table 1 online).

RNA and protein expression analyses provide initial biological context to novel genes

Since we were interested in determining more about ATG16L1's function and biological context, we examined the expression distribution of ATG16L1 and other known autophagy components by quantitative real-time RT-PCR in a variety of epithelial and immune cell lines (see Fig. 2A, B and C). ATG5, 7 and 16 were broadly expressed, with SW480 cells having the lowest overall mRNA abundance relative to GAPDH controls. Apart from the markedly elevated levels of ATG7 seen in THP-1 cells, expression of each of the three components followed a similar pattern in the cells tested, as might be suggested by the stoichiometric relationship between these components in the homeostatic autophagy process. We also investigated the expression of ATG16L1 in primary human immune cells, to gain a better understanding of its role in the association with CD (see Fig. 2D). Amongst samples tested, we found that ATG16L1 was most highly expressed in the T cell compartment, with both CD4+ and CD8+ cells showing high levels (approximately 13- and 10-fold higher than placental control RNA, respectively). CD19+ cells (B cells) also showed almost 5-fold greater expression of ATG16L1 than both mononuclear cells and placental control RNA. Thus it is possible that ATG16L1 and its variants play roles in both the epithelial- and immune-driven aspects of CD.

Open in a separate window

Figure 2

Expression of autophagy components in human cell-lines and primary immune cells

Quantitative real-time PCR was used to determine the expression patterns of ATG5, 7 and 16L1 in a variety of human immune and epithelial cell lines. Panel A - Expression pattern of ATG5 in human cell lines. ATG5 seems broadly expressed in the cell lines studied, with between 1.5-and 4-fold variations in expression level across the cell lines. Panel B - Expression pattern of ATG7 in human cell lines. ATG7 expression levels appeared broadly comparable to those of ATG5, with the exception of THP1 cells, which exhibited a 35-fold higher level than SW480 cells. Panel C - Expression pattern of ATG16L1 in human cell lines. ATG16L1 expression was also broadly similar to those of ATG5 and 7 (with the exception of ATG7 expression in THP1 cells) with moderate variations in expression levels across the cell lines tested (1.5- to 3.5-fold). Panel D - Expression pattern of ATG16L1 in human primary immune cells. In a resting human immune cell RNA panel (Clontech, CA) ATG16L1 showed a clear peak of expression in the T cell compartment, with both CD4+ and CD8+ cells expressing high levels of transcript (15- to 20-fold more than control samples). Modest levels were detected in CD19+ cells, and low levels in CD14+ monocytes and unfractionated mononuclear cells. In all cases, real-time quantitative RTPCR reactions were performed in duplicate and the means plotted, error bars represent 1 standard deviation. Reactions containing no template RNA were performed to control for reaction contamination (water controls). Expression levels were normalized by comparison to GAPDH controls and arbitrary relative expression units plotted where SW480 or placental RNA is equal to 1 (for cell lines, or primary cells respectively). All RNAs were isolated from resting cells in the absence of stimulation or activation.

PHOX2B is known to be expressed in the autonomic nervous system and enteric neurons¹⁶. Attempts to determine expression using real-time RT-PCR were unsuccessful, most likely due to expression being confined to specific cell types and therefore overall transcript levels being low in the tissues studied. To overcome this limitation we obtained a specific antibody and stained mouse and human gut samples (ileum and colon respectively). These clearly showed that the PHOX2B protein was expressed in a specific subset of epithelial cells, possibly neuroendocrine cells (See Supplementary Fig. 3 online).

NCF4 encodes p40phox, a component of the NADPH oxidase complex required for optimal ROS generation in immune cells. Expression is restricted to haematopoietic cells including neutrophils, monocytes, eosinophils, mast cells and basophils (See Supplementary Figure 4 and reference ¹⁹). The p40phox protein serves to enhance delivery of p47 and p67phox to the membrane and thus promotes high-level ROS generation. However, it's role in cells that do not generate high-levels of ROS and lack expression of p47 and p67phox remains unknown.

Genotyping methods

For the genome-wide association study, approximately 750ng of genomic DNA was used to genotype each sample on the Illumina HumanHap300 Genotyping BeadChip (Illumina, San Diego) at the Feinstein Institute for Medical Research. Samples were processed according to the Illumina Infinium 2 assay manual. Briefly, each sample was whole-genome amplified, fragmented, precipitated and resuspended in appropriate hybridization buffer. Denatured samples were hybridized on prepared HumanHap300 beadchips for a minimum of 16 hours at 48°C. Following hybridization, the beadchips were processed for the single base extension reaction, staining and imaging on an Illumina Bead Array Reader. Normalized bead intensity data obtained for each sample was loaded into the Illumina Beadstudio 2.0 software which converted fluorescent intensities into SNP genotypes. The replication study of the all of the putative loci (Table 1), and the additional genotyping of the ATG16L1 region was performed using primer extension chemistry and mass spectrometric analysis (iPlex assay, Sequenom, San Diego) using Sequenom Genetics Services (Sequenom, San Diego). TaqMan MGB technology was used to perform the genotyping of the six variants with nominal replication in replication study #1 using the design available from Applied Biosystems and following the manufacturer's recommendations (Bulletin #4317594) in the Cedars-Sinai Inflammatory Bowel Disease Center cohort. Quality of the SNP data was checked by reproducibility with 5% of samples duplicated.

Analysis of GWA data

Overall, 98.6% of genotypes were called in this experiment; however, before analysis, quality filtering of both samples and SNPs were performed to insure robust association tests. Based on an evaluation of empirical distributions, data quality, and likely introduction of false positive associations, we required that samples pass a 93% genotyping call rate threshold and SNPs pass a 95% call rate threshold in order to be included in the analysis (see Supplementary Fig. 6 online). The data from the genome-wide association studies were used to detect possible relatedness in the case-control cohorts. The HWE test was performed on the genotype data from controls and we investigated the relationship between HWE and genotype yield (call rate). The call rate distribution suggested a 95% and 93% genotype completion threshold for inclusion of samples and SNPs, respectively in the genetic association analyses. At sample call rates of 90% and below, there was an elevation in observed heterozygosity suggesting either a bias in missing data or the presence of false heterozygous calls. To avoid this, we selected a threshold greater than 90 at a point consistent with the tail of successful samples. Below SNP success of 95% there was an excess of markers out of HWE (P-value <0.001 in controls) and a significant bias in missing data between cases and controls. We looked at segments of the data just 15 above and below this cutoff and observed that SNPs with a 90-95% call rate had much more substantial inflation than those in the 95-97% call rate range (genomic control correction 1.34 vs. 1.16) indicative of significant excess of false positives due to lower data quality/biases in missing data so we opted not to lower this threshold further. As a final step, using identity-bystate (IBS) counts from this data set, we identified 8 duplicate samples and 10 additional pairs of samples that shared 10% or more of their genome by descent (first cousins through full siblings) - we eliminated one member of each of these 18 pairs of samples. The data from the Jewish and non-Jewish case-control cohorts were analyzed jointly using a Cochran-Mantel-Haenszel chi-squared test; we used the test as implemented in R (http://www.r-project.org/) with the option that gives an exact p-value.

Analysis of replication study data

For the replication study, the most significantly associated SNP from the top 23 ranking independently associated regions (p-value=5×10⁻⁵ in the combined screening analyses) were selected. Of these, 2 failed in assay design, one failed to type, and 2 were located within previously reported associated loci, leaving 18 SNPs in the replication analysis. In addition to these, 19 additional SNPs were selected to test the variability in the ATG16L1 region. This set of SNPs was evaluated for replication in a family-based ileal-CD cohort composed of 650 mother-father-affected offspring trios. The genotype data from the family-based cohorts were used in determining Mendelian inconsistencies and departures from HWE. After removal of failing SNPs (monomorphic, call rate below 75%, and causing a high number of Mendel errors), families with excess Mendelian inconsistencies and uninformative families (where a parent or the sole affected offspring had a genotyping call rate below 90%), there remained 530 affected trios (from 431 independent nuclear families) with both parents and at least one affected offspring genotyped. Following the removal of additional SNPs causing excess Mendel errors, the final genotyping call rate for each of the 18 replication markers, as well as the 19 ATG16L1 tagging SNPs, ranged from 90-100%, and Mendelian inconsistencies were below 4 for all markers and less than 4 per family. Single marker association tests were performed in the family data using PLINK (http://pngu.mgh.harvard.edu/~purcell/plink/). We then tested the replicated SNPs in the set of 530 independent IBD trios (post-QC) from the NIDDK IBDGC and in the independent Cedars-Sinai cohort of 350 patients with ileal CD and 207 controls, and onetailed p-values are reported. The analysis of ATG16L1 SNPs conditional on Ala197Thr was performed by logistic regression in the family-based sample using the software WHAP (http://pngu.mgh.harvard.edu/~purcell/whap).

mRNA analysis of ATG16L1

RNA templates were purified from cell lines using the RNeasy miniprep kit (Qiagen). In the case of human primary immune cell subsets RNA was purchased from a commercial supplier (Human Blood Fractions MTC Panel, Clontech). In all cases cDNA was made from RNA templates using a standard reverse transcriptase protocol (iScript, BioRad), with 500 ng of RNA per reaction. The resulting template cDNAs were diluted 1:10 in nuclease-free water and 10 ul of this template used for each 40 ul qPCR reaction. Quantitative RT-PCR was performed using real-time SYBR Green method on a BioRad iQ5 thermal cycler, using iQ SYBR Green Supermix (BioRad) and specific primers. Primer sequences used were selected using MGH Primerbank (http://pga.mgh.harvard.edu/primerbank/); all primer sequences are shown in Supplementary Table 2. PCR products were visualized on a 2% agarose gel to confirm correct band sizes. Each reaction was performed in duplicate, with final calculations resulting from means of duplicate wells. Normalization for cDNA quantity was performed with GAPDH control primers for each template and final abundance figures adjusted to yield an arbitrary value of 1 for control templates (SW480 for cell lines or placental control RNA for primary cells) using the delta-delta Ct method.

Immunohistochemistry for localization of PHOX2B in mouse and human gut tissues

Paraffin-embedded tissues from both mouse ileum and human colon biospies were dewaxed, rehydrated, blocked in serum and stained using rabbit anti-PHOX2B antibody (Sigma), incubated overnight at a dilution of 1:500. Tissues were subsequently incubated with biotinylated anti-rabbit secondary antibody, Vector ABC reagent and finally developed with DAB substrate (Vector Laboratories) as has been previously described⁴². Sections were mounted and observed under phase contrast to visualize tissue architecture.

Cell lines and bacterial strain

HeLa cells were obtained from ATCC and cultured in DMEM supplemented with 10% iron-supplemented calf serum (CSFe) (Hyclone) and 20 μg ml⁻¹ gentamicin sulfate (Sigma). Salmonella enterica sbsp. enterica serovar Typhimurium (S. Typhimurium), strain SL1344 was transformed with a DsRed2 expression plasmid (Clontech, CA) (DsRed2 is under the control of the lac promoter, which is constitutive in the absence of lacl) to facilitate bacterial visualization. The plasmid was maintained with 100 μg ml⁻¹ ampicillin selection during culture in LB broth.

siRNA directed against ATG16L1

We obtained modified siRNA duplexes (Stealth siRNA) directed against ATG16L1 (siRNA 1 & 2) as well as a non-targeting scrambled sequence (siRNA control) from a commercial source (Invitrogen). Cells were co-transfected in the absence of antibiotics with siRNAs and plasmids using Lipofectamine (Invitrogen) according to the manufacturer's instructions. We optimized 17 siRNA conditions using a flag-tagged overexpression ATG16L1 construct in HEK293 cells, allowing us to easily assess knockdown by Western blotting. HEK293 cells grown in 12 well plates were transfected with 500 ng Flag-ATG16L1 expression plasmid and 20 pmol siRNA duplex. For knockdown of endogenous ATG16L1, HeLa cells were plated onto 18mm glass coverslips in 12-well plates at a density of 1×10⁵ well⁻¹ and allowed to grow for 24 hours prior to transfection. Each well received 500ng of GFP-LC3 plasmid and 20 pmol siRNA duplex. Four hours later the medium was changed again and the cells were allowed to grow for a further 48 hours. Knockdown was confirmed for each experiment using real-time quantitative RT PCR with specific primers to ATG16L1 and normalized to GAPDH controls, with RNA isolated from specimen wells.

Induction of autophagy and bacterial infection

Autophagy was induced in HeLa cells using serum starvation or rapamycin. Cells were transfected as described above and after 48 hours the medium was changed to either 1% serum, or 10% serum plus 200 nM rapamycin for a further 24 hours. Ammonium chloride treatment was performed with a final concentration of 50 mM for 2 hours in 1% serum medium. S. Typhimurium infections were performed as previously described, with slight modifications⁴³. Briefly, S. Typhimurium SL1344 carrying a DsRed2 expression plasmid was grown overnight in LB broth containing 100 μg ml⁻¹ ampicillin at 37°C with aeration and subcultured at a dilution of 1:33 for a further 3 hours in LB. This culture was further diluted in DMEM 10% CSFe without antibiotics, to yield an m.o.i. of 100. Infections were allowed to proceed for 20 minutes, cells were washed once in complete medium containing 100 μg ml⁻¹ gentamycin sulfate, and then incubated in fresh high-gentamycin medium for 1 hour.

We thank the patients and their families for participating in these studies. We are grateful to all of the clinicians, research nurses and study coordinators for their essential contributions to the work. The authors would like to thank D. Caplan, G. Charron, C. Labbé, C. Lefebvre, and D. Miclaus for their help in the preparation of the manuscript, J. Adams for help with immunohistochemistry, A. Landry expression studies and to D. Altshuler for his critical reading of the manuscript. The National Institute of Diabetes and Digestive and Kidney Diseases IBD Genetics Consortium is funded by the following grants: DK62431 (S.R.B.), DK62420 (R.H.D.), DK62432 (J.D.R.), DK62423 (M.S.S.), DK62413 (K.D.T.), and DK62422 and DK62429 (J.H.C.). A.H.S. is on the Scientific Advisory Boards of Shire Pharmaceuticals, Schering (Canada), and Procter & Gamble Pharmaceuticals. The work on the Cedars-Sinai cohort was supported by project 1 of DK 46763 (J.I.R.). R.J.X. is supported by the following grants AI062773 and DK43351.