Inherited variation in immune genes and pathways and glioblastoma risk.

To determine whether inherited variations in immune function single-nucleotide polymorphisms (SNPs), genes or pathways affect glioblastoma risk, we analyzed data from recent genome-wide association studies in conjunction with predefined immune function genes and pathways. Gene and pathway analyses were conducted on two independent data sets using 6629 SNPs in 911 genes on 17 immune pathways from 525 glioblastoma cases and 602 controls from the University of California, San Francisco (UCSF) and a subset of 6029 SNPs in 893 genes from 531 cases and 1782 controls from MD Anderson (MDA). To further assess consistency of SNP-level associations, we also compared data from the UK (266 cases and 2482 controls) and the Mayo Clinic (114 cases and 111 controls). Although three correlated epidermal growth factor receptor (EGFR) SNPs were consistently associated with glioblastoma in all four data sets (Mantel-Haenzel P values = 1 × 10⁻⁵ to 4 × 10⁻³), independent replication is required as genome-wide significance was not attained. In gene-level analyses, eight immune function genes were significantly (minP < 0.05) associated with glioblastoma; the IL-2RA (CD25) cytokine gene had the smallest minP values in both UCSF (minP = 0.01) and MDA (minP = 0.001) data sets. The IL-2RA receptor is found on the surface of regulatory T cells potentially contributing to immunosuppression characteristic of the glioblastoma microenvironment. In pathway correlation analyses, cytokine signaling and adhesion-extravasation-migration pathways showed similar associations with glioblastoma risk in both MDA and UCSF data sets. Our findings represent the first systematic description of immune genes and pathways that characterize glioblastoma risk.

Introduction

Glioblastoma is an immunosuppressive tumor characterized by a median survival time of only 14 months (1). Although little is known about its etiology, there is evidence that this tumor's immune microenvironment can suppress or promote its development (2). In addition, epidemiological studies have consistently found an inverse association between self-reported allergies and glioma risk (3) and two studies show reduced risk of glioma among non-steroidal anti-inflammatory drug users (4,5). Results of analyses of associations between allergy-related genetic variants and glioma risk are mixed; however, these studies have been restricted to a relatively small number of single-nucleotide polymorphisms (SNPs) or haplotypes (1,6–10). A study of 1397 innate immune system SNPs also suggested that innate immunity may affect glioma risk (11), but none of the results were statistically significant after adjustment for multiple comparisons.

The advantage of single-locus analysis is that because SNPs represent relatively small areas of the genome, results of SNP association studies are less likely to be obscured by misclassification than are results of gene- or pathway-based studies. Furthermore, in contrast to candidate gene association studies, in which candidate markers are selected on the basis of strong prior biological hypotheses, genome-wide analysis (GWA) studies scan markers across the entire genome ‘agnostically’ to identify previously unsuspected risk loci, as was recently demonstrated by the success of breakthrough studies about glioma. Both of these studies (12,13) identified glioma risk loci in or near CDKN2A/B 9p21, 5p15.33 (TERT) and 20q13.33 (RTEL1).

Despite the unparalleled power of GWA studies in generating novel biological hypotheses, there are a few limitations inherent in their design and interpretation. First, GWA studies are specifically designed to identify associations characterized by small p values, which are the result of either strong SNP–disease associations or moderate associations characterized by little variation. Therefore, genetic variants that confer moderate but variable risks may be missed after adjustment for multiple testing, which reduces statistical power. Second, due to etiologic heterogeneity, SNPs identified in a discovery study and confirmed in a replication study are frequently not those most strongly associated with the disease (14,15). Third, gene–gene interactions, which may exert powerful effects on disease risk, are not estimated. These limitations may pose difficulty in validating GWA results and point to potential problems in pursuing only those single-marker associations with the smallest p values.

In contrast to GWA SNP studies, the gene- and pathway-based approaches take advantage of prior knowledge from previous research. Therefore, they can be potentially more powerful (have a higher probability of identifying associations that are present) and at the same time have fewer false-positive findings (16) than studies of the same size that are not based on previous research. In addition, the gene-based approach is less susceptible to erroneous findings due to genetic differences among populations. Given these advantages, gene-based studies are more likely to be replicated than are similar studies based on the single-SNP approach. Gene- and pathway-based analyses also allow detection of additional genes that contribute to susceptibility of complex diseases that might be missed in single-SNP analysis. One such scenario could occur when multiple SNPs in a gene or pathway are important but have small individual effects.

