A genome-wide association study of marginal zone lymphoma shows association to the HLA region

Marginal zone lymphoma (MZL) is the third most common subtype of B-cell non-Hodgkin lymphoma. Here we perform a two-stage GWAS of 1,281 MZL cases and 7,127 controls of European ancestry and identify two independent loci near BTNL2 (rs9461741, P=3.95 × 10−15) and HLA-B (rs2922994, P=2.43 × 10−9) in the HLA region significantly associated with MZL risk. This is the first evidence that genetic variation in the major histocompatibility complex influences MZL susceptibility.

M arginal zone lymphoma (MZL) encompasses a group of lymphomas that originate from marginal zone B cells present in extranodal tissue and lymph nodes. Three subtypes of MZL have been defined, extranodal MZL of mucosaassociated lymphoid tissue (MALT), splenic MZL and nodal MZL, which together account for 7-12% of all non-Hodgkin lymphoma (NHL) cases. Geographic differences in incidence have been observed 1 , and inflammation, immune system dysregulation and infectious agents, such as Helicobacter pylori, have been implicated particularly for the gastric MALT subtype 2 , but little else is known of MZL aetiology.
Here we perform the first two-stage, subtype-specific genomewide association study (GWAS) of MZL and identify two independent single-nucleotide polymorphisms (SNPs) within the HLA region associated with MZL risk. Together with recent studies on other common subtypes of NHL, these results point to shared susceptibility loci for lymphoma in the HLA region.

