Common variants at the MHC locus and at chromosome 16q24.1 predispose to Barrett's esophagus.

Barrett's esophagus is an increasingly common disease that is strongly associated with reflux of stomach acid and usually a hiatus hernia, and it strongly predisposes to esophageal adenocarcinoma (EAC), a tumor with a very poor prognosis. We report the first genome-wide association study on Barrett's esophagus, comprising 1,852 UK cases and 5,172 UK controls in the discovery stage and 5,986 cases and 12,825 controls in the replication stage. Variants at two loci were associated with disease risk: chromosome 6p21, rs9257809 (Pcombined=4.09×10(-9); odds ratio (OR)=1.21, 95% confidence interval (CI)=1.13-1.28), within the major histocompatibility complex locus, and chromosome 16q24, rs9936833 (Pcombined=2.74×10(-10); OR=1.14, 95% CI=1.10-1.19), for which the closest protein-coding gene is FOXF1, which is implicated in esophageal development and structure. We found evidence that many common variants of small effect contribute to genetic susceptibility to Barrett's esophagus and that SNP alleles predisposing to obesity also increase risk for Barrett's esophagus.

Barrett’s Esophagus is one of the most common pre-malignant lesions in the Western world. It affects over 2% of the adult population and, unlike bowel polyps, lacks any proven effective therapy[1]. In the majority of cases, Barrett’s Esophagus is associated with chronic gastro-esophageal reflux disease (GERD), including esophagitis [2,3]. In addition there are structural changes, mainly hiatus hernia, in the lower esophagus in over 80% of patients[4]. This allows both acid and bile to remain immediately adjacent to the esophageal epithelium. The measured annual risk of esophageal adenocarcinoma (EAC) in Barrett’s Esophagus patient’s varies widely but is approximately 0.4-1% [5-7]. Notably the incidence of EAC has been rising by 3% each year for the last 30 years; it is now the fifth commonest cancer in the UK [8]. Despite modern multimodality therapy, the prognosis of EAC remains poor, with a 9-15% 5-year survival [9,10].

The etiology of Barrett’s Esophagus is not well characterised. Environmental factors, such as diet, are weakly associated with GERD, Barrett’s Esophagus and EAC, and obesity is a known risk factor for all three conditions[11]. There is also evidence implicating genetic factors: the relative risks are increased 2-4 fold for GERD, Barrett’s Esophagus and EAC when one first-degree relative is affected[12-17]. A segregation analysis of 881 pedigrees of familial Barrett’s Esophagus supports an incompletely dominant inheritance model with a polygenic component[18]. Extensive candidate gene and linkage searches have, to date, failed to identify genetic variants that are associated with risk of Barrett’s Esophagus[19].

As part of the Wellcome Trust Case Control Consortium 2 (WTCCC2) study of 15 common disorders and traits, we present the results of the first genome-wide association study of Barrett’s Esophagus susceptibility. Using a discovery cohort from the UK (with case samples from Aspirin and Esomeprazole Chemoprevention Trial of Cancer in Barrett’s Oesophagus (AspECT)[20]), and five replication cohorts (including case samples from CHemoprevention Of Premalignant Intestinal Neoplasia (ChOPIN) and Esophageal Adenocarcinoma GenEtics Consortium (EAGLE) studies[9,20]), we identified two variants associated with Barrett’s Esophagus, each with combined evidence at P<5×10−8. The analysis workflow is outlined in Supplementary Figure 1 and characteristics of the case and control samples used can be found in Supplementary Table 1 and Online Methods.

For the discovery analysis, cases with histologically confirmed Barrett’s Esophagus (see methods) were recruited from sites across the UK (Supplementary Table 2). Population controls were taken from the WTCCC2 common set of 1958 Birth Cohort (58C) and National Blood Service (UKBS) samples as previously described[21]. The case individuals were genotyped on the Illumina 660W-Quad array and controls were genotyped on the Illumina custom Human 1.2M-Duo array, with the analysis performed on the overlapping set of SNPs. Following quality control (see Online Methods, Supplementary Note, Supplementary Figure 2 and Supplementary Table 3), a total of 521,744 SNPs typed in 1,852 cases and 5,172 controls (2,499 UKBS and 2,673 58C) were included in the discovery analysis.