In the present investigation, we analyzed data from recent GWA studies to examine inherited variation in immune function SNPs, genes and pathways using previously published immune-related gene and pathway definitions (17). Various statistical methods have been developed for joint analyses of SNPs in a gene or region to identify causal variants. One such well-known method is Fisher's product test, which forms the product of all p values in a gene or region, to assess their joint effect. Similarly, one can form the product of the k most significant p values only (18) or the product of all p values at less than a preset threshold (19,20). These variations on Fisher's product test have proven more powerful than the original test in simulation studies (18). Nevertheless, one of the limits of these tests is that they require arbitrary selection of a truncation point. Another method, the permutation-based minP approach, combines all single-locus tests in a gene into a single test statistic, the minimum p value, and the empirical distribution of this test statistic is then obtained using standard permutation methods (21). Yet a third category of approaches employs multivariate regression to simultaneously model all SNPs in a gene (22).

In the current study, we used the permutation-based minP approach (21) for the gene-level analysis. The minP approach has previously been successfully applied in studying lung and biliary tract cancer risks (23,24). A comprehensive simulation study conducted by Chapman and Whittaker (25) shows that the minP approach performs well over a range of scenarios: on tag-SNP or whole-genome SNP panels, and among low-linkage disequilibrium (LD) or high-LD loci, compared with Fisher's product test and a number of multivariate methods. Given its simplicity and good performance, the minP approach was recommended as the default choice of test statistic for gene- or pathway-level association summaries (25).

We used University of California, San Francisco (UCSF) and MD Anderson (MDA) GWA data sets for SNP-, gene- and pathway-level analyses and then inspected SNP-level data from the UK and Mayo GWA data sets to determine whether SNP-level results were consistent as all four GWA data sets had not been previously used together for SNP-level comparisons (12,13).

Materials and methods

Since current policies regarding sharing large-scale genotyping data vary by country, institution and investigator, data from each site (UCSF, MDA, Mayo and UK) were analyzed separately and results forwarded to UCSF for comparisons across sites. Each site had approval from their local human subjects review board for their study.

Study participants

For a complete description of UCSF subject selection, see Felini et al. (26) and Wrensch et al. (13). Briefly, participants included 525 adults with newly diagnosed histologically verified glioblastoma patients who were identified either through the San Francisco regional population-based registry's rapid case ascertainment program or through the UCSF Neuro-oncology Clinic between 1997 and 2006. Controls from the San Francisco Bay area were identified using random digit dialing and frequency matched to population-based cases on age, sex and self-identified race and ethnicity. The GWA study was limited to people who self-reported being white and provided blood.

The MDA study was based on 531 newly diagnosed histologically verified adult glioblastoma cases (326 males and 205 females) ascertained through the MDA Cancer Center, Texas, between 1990 and 2008 (see ref. 12 for details). Individuals from the Cancer Genetic Markers of Susceptibility Study served as the 1782 controls (15,27).

The UK study included 266 glioblastoma cases ascertained through the INTERPHONE Study (28). Briefly, the INTERPHONE Study was an international multicenter case–control study of primary brain tumors coordinated by the International Agency for Research on Cancer, with material collected between September 2000 and February 2004. UK individuals newly diagnosed with glioblastoma were identified using records from neurosurgery, neuropathology, oncology and neurology centers in the Thames regions of Southeast England and the Northern UK including central Scotland, the West Midlands, West Yorkshire and the Trent area. Cases with previous brain tumors were excluded. To minimize population stratification, cases with self-reported non-western-European ancestry were excluded from the present study. Individuals from the 1958 Birth Cohort served as the source of controls (29).