Results
Stage 1 MZL GWAS. As part of a larger NHL GWAS, 890 MZL cases and 2,854 controls from 22 studies in the United States and Europe (Supplementary Table 1) were genotyped using the Illumina OmniExpress array. Genotype data from the Illumina Omni2.5 was also available for 3,536 controls from three of the 22 studies 3 . After applying rigorous quality control filters (Supplementary Table 2, Methods), data for 611,856 SNPs with minor allele frequency (MAF)41% in 825 cases and 6,221 controls of European ancestry ( Supplementary Fig. 1) remained for the stage 1 analysis (Supplementary Table 3). To discover variants associated with risk, logistic regression analysis was performed on these SNPs adjusting for age, gender and three significant eigenvectors computed using principal components analysis ( Supplementary Fig. 2, Methods). Examination of the quantile-quantile (Q-Q) plot ( Supplementary Fig. 3) showed minimal detectable evidence for population substructure (l ¼ 1.01) with some excess of small P values. A Manhattan plot revealed association signals at the HLA region (Supplementary Fig. 4; 6p21.33:31,061,211-32,620,572) on chromosome 6 reaching genome-wide significance. Removal of all SNPs in the HLA region resulted in an attenuation of the excess of small P values observed in the Q-Q plot, although some excess still remained. To further explore associations within the HLA region and identify other regions potentially associated with risk, common SNPs available in the 1000 Genomes project data release 3 were imputed (Methods).
Stage 2 genotyping. Ten SNPs in promising loci with Pr7.5 Â 10 À 6 in the stage 1 discovery were selected for replication (stage 2) in an additional 456 cases and 906 controls of European ancestry (Supplementary Tables 1 and 3). Of the SNPs selected for replication, two SNPs were directly genotyped on the OmniExpress, while the remaining eight were imputed with high accuracy (median info score ¼ 0.99) in stage 1 (Supplementary  Table 4). Replication was carried out using Taqman genotyping. In the combined meta-analysis of 1,281 cases and 7,127 controls, we identified two distinct loci (Table 1, Fig. 1, Supplementary  Table 4) at chromosomes 6p21.32 and 6p21.33 that reached the threshold of genome-wide statistical significance (Po5 Â 10 À 8 ). These are rs9461741 in the butyrophilin-like 2 (MHC class II associated) (BTNL2) gene at 6p21.32 in HLA class II (P ¼ 3.95 Â 10 À 15 , odds ratio (OR) ¼ 2.66, confidence interval (CI) ¼ 2.08-3.39) and rs2922994 at 6p21.33 in HLA class I (P ¼ 2.43 Â 10 À 9 , OR ¼ 1.64, CI ¼ 1.39-1.92). These two SNPs were weakly correlated (r 2 ¼ 0.008 in 1000 Genomes CEU population), and when both were included in the same statistical model, both SNPs remained strongly associated with MZL risk (rs9461741, P ¼ 2.09 Â 10 À 15 ; rs2922994, P ¼ 6.03 Â 10 À 10 ), suggesting that the two SNPs are independent. Both SNPs were weakly correlated with other SNPs in the HLA region previously reported to be associated with other NHL subtypes or Hodgkin lymphoma (r 2 o0.14 for all pairwise comparisons). None of the previously reported SNPs were significantly associated with MZL risk after adjustment for multiple testing (Po0.0025) in our study, suggesting the associations are subtypespecific (Supplementary Table 5). Another SNP rs7750641 (P ¼ 3.34 Â 10 À 8 ; Supplementary Table 4) in strong linkage disequilibrium (LD) with rs2922994 (r 2 ¼ 0.85) also showed promising association with MZL risk. rs7750641 is a missense variant in transcription factor 19 (TCF19), which encodes a DNA-binding protein implicated in the transcription of genes during the G1-S transition in the cell cycle 4 . The non-HLA SNPs genotyped in stage 2 were not associated with MZL risk (Supplementary Table 4).
MALT versus non-MALT. Heterogeneity between the largest subtype of MZL, namely MALT and other subtypes grouped as non-MALT, was evaluated for the MZL associated SNPs (Supplementary Table 7). The effects were slightly stronger for MALT, but no evidence for substantial heterogeneity was observed (P heterogeneity Z0.05). Studies have suggested that H. pylori infection is a risk factor for gastric MZL 2 . An examination of SNPs previously suggested to be associated with H. pylori infection in independent studies 6 did not reveal any significant association with the combined MZL or the MALT subtype in this study (Supplementary Table 8). Toll-like receptors (TLR) are considered strong candidates in mediating inflammatory immune response to pathogenic insults. A previous study reported 7 a nominally significant association with rs4833103 in the TLR10-TLR1-TLR6 region with MZL risk. After excluding the cases and controls in the previous report 7 , we found little additional support for this locus (MZL: Secondary functional analyses. To gain additional insight into potential biological mechanisms, expression quantitative trait loci (eQTL) analyses were performed in two datasets consisting of lymphoblastoid cell lines (Methods). Significant associations were seen for rs2922994 and rs7750641with HLA-B and HLA-C (Supplementary Table 9) while suggestive associations (false discovery rate, FDRr0.05) for correlated SNPs of rs2922994 (r 2 40.8) in HLA class I and RNF5 (Supplementary Table 10) were observed. No significant eQTL association was observed for rs9461741 or other correlated HLA class II SNPs. Chromatin state analysis (Methods) using ENCODE data revealed correlated SNPs of rs2922994 showed a chromatin state consistent with the prediction for an active promoter (rs3094005) or satisfied the state of a weak promoter (rs2844577) in the lymphoblastoid cell line   Fig. 5). GM12878 is the only lymphoblastoid cell line from which high-quality whole-genome annotation data for chromatin state is readily available. Analyses were also performed with HaploReg (Supplementary Table 11) and RegulomeDB (Supplementary Table 12) that showed overlap of the SNPs with functional motifs, suggesting plausible roles in gene regulatory processes.