Association analysis was carried out under a logistic regression model as implemented in SNPTEST. The genomic over-dispersion factor[22] λ was 1.10 and this was reduced to 1.05 when incorporating the first principal component as a covariate, suggesting that population structure was not a major problem in the discovery analyses (Supplementary Figure 3). For all of the following results presented, unless otherwise stated, the first principal component was used as a covariate.

Following analysis of the genome-wide association results (Figure 1), we adopted a staged approach to replication, outlined below and in Supplementary Figure 1.

Figure 1

Plot of the genome-wide association results after fitting the multiplicative model in SNPTEST. Results shown for the 521,744 SNPs passing quality control filters. Chromosomes are coloured dark blue and light blue alternatively, as labelled on the x-axis. The y-axis shows the −log10 P values. Regions in red show the loci newly identified as associated with BE, as described in Table 1.

Stage 1

100 SNPs that showed evidence of association in the discovery data (at P<5×10−4) were analysed in another UK sample set. This comprised 1,105 cases from ChOPIN and EAGLE and 4,421 controls from the 58C control dataset, all genotyped on the Illumina Immunochip[23] (WTCCC2 contributed SNPs to the Immunochip design to allow for its replication studies, and the set of 100 SNPs followed up in our Stage 1 replication were all on the Immunochip), and a further set of 2,578 UK controls (the People of the British Isles (PoBI) collection [24]) genotyped on the Illumina custom Human 1.2M-Duo array. Results of this first stage of replication are shown in Supplementary Table 4.

Stage 2

The 16 top SNPs (Pcombined<10−5) from meta-analysis of the discovery and Stage 1 replication were replicated in silico in a Dutch collection of 473 cases and 1,780 controls genotyped on the Immunochip[23]. Results from Stage 2 replication are shown in Supplementary Table 5.

Stage 3

Two SNPs with Pcombined<5×10−8 after Stage 2 replication (rs9257809 on chromosome 6p21 and rs9936833 on chromosome 16q24) were studied in three additional replication sample sets. They were directly genotyped in an Irish cohort of 245 cases and 473 controls and a UK cohort of 1,765 cases and 1,586 controls, and data from these SNPs was retrieved from the BEACON consortium for 2,398 cases and 2,167 controls from European, Australian and American individuals with European ancestry.

After these three stages of replication, the two SNPs on chromosome 6p21 and 16q24 showed compelling evidence for association, with combined P values of 4.09×10−9 for rs9257809, OR(95%CI)=1.21(1.13-1.28) and 2.74×10−10 for rs9936833, OR(95%CI) =1.14(1.10-1.19) (Table 1, Figures 2, 3).

Table 1

Loci associated with risk of Barrett’s Esophagus

ChrPosition*	rsID	Riskallele	Discovery(1852/5172)		Stage 1		Stage 2		Stage 3						CombinedP valueOR(95%CI)
					UK replication 1(1105/6819)		Dutch Replication(473/1780)		UK Replication 2(1765/1586)		Irish Replication(245/473)		BEACON(2398/2167)
			RAFCase/Con	P valueOR(95%CI)	RAFCase/Con	P valueOR(95%CI)	RAFCase/Con	P valueOR(95%CI)	RAFCase/Con	P valueOR(95%CI)	RAFCase/Con	P valueOR(95%CI)	RAFCase/Con	P valueOR(95%CI)
6p2129464310	rs9257809	A	0.90/0.87	2.78×10−71.36(1.21-1.53)	0.89/0.87	0.06061.14(0.99-1.31)	0.91/0.87	9.16×10−41.45(1.16-1.80)	0.88/0.87	0.1511.11(0.96-1.29)	0.85/0.86	0.7670.96(0.71-1.29)	0.91/0.90	0.0831.13(0.98-1.31)	4.09×10−91.21(1.13-1.28)
16q2484960619	rs9936833	C	0.42/0.38	8.18×10−61.20(1.11-1.29)	0.42/0.37	6.24×10−51.21(1.10-1.32)	0.35/0.34	0.4021.07(0.92-1.25)	0.39/0.39	0.8101.01(0.92-1.12)	0.41/0.39	0.4681.09(1.06-1.11)	0.40/0.36	5.13 ×10−41.16(1.07-1.27)	2.74×10−101.14(1.10-1.19)