The Mayo Clinic study included 114 glioblastoma patients newly diagnosed between 2005 and 2008. Cases were identified within 24 h of diagnosis, except for those who were initially diagnosed elsewhere and later had their diagnosis verified at the Mayo Clinic. Pathologic diagnosis was confirmed by review of the primary surgical material for all cases by two Mayo Clinic neuropathologists based on surgically resected material. The control group consisted of consented individuals who had a general medical exam at the Mayo Clinic. Matching variables were sex, date of birth (within two and one half years), race and residence. Ninety-eight percent of the cases described their race as ‘White’. Geographic region of residence was matched in three zones based on the distance to the Mayo Clinic Rochester: Olmsted County; the rest of Minnesota, Wisconsin, Iowa, North Dakota and South Dakota and the rest of the USA and Canada. Excluded were individuals under the age of 18 years and those with a history of brain tumor. The Mayo Clinic case and control enrollment research protocol was approved by the Mayo Clinic's Institutional Review Board (13).

Genotyping

Genotyping for UCSF subjects was conducted by deCODE genetics using Illumina's HumanCNV370-duo BeadChip. Mayo DNA was genotyped using Illumina 610Quad SNP arrays according to the manufacturer's recommendations. For complete details of methods and quality control procedures that were used for UCSF and Mayo genotyping, see Wrensch et al. (13). Genotyping MDA and UK subjects was conducted by Illumina Service laboratory using the Illumina Infinium Human610-Quad BeadChips according to Illumina protocols. Further details of genotyping procedures are available in Shete et al. (12).

Selection of allergy and immune function pathways, genes and SNPs

Loza et al. (17) identified 17 immune function pathways containing 1027 genes and 12 011 SNPs that characterize allergy, asthma and inflammation in a European population. We used this catalog to guide our selection of immune function SNPs from the UCSF GWA data set. Exact matching of genes was not possible because some of the genes were not assigned standard gene symbols and we therefore manually curated the lists changing 15 names to standard gene names.

Also, to standardize coverage for the genes, we identified all SNPs in the GWA data set that were within 5 kb upstream or 2 kb downstream of any of the immune function genes (17). We mapped 6629 SNPs to 911 of the 1027 genes previously defined by Loza et al. (17). If a SNP was within the vicinity of multiple genes, then it was matched to each of these genes. Supplementary Figure 1 (available at Carcinogenesis Online) shows the distribution of the number of SNPs representing each gene in the UCSF data set. UCSF investigators then sent the list of 6629 SNPs together with gene and pathway annotations to MDA investigators who identified 6029 SNPs on 893 genes belonging to 17 immune pathways using their GWA data. No gene was included in more than one pathway.

Because of the relatively small numbers of glioblastoma cases in the UK and Mayo Clinic data sets, we did not attempt full gene- and pathway-level analyses in these data sets. Instead, we used these GWA data to assess consistency of 20 SNPs that had Cochran–Armitage trend test p values (ptrend) <0.05 for association with glioblastoma in both the UCSF and the MDA data. Seventeen of these 20 SNPs were available in the Mayo Clinic GWA data set.

Statistical methods

SNP-level associations.

For each of the 6629 SNPs in the UCSF data set and the 6029 SNPs in the MDA data set, we calculated the Cochran–Armitage ptrend statistic. We next identified 20 SNPs associated with glioblastoma (ptrend ≤ 0.05) in both data sets and computed ptrend for each of the 20 SNPs in the UK data set and for the 17 available SNPs in the Mayo Clinic data set. Results from the four studies were summarized using a Bonferroni-adjusted Mantel–Haenszel statistic (adjusted for 20 hypothesis tests).

Gene-level associations: the minP approach.

To combine the multiple SNP-level association p values into one number that represents the gene-level association, and at the same time, to control for the problem that genes with more SNPs are more likely to have lower p values than genes with fewer SNPs, we used the minP approach, first proposed by Westfall and Young (21) to adjust for family wise (gene level) type I error. The outcome of this analysis is a set of minP statistics that quantify associations with glioblastoma at the gene level and are directly comparable among genes. The minP approach performs well in simulations using a variety of scenarios (25) and has several advantages. First, the empirical null distribution generated via permutation of case–control labels preserves the correlation among SNPs in the same gene, implicitly taking account of LD among SNPs. It is therefore more statistically powerful than tests, which ignore these correlations (30). Second, the derived gene-level association p values control for different numbers of SNPs among genes. A conceptually similar method was used in a paper on associations between DNA repair genes and lung cancer risk (23).