Discussion
The most statistically significant SNP associated with MZL, rs9461741, is located in HLA class II in the intron between exons 3 and 4 of the BTNL2 gene. BTNL2 is highly expressed in lymphoid tissues 8 and has close homology to the B7 costimulatory molecules, which initiate lymphocyte activation as part of antigen presentation. Evidence is consistent with BTNL2 acting as a negative regulator of T-cell proliferation and cytokine production 8,9 and attenuating T-cell-mediated responses in the gut 10 . We were unable to statistically differentiate the effects of rs9461741 from HLA-DRB1*01:02 and, thus, our observed association could be due to HLA-DRB1. HLA-DRB1 has been shown to be associated with other autoimmune diseases, including rheumatoid arthritis 11 and selective IgA deficiency 12 .
Similarly, rs2922994 is located 11 kb upstream of HLA-B, which is known to play a critical role in the immune system by presenting peptides derived from the endoplasmic reticulum lumen. rs7750641, a missense variant in TCF19, was previously associated with pleiotropic effects on blood-based phenotypes 13 and is highly expressed in germinal center cells and up-regulated in human pro-B and pre-B cells 14 . Autoimmune diseases, such as Sjögren's syndrome and systemic lupus erythematosus, are established risk factors for developing MZL, with the strongest associations seen between Sjögren's syndrome and the MALT subtype 15 . Of note, SNPs in HLA-B and the classical alleles HLA-DRB1*01:02 are strongly associated with Sjögren's syndrome 16 , while HLA-DRB1*03 has been associated with rheumatoid arthritis 17 . The multiple independent associations in the HLA region and their localization to known functional autoimmune and B-cell genes suggest possible shared genetic effects that span both lymphoid cancers and autoimmune diseases. Chronic autoimmune stimulation leading to overactivity and defective apoptosis of B cells, and secondary inflammation events triggered by genetic and environmental factors are biological mechanisms that may contribute to the pathogenesis of MZL. We have performed the largest GWAS of MZL to date and identified two independent SNPs within the HLA region that are robustly associated with the risk of MZL. In addition to the known diversity in etiology and pathology, there is mounting evidence of genetic heterogeneity across the NHL subtypes of lymphoma. However, the HLA region appears to be commonly associated with multiple major subtypes, such as MZL, CLL 18 , DLBCL 19 and FL [20][21][22][23] . Further studies are needed to identify biological mechanisms underlying these relationships and advance our knowledge regarding their interactions with associated environmental factors that may modulate disease risks. Cases were ascertained from cancer registries, clinics or hospitals or through self-report verified by medical and pathology reports. To determine the NHL subtype, phenotype data for all NHL cases were reviewed centrally at the International Lymphoma Epidemiology Consortium (InterLymph) Data Coordinating Center and harmonized using the hierarchical classification proposed by the InterLymph Pathology Working Group 24,25 based on the World Health Organization (WHO) classification 26 .