Discovery and replication results at the lead SNPs at the two loci for which there is combined evidence of P< 5×10−8. P values are two-sided. ‘RAF’- Risk allele frequency, *NCBI Build 36. The number of cases and controls, respectively, in each cohort is shown under the title of the cohort.

Figure 2

Regional association plot of the associated loci as detailed in Table 1, showing the signal at the lead SNP. The −log10 P values for the SNPs are shown on the upper part of each plot. SNPs are coloured based on their r2 with the labelled hit SNP which has the smallest P value in the region. r2 is calculated from the 58C data. The bottom section of each plot shows the fine scale recombination rates estimated from individuals in the HapMap population, and genes are marked by horizontal blue lines.

Figure 3

Forest plots showing evidence for association at each of the loci described in Table 1. The effect size and 95% CI are shown to the right of the cohort label for the discovery and replication cohorts and for the fixed effects meta-analysis. The red dashed line marks the effect size calculated from the fixed effects meta-analysis. P values for each cohort are shown at the right of the plot and the meta-analysis P value is also given, all P values are two-sided.

We performed tests for pair-wise interaction (see Supplementary Note) in the discovery data between all pairs of the 16 SNPs taken forward to Stage 2 replication (Supplementary Table 5), but no significant interactions (P<0.01) were found.

Imputation was carried out in the discovery data for the chromosome 6p21 and 16q24 regions, using the 1000 Genomes June 2010 CEU reference panel. In each case, rs9257809 and rs9936833 respectively remained as the strongest signal of association in each region (Supplementary Figure 4).

The lead SNP on 16q24, rs9936833, maps 24kb from the spliced, non-coding transcript LOC732275. The closest coding gene, 141kb towards the telomere, is FOXF1, a forkhead family transcription factor that acts in the hedgehog signaling pathway. FOXF1 is known to have a role in the development of the gastrointestinal tract and has been reported to cause esophageal structural alterations, especially atresia, when inactivated [25]. The region around rs9936833 contains multiple binding sites for specific transcription factors, such as FOXP2, that are known to control FOXF1 expression (assessed using ENCODE data, see URLs).

The lead SNP on 6p21, rs9257809, lies on the telomeric edge of the major histocompatibility complex (MHC) region between olfactory receptor genes OR2D12 and OR2D13. It is in strong long-range linkage disequilibrium (r2>0.6 calculated in the control data) with SNPs over 1 Mb away, including two at which Stage 2 replication was attempted, rs13211507 (Pcombined =8.77×10−9) and rs9262143 (Pcombined=2.18×10−8). When conditioning on rs9257809, no other SNP in the MHC was significant at P<10−5.

To further investigate the SNP signal in the MHC region, we took two approaches: GENECLUSTER, which is a Bayesian tree building method[26,27]; and HLA*IMP, which is a method for imputing classical HLA alleles from SNP data[28]. Both methods provided evidence of association in the discovery data for reduced risk conferred by three classical HLA alleles that are in strong LD with each other (HLA-C*07:01, HLA-A*01:01 and HLA-B*08:01), see Supplementary Table 6. However, conditional analysis suggested that rs9257809 better captures the association in our discovery data and none of these three classical HLA alleles showed an association signal in the replication data (P>0.1, Supplementary Table 6).

We used standard UK criteria, in accordance with the British Society of Gastoenterology, for diagnosis of Barrett’s Esophagus. However some countries use the American College of Gastroenterology criteria that require the presence of intestinal metaplasia for the diagnosis of Barrett’s Esophagus. To investigate this, we analysed the two replicated loci using only the subset of discovery and replication cases (86%) with histological evidence of intestinal metaplasia. Both signals remained significant, with combined evidence across discovery and all stages of replication of P<5×10−8 (Supplementary Table 7A and 7B).

We also investigated associations with the related quantitative traits of circumferential extent (C) and maximal extent (M) of the length of Barrett’s segment. In the discovery cohort, the C measurement was available for 1,744 cases, and the M measurement for 1,618 cases. In a linear regression analysis of cases, neither SNP showed evidence of association with C or M status (for rs9936833, P=0.63 and P=0.87 respectively; for rs925809, P=0.10 and P=0.09 respectively). We then extended the C and M analysis genome-wide. No SNP reached P<10−6 in the analysis of C. One SNP (rs1023313) reached P<10−6 in the analysis of M, but this association was not confirmed in Stage 1 or Stage 2 replication (see Supplementary Table 8).