The specific procedures for the minP analysis are as follows. For each SNP, we calculated its marginal p value for association with glioblastoma risk (ptrend). For each gene, we combined the ptrend values of all SNPs mapped to that gene (as described above) by obtaining their minimum, and this minimum ptrend value became the gene-level test statistic. Suppose gene G has N SNPs and denote the marginal association ptrend value for the jth (j = 1, …, N) SNP located in the gene to be p, the test statistic for gene-level association is defined as . To evaluate the significance of minP0, we used permutations to generate its null distribution by randomly shuffling case–control status a thousand times. Let p be the p value for the jth SNP assigned to gene G in the bth permuation, the permuted minP statistic is Then the permutation-adjusted p value for the minP statistic can be calculated as the proportion of {minP,b = 1, …,B} that are smaller or equal to the observed minP statistic minP0,. The procedures described above were applied to the UCSF data and then repeated using the 6029 SNPs on 893 genes from the MDA data set. (Note that throughout the present manuscript, we use the lower case ptrend to represent the p value associated with the Cochran–Armitage trend test on a single SNP, whereas the upper case P in minP is used together with a subscript (minP0) to represent the smallest ptrend among all SNPs matched to that gene. The symbol minP with no subscript represents the permutation-adjusted minimum p value for each gene. With each data set, UCSF and MDA, we used the binomial test to evaluate the probability that more genes than expected by chance had minP ≤ 0.05 (31).

Pathway-level associations.

We applied two approaches to examine correlative patterns between the minP values from the UCSF and MDA data sets within each of the 17 pathways. These approaches address the question: within each pathway, do the MDA and UCSF data sets show that similar genes are associated with glioblastoma risk? In the first approach, for each pathway, we used the binomial test (31) to calculate the probability of observing an equal or greater number of genes that were significant (minP ≤ 0.05) in both data sets. Let n be the number of genes in a specific pathway, p1 and p2 be the percentages of significant genes (minP ≤ 0.05) in the UCSF and MDA data sets, respectively, and k be the number of genes that are significant in both datasets. The binomial p value equals the probability of having k or more significant genes among the n genes given that the background probability is p1Xp2. In the second approach, we calculated Pearson correlation coefficients (31) between −log10(minP) values from both sites. The significance of the correlation coefficients was determined by 1000 permutations that disrupted the pairings of genes between the two sites. The first approach relies on setting an (arbitrary) significance threshold to classify genes into two categories and can therefore be less powerful than the second approach, which examines the overall trend of concordance in the entire pathway. We used one-sided right-tailed tests for both approaches. Results from the two tests can be used together to determine the likelihood that genes in a pathway are similarly related to glioblastoma risk in the two data sets. If both tests are significant, it suggests that in both studies, more of the same genes are significantly associated with glioblastoma and that there is general linear correlation of minP values across the two data sets. If the correlation test p value is significant but the binomial test is not, it suggests that although there is concordant trend in p values between the two groups, there are not enough statistically significant genes in both sites for binomial p value significance.

Using pathway enrichment analyses methods similar to those used for pathway analysis in expression microarrays (32), we also conducted preliminary pathway enrichment analyses. These analyses address the question: does a particular pathway have more genes than expected by chance that are significantly associated with glioblastoma risk? For these analyses, we used Fisher's exact and Wilcoxon tests (32,33) to evaluate whether each pathway is enriched with genes that are associated with glioblastoma risk (those are the genes that had low minP values). Unlike the pathway correlative analysis discussed above, they were carried out separately for the UCSF and MDA data sets at each site. However, since the UCSF (525 cases and 602 controls) and MDA (531 cases and 1782 controls) data sets did not have sufficient numbers of observations to produce definitive pathway enrichment findings, we only include these results in the supplementary materials (see supplementary Table 1 is available at Carcinogenesis Online).

Results

SNP-level associations

There were more SNPs with nominal associations ptrend ≤ 0.05 from both UCSF and MDA sites than expected by chance (n = 20, Fisher's exact test, p = 0.04; see supplementary Table 2, available at Carcinogenesis Online, for a complete list). Four SNPs had odds ratios that were consistent across all four data sets (UCSF, MDA, Mayo and UK); these included three epidermal growth factor receptor (EGFR) SNPs (rs6969537, rs1015793 and rs11979158) and one MICB SNP with Bonferroni-adjusted Mantel–Haenszel p values of 0.0042, 3.04 × 10−5, 1.19 × 10−5 and 0.0013, respectively (Table I).

Table I.

Consistent associations of glioblastoma risk with allergy- and inflammation-related SNPs in the UCSF, MDA, UK and Mayo Clinic Studies

Gene/SNPs	Minor allele	UCSF Study (525/602)a		MDA Study (531/1782)a		UK Study (266/2482)a		Mayo Study Clinic (114/111)a		Mantel–Haenszel combined
		Case/control minor allele frequency	Additive odds ratio (95% CI)bc	Case/control minor allele frequency	Additive odds ratio (95% CI)bc	Case/control minor allele frequency	Additive odds ratio (95% CI)bc	Case/control minor allele frequency	Additive odds ratio (95% CI)bc	Mantel–Haenszel P value	Bonferroni-adjusted Mantel–Haenszel P valued
EGFR
rs6969537	A	0.12/0.15	0.80 (0.60–0.99)	0.11/0.15	0.68 (0.55–0.85)	0.14/0.15	0.89 (0.67–1.17)	0.13/0.16	0.75 (0.44–1.26)	2.0 × 10−4	4.2 × 10−3
rs1015793	G	0.11/0.15	0.70 (0.55–0.90)	0.11/0.15	0.70 (0.56–0.87)	0.12/0.15	0.78 (0.58–1.03)	0.11/0.18	0.55 (0.32–0.95)	1.5 × 10−6	3.0 × 10−5
rs11979158	G	0.12/0.17	0.70 (0.55–0.89)	0.13/0.17	0.76 (0.71–0.93)	0.13/0.17	0.68 (0.51–0.89)	0.11/0.19	0.50 (0.28–0.87)	6.0 × 10−7	1.2 × 10−5
TGF-alpha
rs3755377	C	0.48/0.43	1.26 (1.06–1.50)	0.48/0.43	1.40 (1.14–1.73)	0.42/0.44	0.87 (0.62–1.12)	0.52/0.48	1.26 (0.84–1.91)	1.2 × 10−3	2.4 × 10−2
HLA-B
rs2596503	T	0.21/0.18	1.25 (1.01–1.54)	0.23/0.20	1.21 (1.01–1.46)	0.22/0.21	1.08 (0.84–1.38)	0.23/0.22	1.15 (0.74–1.79)	9.0 × 10−4	1.7 × 10−2
MICB
rs3130922	A	0.37/0.32	1.23 (1.04–1.46)	0.39/0.33	1.49 (1.23–1.80)	0.34/0.32	1.09 (0.86–1.40)	0.35/0.32	1.11 (0.75–1.64)	6.5 × 10−5	1.3 × 10−3

CI, confidence interval.

Glioblastoma cases/controls.

Odds ratios are for 0, 1, 2 minor alleles.

Twenty SNPs were selected if ptrend (trend test) <0.05 for both UCSF and MDA data sets. The SNPs in this table are a subset of 20 that have similar odds ratios across all four studies and Bonferroni-adjusted Mantel–Haenszel p < 0.05.

Adjusted for 20 hypothesis tests.

Gene-level associations

As shown in supplementary Figure 2 (available at Carcinogenesis Online), although, as expected, the observed minimum P value (minP0) for gene–disease risk was inversely associated with the number of SNPs per gene, permutation-adjusted minimum p values (minP) were independent of number of SNPs per gene and thus this potential bias was removed from the following results.

In the gene-based analyses, 59 out of 911 and 66 out of 893 genes were significantly associated (minP ≤ 0.05) with glioblastoma risk in the UCSF and MDA datasets, respectively. There is significant enrichment of glioblastoma-associated genes compared with what is expected under the null hypothesis in both datasets (UCSF binomial p = 0.048 and MDA binomial p = 0.009). Figure 1 displays the scatter plot of the minP values for all genes from both data sets. Eight genes were significantly associated with glioblastoma risk in both data sets: IL-2RA, CCL15, ITGAM, JAK1, IFNAR1, MAPK11, ITGAD and IL-18. Note that the IL-2RA gene has relatively small minP values in both analyses (UCSF minP = 0.01 and MDA minP = 0.001).

Fig. 1.

Scatter plot of minP values of allergy- and inflammation-related genes (893) present in both UCSF Adult Glioma Study and MDA Study data sets. Each gene from each study site is assigned a minP value, which represents the results of a test of the association of that gene with glioblastoma. The minP is adjusted for multiple testing and number of SNPs per gene.

Pathway-level associations

As shown in Table II and Figure 2a and b, the cytokine signaling and the adhesion–extravasation–migration pathways have statistically significant binomial and correlation coefficient p values suggesting a similar relation to glioblastoma risk from both sites. In addition, for the cytokine signaling pathway, both sites have significantly more genes than expected that have minP values < 0.05 (14 out of 147 genes in the UCSF dataset, p = 0.007 and 13 out of 146 genes in the MDA dataset, p = 0.03; supplementary Table 1 is available at Carcinogenesis Online), further suggesting an association to glioblastoma risk for the cytokine signaling pathway from both studies.

Table II.

Correlation analyses between minPa values for gene-based associations of glioblastoma risk from UCSF Adult Glioma Study and MDA Study cases and controls within 17 immune or inflammation pathways identified by Loza et al. (17)

Pathways	Number of genes	Binomial pb	Correlation coefficient rc	Correlation coefficient pc
Natural killer cell signaling	22	1	−0.03	0.484
Cytokine signaling	146	0.029	0.15	0.003
Apoptosis signaling	65	1	−0.17	0.927
Complement caspase	37	1	−0.20	0.927
G protein-coupled receptor signaling	41	1	0.07	0.258
Glucocorticoid/PPAR signaling	19	1	0.43	0.003
ROS/glutathione/cytotoxic granules	18	1	0.03	0.410
Adhesion–extravasation–migration	119	0.003	0.30	<0.001
Eicosanoid signaling	31	1	−0.11	0.701
Innate pathogen detection	37	1	−0.03	0.549
MAPK signaling	107	0.209	0.12	0.116
Nuclear factor-kappaB signaling	31	1	−0.01	0.504
Phagocytosis-Ag presentation	32	1	−0.19	0.857
Leukocyte signaling	105	1	−0.07	0.791
Calcium signaling	14	1	−0.31	0.878
PI3K/AKT signaling	35	1	−0.03	0.528
Tumor necrosis factor superfamily signaling	34	1	0.10	0.294

ROS, reactive oxygen species.

minP is a permutation-based association p value for the association of genes with glioblastoma risk.

Binomial test p value calculates a one-sided (right-tailed) probability of observing the number of genes that are significant (minP ≤ 0.05) in both studies.

The Pearson correlation coefficient describes the association between −log10(minP) values from each pathway at both sites. Significance of the observed correlation coefficients indicates the probability of finding an equal or higher correlation coefficient from 1000 random permutations.

Fig. 2.

(a and b) Scatter plots of minP of genes by pathway from UCSF Adult Glioma Study and MDA Study data sets. Each gene from each study site is assigned a minP value, which represents the results of a test of the association of that gene with glioblastoma. The minP is adjusted for multiple testing and number of SNPs per gene.

Consistent with these findings, all but one of the eight genes having minP ≤ 0.05 in Figure 1 are found on the cytokine signaling or adhesion–extravasation–migration pathways. The exception is the MAPK11 (p38 MAP kinase) gene on the MAPK signaling pathway. In addition, the pathway enrichment analysis (supplementary Table 1 is available at Carcinogenesis Online) provides some support for a role for the cytokine signaling pathway. Using Fisher's exact test, the cytokine signaling pathway had the second and first lowest p values in the UCSF and MDA analyses, respectively.

Although the glucocorticoid/PPAR signaling pathway correlation coefficient is also statistically significant (Table II), its binomial p value is 1.0, indicating a lack of significantly low p values for gene–glioblastoma associations for this pathway.

Discussion

Our gene- and pathway-level analyses of immune function genes suggest a role for eight immune function genes and two immune function pathways: the cytokine signaling and adhesion–extravasation–migration pathways (Table II, Figure 2a and b). The cytokine signaling pathway consists of immunomodulatory proteins that act as short-range signaling molecules between cells. There is extensive evidence for a central role of cytokines in glioma growth and angiogenesis (34). This pathway includes the IL-2RA (CD25) receptor gene (Figure 1, Figure 2a), expressed on the surface of a subset of regulatory T cells (Tregs) that contribute to immunosuppression that is characteristic of the glioblastoma microenvironment (35–45).

The adhesion–extravasation–migration pathway, also associated with glioblastoma risk in both data sets, includes the CCL15 gene (Figures 1 and 2b), a member of a class of genes that regulates leukocyte trafficking across the blood–brain barrier (46). Although there is no evidence for a specific association between the CCL15 gene and glioblastoma, this chemokine may mediate Treg infiltration of melanoma (47) and also induces cell migration and differentiation of human eosinophilic leukemia (48). It is therefore possible that it may also play an as yet unidentified role in glioblastoma risk.

Limitations of our gene and pathway findings rest first on the limits of GWA data. These data do not include comprehensive coverage of all genes and pathways that may be important in understanding associations between inherited variation in immune function and glioblastoma risk. However, reanalysis of GWA is useful for illuminating potentially important areas for future, in-depth, candidate gene and pathway studies. Another potential source of error comes from misspecification of immune function pathways. Although Loza et al. (17) assigned genes to each of the 17 immune function pathways, genes appropriately included in each pathway may vary with each pathway-specific function, its physiological environment and location. In addition, although population stratification is minimized using genetic- and pathway-level analyses, it may not be eliminated. Further potential errors in interpretation may arise from the heterogeneity of glioblastoma (49) with the possibility of the immune system playing different roles in different subtypes. In the present study, to maintain adequate statistical power, we have treated this tumor as a single entity thus possibly obscuring associations that may exist between immune function genes or pathways for specific tumor subtypes. In fact, a larger sample may have allowed us to identify more genes or pathways related to glioblastoma risk. Finally, we did not evaluate gene–gene or gene–environment interactions although both undoubtedly play central (although presently unknown) roles in glioblastoma development (50,51).

Although there are many potential sources of error and bias in our study, evidence for the plausibility of our gene and pathway results comes from the extensive literature on the role of the IL-2RA (CD25) receptor in carcinogenesis. Specifically, this receptor, expressed on the surface of immunosuppressive CD4+ CD25+ Tregs, is found in proportionately higher levels in the peripheral circulation of cancer patients, including patients with glioma, than in that of controls (35–45). This excess has been attributed to the antitumor immunity-inhibiting role of these cells and has been documented in both glioma and glioblastoma (52). Learn et al. (53) compared differences in T-cell gene expression profiles in individuals with and without glioma. They found that genes in glioma patients involved in T-cell receptor ligation were downregulated, whereas genes associated with Tregs and their immunosuppressive cytokines were upregulated. Grauer et al. (54) showed that Tregs gradually accumulate in murine gliomas and suppress antitumor immunity. El Andalossui and Lesniak (55) confirmed this observation in humans noting that FOXP3-expressing Tregs increase during human glioma progression and that this increase is correlated with tumor grade. The association between Tregs and glioma, however, is complex as an messenger RNA expression study of glioblastoma tissue conducted by Schwartzbaum et al. (56) illustrates. These authors found downregulation of Treg-associated immunosuppressive cytokines (e.g. IL-10) with increased expression of CD133, an indicator of tumor progression.

With respect to single-SNP–glioma risk associations, odds ratios for three SNPs in EGFR results were consistent across the four data sets (Table I). EGFR has a direct role in inflammation and subsequent immunosuppression and both its amplification and mutation have been frequently observed in glioblastoma tumors (57). Elevated EGFR expression is also a negative prognostic indicator (58,59). In addition, three reports suggest that inherited variation in EGFR is related to glioblastoma risk (60–62) with three EGFR SNPs reported in the present study being highly correlated (0.47 ≤ r2 ≤ 1.0) with four other EGFR SNPs (rs759171, rs17172430, rs17172432 and rs17172433) found to be significantly associated with glioma risk (62). The latter paper includes UK cases and controls presented in aggregate with data from the Swedish, Danish and Finnish INTERPHONE glioma case–control studies. In the present study, we report the UK data separately. Because there is an overlap between subjects in our report and that of Andersson et al. (62) and our results do not achieve genome-wide significance, further replication of the association among these EGFR variants and glioma risk is warranted.

There are apparent differences between the SNP- and gene-level associations. Specifically, although three EGFR SNPs were statistically significant, we found only marginal significance of EGFR at the gene level (minP = 0.11 and 0.06, respectively) in UCSF and MDA data sets. This difference can be attributed, in part, to the way we estimated associations at each level of analysis. In the SNP-level analysis, we evaluated whether identical SNPs were related to glioblastoma across data sets, whereas the gene-level comparisons were based on the SNP associated with the smallest adjusted p value (minP) for each gene independently in each data set. There are at least two additional factors, besides a low ptrend that determine whether a SNP will attain gene-level significance. The first is the LD structure of SNPs within the gene and the second is the number of SNPs in the gene that are genotyped. The three EGFR SNPs (rs6969537, rs1015793 and rs11979158) that we report in Table I don't represent independent association signals; in fact, they are linked with each other, with their LD r2 ranging from 0.58 to 0.74. In addition, EGFR is a fairly long gene, spanning >130 Kb, and has 43 SNPs on the UCSF genotype panel, whereas an average immune gene that we queried has only 7 SNPs. The minP approach takes into account the correlation structure among SNPs so that closely linked signals are not erroneously treated as independent and given additional weights. Furthermore, the multiple testing procedures prevent long genes from attaining higher significance purely on the basis of their larger number of SNPs. As a result of the two factors, the EGFR gene did not obtain significance at the gene level. This example aptly demonstrates the advantages of the minP approach and the stringency that we applied to our study to avoid reporting false-positive findings. Nonetheless, the marginal significance that we observed at the SNP level may have given rise to statistical significance had we been able to combine the two data sets (see Methods ) at the gene level. In any case, further independent replication of the association between these and linked SNPs in the EGFR gene and glioblastoma is now warranted.

Results for IL-2RA are also different for SNP- and gene-level analyses. The IL-2RA gene is represented by rs4749926 in UCSF data and rs791589 in MDA data. Neither one of these SNPs is significant in both data sets thus explaining, in part, the absence of SNPs from this gene in Table I.

Results from previous studies of asthma- and allergy-related SNPs and haplotypes have been inconsistent (1,6–9). None of the 13 allergy- or asthma-related SNPs (11 in IL-4Ralpha and 2 in IL-13 SNPs) in the UCSF or the MDA data set were statistically significantly (ptrend < 0.05) associated with glioblastoma. Nonetheless, our findings suggest that this relation may result from the relative absence of immunosuppression among people with allergies. Specifically, allergies are associated with deficiency in CD4+ CD25+ Treg function (63,64), whereas, as previously noted, glioma patients have a relatively high proportion of CD4+ CD25+ Tregs in their peripheral blood (45).

Previous work has also suggested that the innate immune system may be related to glioblastoma risk. Rajaraman et al. (11) reported results of a study of 1397 SNPs related to innate immunity and both glioma and glioblastoma risk. Although both the UCSF and the MDA data sets included five of the eight genes on which these investigators found SNPs significantly associated with glioblastoma in their data set, none of these genes were statistically significantly related to glioblastoma.

To summarize, we note different limitations and advantages of both genome-wide and candidate gene approaches. Our analyses indicate the importance of inherited variation in the cytokine signaling and adhesion–extravasation–migration pathways (most notably the IL-2RA gene on the cytokine pathway) in glioblastoma risk. They also possibly suggest that inherited variation in EGFR may influence risk, supporting results recently published by Andersson et al. (62). IL-2RA and other genes on the cytokine signaling pathway may contribute to an immunosuppressive microenvironment necessary for tumor growth (2), whereas EGFR can act as an inflammation signaling hub. Since GWA data are not designed to examine specific candidate genes, our results suggest that more comprehensive candidate gene and pathway analyses are now warranted for studies for IL-2RA and the cytokine signaling pathway.

Supplementary material

Supplemental Figures 1 and 2 and Tables 1 and 3 can be found at http://carcin.oxfordjournals.org/.