Methods
Genotyping and quality control. All MZL cases with sufficient DNA (n ¼ 890) and a subset of controls (n ¼ 2,854) frequency matched by age, sex and study to the entire group of NHL cases, along with 4% quality control duplicates, were genotyped on the Illumina OmniExpress at the NCI Core Genotyping Resource (CGR). Genotypes were called using Illumina GenomeStudio software, and quality control duplicates showed 499% concordance. Monomorphic SNPs and SNPs with a call rate of o95% were excluded. Samples with a call rate of r93%, mean heterozygosity o0.25 or 40.33 based on the autosomal SNPs or gender discordance (45% heterozygosity on the X chromosome for males and o20% heterozygosity on the X chromosome for females) were excluded. Furthermore, unexpected duplicates (499.9% concordance) and first-degree relatives based on identity by descent sharing with Pi-hat 40.40 were excluded. Ancestry was assessed using the Genotyping Library and Utilities (GLU-http://code.google.com/p/glu-genetics/) struct.admix module based on the method by Pritchard et al. 27 and participants with o80% European ancestry were excluded (Supplementary Fig. 1). After exclusions, 825 cases and 2,685 controls remained (Supplementary Table 2). Genotype data previously generated on the Illumina Omni2.5 from an additional 3,536 controls from three of the 22 studies (ATBC, CPS-II and PLCO) were also included 3 , resulting in a total of 825 cases and 6,221 controls for the stage 1 analysis (Supplementary Table 3). Of these additional 3,536 controls, 703 (B235 from each study) were selected to be representative of their cohort and cancer free 3 , while the remainder were cancer-free controls from an unpublished study of prostate cancer in the PLCO. SNPs with call rate o95%, with Hardy-Weinberg equilibrium Po1 Â 10 À 6 , or with a MAF o1% were excluded from analysis, leaving 611,856 SNPs for analysis. To evaluate population substructure, a principal components analysis was performed using the Genotyping Library and Utilities (GLU), version 1.0, struct.pca module, which is similar to EIGENSTRAT 28 -http://genepath. med.harvard.edu/Breich/Software.htm. Plots of the first five principal components are shown in Supplementary Fig. 2. Genomic inflation factor was computed prior (l ¼ 1.014) and after removal of SNPs in the HLA loci (l ¼ 1.010). Association testing was conducted assuming a log-additive genetic model, adjusting for age, sex and three significant principal components. All data analyses and management were conducted using GLU.

Imputation of variants.
To more comprehensively evaluate the genome for SNPs associated with MZL, SNPs in the stage 1 discovery GWAS were imputed using IMPUTE2 (ref. 29)-http://mathgen.stats.ox.ac.uk/impute/impute_v2.html and the 1000 Genomes Project (1kGP-http://www.1000genomes.org/) version 3 data 29,30 . SNPs with a MAF o1% or information quality score (info) o0.3 were excluded from analysis, leaving 8,478,065 SNPs for association testing. Association testing on the imputed data was conducted using SNPTEST 31 -https://mathgen.stats.ox.ac. uk/genetics_software/snptest/snptest.html version 2, assuming dosages for the genotypes and adjusting for age, sex and three significant principal components. In a null model for MZL risk, the three eigenvectors EV1, EV3 and EV8 were nominally associated with MZL risk and hence were included to account for potential population stratification. Heterogeneity between MZL subtypes was assessed using a case-case comparison, adjusting for age, sex and significant principal components.
Stage 2 replication of SNPs from the GWAS. After ranking the SNPs by P value and LD filtering (r 2 o0.05), 10 SNPs from the most promising loci identified from stage 1 after imputation with Po7.5 Â 10 À 6 were taken forward for de novo replication in an additional 456 cases and 906 controls (Supplementary Tables 1  and 4). Wherever possible, we selected either the best directly genotyped SNP or ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6751 the most significant imputed SNP for the locus. Only imputed SNPs with an information score 40.8 were considered for replication. Only SNPs with MAF 41% were selected for replication, and no SNPs were taken forward for replication in regions where they appeared as singletons or obvious artifacts. For the HLA region, we selected one additional SNP (rs7750641) that was highly correlated with rs2922994 for additional confirmation. Genotyping was conducted using custom TaqMan genotyping assays (Applied Biosystems) validated at the NCI Core Genotyping Resource. Genotyping was done at four centres. HapMap control samples genotyped across two centres yielded 100% concordance as did blind duplicates (B5% of total samples). Due to the small number of samples, the MD Anderson, Mayo and NCI replication studies were pooled together for association testing; however, MSKCC samples were analysed separately to account for the available information on Ashkenazi ancestry. Association results were adjusted for age and gender and study site in the pooled analysis. The results from the stage 1 and stage 2 studies were then combined using a fixed effect meta-analysis method with inverse variance weighting based on the estimates and s.e. from each study. Heterogeneity in the effect estimates across studies was assessed using Cochran's Q statistic and estimating the I 2 statistic. For all SNPs that reached genome-wide significance in Table 1, no substantial heterogeneity was observed among the studies (P heterogeneity Z0.1 for all SNPs, Supplementary Table 4).
Technical validation of imputed SNPs. Genotyping was conducted using custom TaqMan genotyping assays (Applied Biosystems) at the NCI Cancer Genomics Research Laboratory on a set of 470 individuals included in the stage 1 MZL GWAS. The allelic dosage r 2 was calculated between the imputed genotypes and the technical validation done using assayed genotypes which showed that both SNPs were imputed with high accuracy (INFO Z0.99) and a high correlation (r 2 Z0.99) between dosage imputation and genotypes obtained by Taqman assays.
HLA imputation and analysis. To determine if specific coding variants within HLA genes contributed to the diverse association signals, we imputed the classical HLA alleles (A, B, C, DQA1, DQB1, DRB1) and coding variants across the HLA region (chr6:20-40 Mb) using SNP2HLA 5 -http://www.broadinstitute.org/mpg/ snp2hla/. The imputation was based on a reference panel from the Type 1 Diabetes Genetics Consortium (T1DGC) consisting of genotype data from 5,225 individuals of European descent who were typed for HLA-A, B, C, DRB1, DQA1, DQB1, DPB1, DPA1 4-digit alleles. Imputation accuracy of HLA alleles was assessed by comparing HLA alleles to the HLA sequencing data on a subset of samples from the NCI 32 . The concordance rates obtained were 97.32, 98.5, 98.14 and 97.49% for HLA-A, B, C and DRB1, respectively, in the NCI GWAS suggesting robust performance of the imputation method. Due to the limited number of SNPs (7,253) in the T1DGC reference set, imputation of HLA SNPs was conducted with IMPUTE2 and the 1kGP reference set as described above. A total of 68,488 SNPs, 201 classical HLA alleles (two-and four-digit resolution) and 1,038 AA markers including 103 AA positions that were 'multi-allelic' with three to six different residues present at each position, were successfully imputed (info score 40.3 for SNPs or r 2 40.3 for alleles and AAs) and available for downstream analysis. Multiallelic markers were analysed as binary markers (for example, allele present or absent) and a meta-analysis was conducted where we tested SNPs, HLA alleles and AAs across the HLA region for association with MZL using PLINK 33 or SNPTEST 31 as described above. eQTL analysis. We conducted an eQTL analysis using two independent datasets: childhood asthma 34 and HapMap 35 . As described previously 34 for the childhood asthma data set 35 , peripheral blood lymphocytes were transformed into lymphoblastoid cell lines for 830 parents and offspring from 206 families of European ancestry. Data from 405 children were used for the analysis as follows: using extracted RNA, gene expression was assessed with the Affymetrix HG-U133 Plus 2.0 chip. Genotyping was conducted using the Illumina Human-1 Beadchip and Illumina HumanHap300K Beadchip, and imputation performed using data from 1kGP. All SNPs selected for replication were tested for cis associations (defined as gene transcripts within 1 Mb), assuming an additive genetic model, adjusting for non-genetic effects in the gene expression value. Association testing was conducted using a variance component-based score test 36 in MERLIN 37 , which accounts for the correlation between siblings. To gain insight into the relative importance of associations with our SNPs compared with other SNPs in the region, we also conducted conditional analyses, in which both the MZL SNP and the most significant SNP for the particular gene transcript (that is, peak SNP) were included in the same model. Only cis associations that reached Po6.8 Â 10 À 5 , which corresponds to a FDR of 1% are reported (Supplementary Table 9).
The HapMap data set consisted of a publicly available RNAseq data set 35 from transformed lymphoblastoid cell lines from 41 CEPH Utah residents with ancestry from northern and western Europe (HapMap-CEU) samples available from the Gene Expression Omnibus repository (http://www.ncbi.nlm.nih.gov/geo) under accession number GSE16921. In this data set, we examined the association between the two reported SNPs in the HLA region, rs2922994 and rs9461741, as well as all SNPs in LD (r 2 40.8 in HapMap-CEU release 28) and expression levels of probes within 1 Mb of the SNPs. As rs9461741 was not genotyped in HapMap, we selected rs7742033 as a proxy as it was the strongest linked SNP available in HapMap (r 2 ¼ 0.49 in 1kGP-CEU). Genotyping data for these HapMap-CEU individuals were directly downloaded from HapMap (www.hapmap.org). Correlation between expression and genotype for each SNP-probe pair was tested using the Spearman's rank correlation test with t-distribution approximation and estimated with respect to the minor allele in HapMap-CEU. P values were adjusted using the Benjamini-Hochberg FDR correction and eQTLs were considered significant at an FDRo0.05 (Supplementary Table 10).
Bioinformatics ENCODE and chromatin state dynamics. To assess chromatin state dynamics, we used Chromos 38 , which has precomputed data from ENCODE on nine cell types using Chip-Seq experiments 39 . These consist of B-lymphoblastoid cells (GM12878), hepatocellular carcinoma cells (HepG2), embryonic stem cells, erythrocytic leukaemia cells (hK562), umbilical vein endothelial cells, skeletal muscle myoblasts, normal lung fibroblasts, normal epidermal keratinocytes and mammary epithelial cells. These precomputed data have genome-segmentation performed using a multivariate hidden Markov model to reduce the combinatorial space to a set of interpretable chromatin states. The output from Chromos lists data into 15 chromatin states corresponding to repressed, poised and active promoters, strong and weak enhancers, putative insulators, transcribed regions and large-scale repressed and inactive domains ( Supplementary Fig. 5).