There is an established sex bias in BE susceptibility, with men at greater risk than women[3,29]. The ratio of males to females is 4:1 in our case discovery data. To see whether there might be sex-specific effects of any predisposition SNPs, we performed a sex-stratified analysis for the 16 SNPs analysed in Stage 2 (Supplementary Table 9). The SNP showing the most evidence for a sex-specific effect from the combined discovery and Stage 1 and 2 replication was rs9257809. The association signal was stronger in males than females (uncorrected P=0.01 for difference of effects between sexes), corresponding to a male odds ratio of 1.38 (95%CI 1.25-1.53, Pcombined=1.71×10−10) and a female odds ratio of 1.11 (95%CI 0.95-1.30, Pcombined=0.19), see Supplementary Note for further details. This finding warrants further investigation.

Previous genome-wide association studies of common diseases or phenotypes have found evidence for a model where many common variants of small effect influence risk [30,31]. We looked for these en masse effects in Barrett’s Esophagus using two methods (see Online Methods). Firstly, taking the top K SNPs (for different values of K) in independent regions in the discovery data, we performed a sign test to see whether there was an excess (over the proportion expected under the null of 50%) of SNPs for which the effect was in the same direction in the Stage 1 replication data. Secondly a disease-score test analysis was undertaken, as described by the International Schizophrenia Consortium [30]. Both methods found evidence of an excess of SNPs that have the same risk allele in both cohorts. The strongest evidence in the sign test was for the top 1,100 SNPs, for which the sign test gave Puncorrected=2.30×10−5 (Supplementary Figure 5). From the disease-score analysis, the strongest evidence was for the top 1,710 SNPs, for which Puncorrected=7.07×10−11 (Supplementary Figure 6). Both analyses thus implicate a large number of common SNPs of small effect in susceptibility to Barrett’s Esophagus.

There is a well-established link between Barrett’s Esophagus and obesity[32,33]. To investigate whether this may in part reflect genetic effects, we repeated the sign test at 40 of the SNPs that have been found to be associated with either Body Mass Index (BMI) or Waist Hip Ratio (WHR), where genotype data or tag SNPs were available in our discovery samples [34-38]. In our discovery data, a total of 29 out of 40 BMI/WHR-associated SNPs (14 genotyped, 15 tagging, Supplementary Tables 10A and 10B) shared the same risk alleles in Barrett’s Esophagus as they did for BMI/WHR (P=6.42×10−3).

Our results provide direct evidence that Barrett’s Esophagus aetiology has a genetic component. Inference as to the underlying genes must be cautious, especially for the variant (tagged by rs9257809) in the gene-rich MHC region in which linkage disequilibrium is long-range and complex. However, the location of the other associated SNP, rs9936833, near FOXF1 suggests a role for structural factors in the esophagus and stomach as a disease-predisposing factor, consistent with the evidence that changes such as hiatus hernia are known to be strongly associated with Barrett’s Esophagus. We also found evidence to show that body weight SNPs are more likely than by chance to show effects in the same direction in Barrett’s Esophagus, suggesting that genetic effects may in part underpin the epidemiological observation that BMI is a risk factor for Barrett’s Esophagus [39]. Given that Barrett’s Esophagus is an accepted status as a precursor lesion, the SNPs that we have identified could also be de facto risk factors for esophageal adenocarcinoma and may give clues as to the biology of both of these important phenotypes.

Online Methods

Samples

Cases from Discovery, Stages 1 and 2 replication, and Stage 3 UK and Irish

For the discovery, we ascertained cases of histologically-confirmed Barrett’s Esophagus through the United Kingdom-based ASPECT clinical trial of proton pump-inhibitor (esomeprazole) and aspirin as preventive agents for progression of Barrett’s Esophagus to EAC[20]. UK, Irish and Dutch replication cases were from the Chemoprevention of Premalignant Intestinal Neoplasia (ChOPIN) genetic study and the Esophageal Adenocarcinoma GenE (EAGLE) consortium[9]. Replication cases were diagnosed with Barrett’s Esophagus with lengths of at least 1cm (CIMI) circumferential Barrett’s Esophagus or at least a 2cm tongue (C0M2) according to the Prague criteria[40]. Case collection was in accordance with the British Society of Gastoenterology criteria[41], the standard practice for collaborating Histopathologists in the UK and much of Europe. We found that 90% of our discovery samples (for which full clinico-pathological data were available) had evidence of intestinal metaplasia and therefore also met the American College of Gastroenterology criteria that are widely used in the USA[42] . For full details of the ethnicity, age and sex distributions and Prague criteria measurements of the cases see Supplementary Table 1.

Discovery

The full data set comprised of 1,991 cases and 5,667 controls. After QC, 1,852 cases and 5,172 controls were analysed. Controls were taken from the WTCCC2 set, made up of samples from the 1958 British Birth Cohort (58C) and the National Blood Service collection (UKBS). Samples were genotyped at the Wellcome Trust Sanger Institute (WTSI), cases on the Illumina Human660W-Quad array, and controls on the Illumina custom Human 1.2M-Duo. The primary analysis was performed on the overlapping set of SNPs.

Stage 1

After QC, the UK replication totalled 1,105 cases and 6,819 controls. The controls were from the PoBI cohort (2,578)[24] and 58C (4,241) samples that were not genotyped in the discovery phase. The case and 58C control samples were genotyped on the Illumina Immunochip and the PoBI samples were genotyped on the Illumina custom Human 1.2M-Duo array. The Immunochip is a custom-designed chip containing 196,524 SNPs in total, of which ~2,400 were selected on the basis of our discovery GWAS study.

Stage 2

The Dutch replication cohort consisted of 473 cases and 1,780 controls. These samples were all genotyped on the Illumina Immunochip but in two separate locations; the case samples were genotyped at WTSI and the control samples were genotyped as described in a previous report[43].

See Supplementary Note for information on DNA sample preparation.

Circumferential and Maximal Extent Phenotypes

Length of the Barrett’s segment was available for a subset of discovery and replication phase samples. Where baseline measurements were not available, the earliest measurement taken after baseline was used. A small number of cases were excluded on the basis of reporting errors (if C >M or if either value exceeded 25cm). Of the discovery phase individuals after quality control, 1,744 had C measurements and 1,618 had M measurements, C mean=4.05 (range 0-22); M mean=4.60 (range 1-24). M measurements were available for 1,015 of the Stage 1 replication (M mean=4.66 (range 1-23)) and for 240 of the Stage 2 replication (M mean=4.44 (range 1-15)). Both C and M phenotypes were square-root transformed prior to analysis, to improve the fit of the linear regression model.

Stage 3

Irish Replication

The Irish replication cohort consisted of 245 cases and 473 controls. Cases were provided by St James’s Hospital and Mater Misericordiae University Hospital Dublin as part of EAGLE. Controls were provided by Trinity Biobank. 168 cases were genotyped on the Illumina Immunochip at WTSI. rs9257809 and rs9936833 were genotyped in 77 cases and all controls using competitive allele-specific PCR KASPar chemistry (KBiosciences Ltd, Hertfordshire, UK). Primers, probes and conditions used are available on request. Genotyping quality control was tested using duplicate DNA samples within studies and SNP assays, together with direct sequencing of subsets of samples to confirm genotyping accuracy. For all SNPs, >99% concordant results were obtained.

UK Replication 2

1,765 cases were ascertained using the diagnostic criteria and sampling from ASPECT as described above for discovery. 1586 controls were collected as part of the Colorectal Tumour Gene Identification (CoRGI) consortium[44]. Controls were spouses or partners unaffected by cancer and without a personal family history (to 2nd degree relative level) of colorectal neoplasia. All were of white UK ethnic origin, 45% male; mean age 45.1 years, SD±15.9. All samples were genotyped using KASPar competitive allele-specific PCR as described above.

BEACON Replication

2,398 cases and 2,167 controls were analyzed. Samples were collected as part of a GWAS study (BEAGESS) undertaken by the BEACON collaboration. Samples were collected from sites in Australia (cases n=325, controls n=561), Europe (England, Ireland, Sweden; cases n=363, controls n=333), and North America (Canada, United States, cases n=1,710, controls n=1,273). Samples were genotyped at the Fred Hutchinson Cancer Research Center (FHCRC) on the Illumina Omni1M Quad.

Quality control

Samples

As previously described[21,45], we identified and removed samples whose genome-wide patterns of diversity differed from those of the collection at large, interpreting them as likely to be due to biases or artefacts. See Supplementary Note for further details. Following sample quality control our final discovery dataset consisted of 1,852 cases and 5,172 controls (Supplementary Table 3).

SNPs

For all arrays, normalised probe intensities were exported using the BeadStudio program and genotypes were called at the WTSI using Illuminus[46]. SNPs were excluded from analysis if in any of the data sets (58C, UKBS or cases) they had: a very low minor allele frequency (defined as <0.01%); extreme departures from Hardy-Weinberg equilibrium (P<10−20); showed a strong plate effect (P<10−6). SNPs were also excluded if the observed statistical (Fisher) information about the allele frequency was less than 98% of the information contained in a hypothetical sample of the same size and expected MAF but with no missing data. 45 SNPs were removed following visual inspection of cluster plots. In total 521,744 autosomal SNPs were available for analysis after quality control.

To confirm genotyping accuracy of the different platforms used in the study, 5% of the UK, Irish and Dutch samples typed on each platform were re-genotyped at rs9257809 and rs9936833 using competitive allele-specific PCR KASPar. Concordance was >99% (Supplementary Table 11) suggesting genotyping robustness across platforms.

HLA Imputation

Classical HLA alleles were imputed using HLA*IMP [28]. Further details of this can be found in Supplementary Note.

Statistical analysis

Genome-wide case-control analysis was performed using frequentist tests, under a missing data logistic regression model, as implemented in SNPTEST. Unless otherwise stated, we assumed a multiplicative model for allelic risk by encoding the genotypes at each SNP as a discrete explanatory variable with an indicator of case status as the binary response and the first principal component as a covariate (see Supplementary Note). Quantitative C and M measurements were analysed using frequentist tests under a missing data linear regression model, as implemented in SNPTEST. To combine the evidence of association across the discovery and replication datasets we conducted an inverse-variance weighted fixed effect meta-analysis in the statistical package R (see Supplementary Note). To test for interactions (see Supplementary Note), between SNPs, or between a SNP and sex, and to compare models which include additional SNPs or classical HLA alleles as predictors, we used logistic regression models implemented in R. These analyses used thresholded (posterior probability > 0.9) genotype calls.

SNP Imputation was performed using IMPUTE2[47], which adopts a two-stage approach using both a haploid reference panel and a diploid reference panel.

BEACON data was analysed under an additive logistic regression model including the first four principal components as covariates (see Supplementary Note). Genomic inflation λ was 1.037.

En Masse analysis was carried out on the discovery and Stage 1 data. In order to reduce possible population structure (such analyses are sensitive to this), we restricted the Stage 1 control set to the 58C individuals. SNPs with MAF > 0.01 which were genotyped in both the discovery (Illumina 670K and Illumina custom Human 1.2M-Duo) and the replication (Illumina Immunochip) were pruned to remove strong linkage disequilibrium. This was done by ranking the SNPs by Bayes factor calculated under an additive model in SNPTEST, and successively selecting SNPs from the top so that they were at least 0.125cM plus 25kb away from any SNPs that had already been selected. We obtained 7,673 SNPs from a total of 28,972 (after quality control) that were typed in discovery and UK Immunochip data. For the K SNPs showing the strongest signal of association, the sign test compares the direction of effect of each SNP in the discovery and replication samples. Using a likelihood-ratio test we compared the null model where the probability of the same direction of effect is assumed to be a half, to a model where the probability is not a half (two sided).

The disease-score test aims to measure indirectly the collective effect of many weakly associated alleles. We determined the risk allele and odds ratio for each pruned SNP from the discovery data as described above. Then, we used the top K SNPs to calculate the “score” for each individual in the replication data as the number of risk alleles carried by each individual weighted by the log of the odds ratio estimated from the discovery data. Under the null hypothesis, the risk alleles and odds ratios in the discovery and replication samples are independent. We tested a logistic regression model of disease status on the score in the replication data, conditioning on the first principal component, to control for population structure, and the number of missing genotypes (called with maximum probability < 0.9), to control for potential differences in genotyping rate, as covariates.