Abstract
Genome-wide association studies (GWAS) have so far identified 25 loci associated with glioma risk, with most showing specificity for either glioblastoma (GBM) or non-GBM tumors. The majority of these GWAS susceptibility variants reside in noncoding regions and the causal genes underlying the associations are largely unknown. Here we performed a transcriptome-wide association study to search for novel risk loci and candidate causal genes at known GWAS loci using Genotype-Tissue Expression Project (GTEx) data to predict cis-predicted gene expression in relation to GBM and non-GBM risk in conjunction with GWAS summary statistics on 12,488 glioma cases (6,183 GBM and 5,820 non-GBM) and 18,169 controls. Imposing a Bonferroni-corrected significance level of P < 5.69 × 10−6, we identified 31 genes, including GALNT6 at 12q13.33, as a candidate novel risk locus for GBM (mean Z = 4.43; P = 5.68 × 10−6). GALNT6 resides at least 55 Mb away from any previously identified glioma risk variant, while all other 30 significantly associated genes were located within 1 Mb of known GWAS-identified loci and were not significant after conditioning on the known GWAS-identified variants. These data identify a novel locus (GALNT6 at 12q13.33) and 30 genes at 12 known glioma risk loci associated with glioma risk, providing further insights into glioma tumorigenesis.
This study identifies new genes associated with glioma risk, increasing understanding of how these tumors develop.
Introduction
Diffuse gliomas are the most common malignant primary brain tumor affecting adults (1). Gliomas can be broadly classified into glioblastoma (GBM) and low-grade non-GBM tumors. Gliomas typically have a poor prognosis irrespective of medical care, with the most common form, glioblastoma multiforme (GBM), having a median overall survival of only 10–15 months (1). While the glioma subtypes have distinct molecular profiles resulting from different etiologic pathways, no environmental exposures have consistently been linked to risk except for ionizing radiation, which only accounts for a very small number of cases (1). Inherited genetic factors do, however, play an important role in the etiology of glioma and genome-wide association studies (GWAS) have so far identified common variants at 25 loci influencing disease risk (2). Perhaps not surprisingly given differences in the molecular profile of GBM and non-GBM tumors, subtype-specific associations are shown for a number of the risk variants (3, 4). Collectively, the known risk loci only account for around a third of the familial risk of both GBM and non-GBM glioma (2), indicating that additional susceptibility variants remain to be identified.
Many of the GWAS risk variants are likely to have a small effect size, and thus are difficult to identify in individual SNP-based GWAS, even with large sample numbers (2). Applying gene-based approaches that aggregate the effects of multiple variants into a single testing unit is thus attractive and offers the prospect of increasing study power. Most GWAS risk variants reside in noncoding regions and are primarily located in active chromatin regions, which are highly enriched with expression quantitative trait loci (eQTL; ref. 5). Hence transcriptome-wide association studies (TWAS) that systematically investigate the association of genetically predicted gene expression with disease risk offers a potentially attractive strategy to identify novel susceptibility genes for glioma (6, 7).
Herein, we report results from a TWAS of glioma implementing the MetaXcan (8) methodology to analyze summary statistics data from 12,488 cases and 18,169 controls of European descent. We identify 31 genes at 13 loci associated with glioma risk, and provide additional evidence of a potential role for a number of genes that are dysregulated in glioma tumorigenesis.
Materials and Methods
Ethics
A TWAS was undertaken using previously reported GWAS data (2). Ethical approval was not sought for this specific project because all data came from the summary statistics from the published GWAS, and no individual-level data were used.
GWAS data
Glioma genotyping data were derived from the most recent meta-analysis of GWAS in glioma, which related >6 million genetic variants (after imputation) to glioma, in 12,488 patients and 18,169 controls from eight independent studies of individuals of European descent (Supplementary Figs. S1 and S2; Supplementary Table S1). Comprehensive details of the genotyping and quality control of these GWAS have been reported previously (2). Gliomas are heterogeneous and different tumor subtypes, defined in part by malignancy grade [e.g., pilocytic astrocytoma World Health Organization (WHO) grade I, diffuse “low-grade” glioma WHO grade II, anaplastic glioma WHO grade III, and GBM WHO grade IV] can be distinguished. For the sake of brevity we considered gliomas as being either GBM or non-GBM tumors.
Association analysis of predicted gene expression with glioma risk
Associations between predicted gene expression and glioma risk were examined using MetaXcan (8), which combines GWAS and eQTL data, accounting for LD-confounded associations. Briefly, genes likely to be disease causing were prioritized using S-PrediXcan that uses GWAS summary statistics and prespecified weights to predict gene expression, given covariances of SNPs. SNP weights and their respective covariance for 13 brain tissues (amygdala, anterior cingulate cortex, caudate basal ganglia, cerebellar hemisphere, cerebellum, cortex, frontal cortex, hippocampus, hypothalamus, nucleus accumbens basal ganglia, putamen basal ganglia, spinal cord, and substantia nigra) from 80 to 154 individuals were obtained from predict.db (http://predictdb.org/; ref. 8), which is based on GTEx version 7 eQTL data (Supplementary Table S2). To combine S-PrediXcan data across the different brain tissues taking into account tissue–tissue correlations, we used S-MultiXcan.
To determine whether associations between genetically predicted gene expression and glioma risk were influenced by variants previously identified by GWAS, we performed conditional analyses adjusting for sentinel GWAS risk SNPs (Supplementary Table S3) using GCTA-COJO (9, 10). Adjusted output files were provided as the input GWAS summary statistics for S-PrediXcan analyses as above.
For all significant genes identified by S-MultiXcan analyses, we additionally considered the effect of the top eSNP on glioma risk. For each identified gene, the most significant eSNP for each brain tissue was identified from GTEx v7 “allpairs.txt.gz” files. Glioma GWAS summary statistics for the surrounding region were estimated after conditioning on identified significant eSNP/s using GCTA-COJO (9, 10), using “—cojo-slct” and “—cojo-p 0.05″ to select independent eSNPs and avoid collinearity in association testing.
To account for multiple comparisons, we first considered a simple Bonferroni-corrected P value threshold of 3.45 × 10−6 (i.e., 0.05/14,486 genes) to determine a statistically significant association. This is, however, inherently conservative because expression of genes can be correlated. To identify highly correlated genes, we performed a weighted correlation network analysis (WGCNA) using WGCNA v1.63 (11). Plots of soft threshold against the scale-free topology model fit were used to determine the threshold preserving 90% of topology (Supplementary Table S4). Dendrograms and heatmaps were generated to visualize coexpression of genes (Supplementary Figs. S3–S5). The number of clusters reflects the number of independent gene sets. We examined the comparability of gene clustering across brain tissues by dendrogram Z-values; with a Z-value of 5–10 corresponding to moderate preservation and a Z-value > 10 being indicative of strong preservation (Supplementary Table S5). To estimate the number of independently expressed genes per brain tissue, we assessed gene–gene adjacency (i.e., correlation) values. Significantly correlated gene–gene pairs were identified as those with adjacency values greater than three SDs from the mean. Removing at random one correlated gene from each pair left an estimate of the number of “independent genes” (Supplementary Table S4). The median number of independent genes was 8,781, which defined the TWAS Bonferroni-correct threshold as P < 5.69 × 10−6.
S-PrediXcan analyses were additionally carried out on 922 whole-blood samples from Depression Genes and Networks (DGN), to compare associations at genes identified as significant from S-MultiXcan analyses in brain, and aid interpretation of potential tissue-specific and generic eQTL effects.
Identified genes were annotated by their potential presence in the v87 COSMIC cancer gene census (https://cancer.sanger.ac.uk/cosmic/), as well as their potential overlap with copy number gains and losses as annotated in CosmicCompleteCNA.tsv.gz.
Statistical power for association tests
To assess the power of our TWAS to identify associations, we performed a simulation analysis adopting a similar strategy to Wu and colleagues (6). We set the number of cases and controls as 12,488 (6,183 GBM, 5,820 non-GBM) and 18,169, respectively. Glioma prevalence estimates were obtained from CBTRUS 2017 (12), assuming an overall incidence of primary brain and central nervous system tumors to be 22.6 per 100,000, of which 27% are gliomas and 56% of gliomas are GBM. We generated the gene expression levels from the empirical distribution of gene expression levels in GTEx normalized expression dataset for each brain tissue. We calculated statistical power at P < 5.69 × 10−6, corresponding to the TWAS genome-wide significance level, according to various cis-heritability (h2) thresholds that are assumed to be equivalent to gene expression prediction models (R2). The results, based on 1,000 replicates, are summarized in Supplementary Fig. S6.
Availability of data and materials
Genotype data from the Glioma International Case-Control Consortium Study GWAS are available from the database of Genotypes and Phenotypes (dbGaP) under accession number phs001319.v1.p1. In addition, genotypes from the GliomaScan GWAS can be accessed through dbGaP (accession number phs000652.v1.p1). Summary statistics from the glioma GWAS meta-analysis are available from the European Genome-phenome Archive (EGA, http://www.ebi.ac.uk/ega/) under accession number EGAS00001003372.
Results
We evaluated the association between predicted gene expression levels and glioma risk using MetaXcan with summary statistics for GWAS SNPs in 12,488 glioma cases and 18,169 controls. In view of associations for glioma being strongly subtype-specific (2), we analyzed TWAS results for GBM and non-GBM cases. Figure 1 shows Manhattan plots for respective TWAS associations. Quantile–quantile plots of TWAS association statistics did not show evidence of systematic inflation (Supplementary Fig. S1).
Manhattan plots of gene genomic coordinates against –log10 (P) of TWAS results. A, GBM glioma. B, Non-GBM glioma. The line represents the Bonferroni-corrected threshold of P ≤ 5.69 × 10−6.
Manhattan plots of gene genomic coordinates against –log10 (P) of TWAS results. A, GBM glioma. B, Non-GBM glioma. The line represents the Bonferroni-corrected threshold of P ≤ 5.69 × 10−6.
In total, the expression levels of 14,485 genes were tested for an association with glioma. To establish the threshold for assigning genome-wide statistical significance taking into account correlations between gene expression, we carried out WGCNA (11) analysis to determine the number of independent gene sets (Supplementary Table S4). On the basis of an estimated number of uncorrelated genes of 8,781 we imposed a Bonferroni multiple-testing threshold of P < 5.69 × 10−6 to declare significant associations.
Applying this threshold, we identified 23 genes associated with GBM (Fig. 1A), and 8 with non-GBM glioma (Fig. 1B; Table 1; Supplementary Table S6 and S7). All identified genes but 1 were within 1 Mb of previously reported glioma risk SNPs. After conditioning on the nearby GWAS glioma-risk SNP in each case gene associations were severely abrogated, consistent with the TWAS associations reflecting the previously identified GWAS associations. The exception was GALNT6 at 12q13.13, which did not map within 1 Mb of a previously identified GWAS risk SNP, and was significantly associated with GBM. The risk allele (T) of sentinel SNP rs3782473 at 12q13.13 had an association P of 9.08 × 10−8 (OR 1.15; 95% confidence interval 1.09–1.21) with GBM (Fig. 2). After conditioning on rs3782473, there were no significant TWAS associations at 12q13.13, consistent with the association signal defined by rs3782473 underlying the association with GALNT6 (Supplementary Table S6). In 9 of 13 brain regions, there was a significant association between the risk allele (T) of rs3782473 and increased expression of GALNT6 (Supplementary Fig. S7).
Genes significantly associated with risk of GBM and non-GBM glioma
| Locus | Gene | P | N/Nindep | Z-score min/max | Z-score mean | Z-score SD | Within 1 Mb of glioma risk SNP | SNP/s adjusting for | P after SNP adjustment |
|---|---|---|---|---|---|---|---|---|---|
| GBM | |||||||||
| 20q13.33 | ZGPAT | 6.85 × 10−45 | 3/3 | −0.07/14.3 | 6.69 | 7.21 | YES | rs2297440 | 8.39 × 10−3 |
| 20q13.33 | SLC2A4RG | 4.90 × 10−39 | 1/1 | 13.1/13.1 | 13.1 | — | YES | rs2297440 | 0.09 |
| 20q13.33 | ARFRP1 | 1.93 × 10−30 | 3/3 | 8.63/11.5 | 10.5 | 1.66 | YES | rs2297440 | 0.77 |
| 20q13.33 | STMN3 | 4.54 × 10−27 | 4/4 | −10.9/−0.88 | −7.70 | 4.60 | YES | rs2297440 | 0.62 |
| 5p15.33 | TERT | 5.63 × 10−26 | 2/2 | −0.43/10.7 | 5.12 | 7.86 | YES | rs10069690 | 0.63 |
| 20q13.33 | GMEB2 | 3.05 × 10−16 | 2/2 | −8.26/−8.16 | −8.21 | 0.07 | YES | rs2297440 | 0.55 |
| 5p15.33 | NKD2 | 9.46 × 10−16 | 6/4 | −0.08/4.85 | 1.49 | 1.87 | YES | rs10069690 | 1.36 × 10−4 |
| 20q13.33 | LIME1 | 3.60 × 10−13 | 2/2 | −6.64/5.24 | −0.70 | 8.40 | YES | rs2297440 | 5.11 × 10−3 |
| 16q12.1 | HEATR3 | 3.48 × 10−10 | 13/1 | 4.86/6.73 | 6.03 | 0.52 | YES | rs10852606 | 0.82 |
| 22q13.1 | BAIAP2L2 | 8.61 × 10−9 | 1/1 | 5.76/5.76 | 5.76 | — | YES | rs2235573 | 0.27 |
| 7p11.2 | EGFR | 1.35 × 10−8 | 2/2 | −4.70/−4.44 | −4.57 | 0.18 | YES | rs723527,rs75061358 | 0.46 |
| 9p21.3 | CDKN2B | 3.11 × 10−8 | 1/1 | −5.53/−5.53 | −5.53 | — | YES | rs634537 | 0.38 |
| 22q13.1 | SLC16A8 | 4.88 × 10−8 | 3/3 | 5.45/5.54 | 5.48 | 0.05 | YES | rs2235573 | 0.84 |
| 20q13.33 | HAR1A | 2.33 × 10−7 | 11/5 | −1.48/4.08 | 0.32 | 1.48 | YES | rs2297440 | 0.90 |
| 20q13.33 | OPRL1 | 6.97 × 10−7 | 2/2 | −3.99/−2.02 | −3.00 | 1.39 | YES | rs2297440 | 0.03 |
| 1p31.3 | JAK1 | 9.29 × 10−7 | 4/3 | 4.11/5.36 | 4.87 | 0.56 | YES | rs12752552 | 0.16 |
| 20q13.33 | PCMTD2 | 1.07 × 10−6 | 5/5 | −1.85/3.17 | 0.91 | 2.34 | YES | rs2297440 | 0.02 |
| 22q13.1 | CTA-228A9.3 | 1.38 × 10−6 | 4/4 | 1.24/5.04 | 3.80 | 1.76 | YES | rs2235573 | 0.44 |
| 22q13.1 | PICK1 | 1.90 × 10−6 | 7/5 | 3.12/5.78 | 4.77 | 1.00 | YES | rs2235573 | 0.38 |
| 20q13.33 | DIDO1 | 2.11 × 10−6 | 4/3 | −2.10/3.16 | 0.38 | 2.16 | YES | rs2297440 | 0.92 |
| 5.08 × 10−6 | −5.25/−1.75 | −3.32 | 1.61 | YES | rs2562152 (GBM) | 0.36 | |||
| 16p13.3a | IL9RP3 | 4/4 | rs3751667 (non-GBM) | 9.42 × 10−6 | |||||
| 20q13.33 | TCEA2 | 5.45 × 10−6 | 3/3 | 1.68/5.10 | 3.55 | 1.73 | YES | rs2297440 | 0.42 |
| 12q13.13 | GALNT6 | 5.68 × 10−6 | 10/3 | 3.10/5.26 | 4.43 | 0.68 | NO | rs3782473 | 0.82 |
| Non-GBM | |||||||||
| 11q23.3 | PHLDB1 | 4.08 × 10−32 | 2/2 | 0.02/12.0 | 6.02 | 8.49 | YES | rs12803321 | 3.71 × 10−4 |
| 11q23.3 | TREH | 1.90 × 10−16 | 1/1 | 8.23/8.23 | 8.23 | — | YES | rs12803321 | 1.20 × 10−4 |
| 20q13.33 | ZGPAT | 1.18 × 10−11 | 3/3 | 0.04/6.20 | 3.90 | 3.36 | YES | rs2297440 | 2.17 × 10−4 |
| 20q13.33 | SLC2A4RG | 4.44 × 10−11 | 1/1 | 6.59/6.59 | 6.59 | — | YES | rs2297440 | 0.09 |
| 11q23.3 | RPL5P30 | 2.09 × 10−9 | 2/2 | −6.32/−3.99 | −5.16 | 1.65 | YES | rs12803321 | 0.55 |
| 5p15.33 | TERT | 5.09 × 10−7 | 2/2 | 0.50/5.38 | 2.94 | 3.45 | YES | rs10069690 | 0.96 |
| 20q13.33 | LIME1 | 3.78 × 10−6 | 2/2 | −3.66/4.24 | 0.29 | 5.58 | YES | rs2297440 | 3.32 × 10−3 |
| 11q23.3 | TMEM25 | 5.15 × 10−6 | 2/2 | 4.74/4.91 | 4.83 | 0.12 | YES | rs12803321 | 0.23 |
| Locus | Gene | P | N/Nindep | Z-score min/max | Z-score mean | Z-score SD | Within 1 Mb of glioma risk SNP | SNP/s adjusting for | P after SNP adjustment |
|---|---|---|---|---|---|---|---|---|---|
| GBM | |||||||||
| 20q13.33 | ZGPAT | 6.85 × 10−45 | 3/3 | −0.07/14.3 | 6.69 | 7.21 | YES | rs2297440 | 8.39 × 10−3 |
| 20q13.33 | SLC2A4RG | 4.90 × 10−39 | 1/1 | 13.1/13.1 | 13.1 | — | YES | rs2297440 | 0.09 |
| 20q13.33 | ARFRP1 | 1.93 × 10−30 | 3/3 | 8.63/11.5 | 10.5 | 1.66 | YES | rs2297440 | 0.77 |
| 20q13.33 | STMN3 | 4.54 × 10−27 | 4/4 | −10.9/−0.88 | −7.70 | 4.60 | YES | rs2297440 | 0.62 |
| 5p15.33 | TERT | 5.63 × 10−26 | 2/2 | −0.43/10.7 | 5.12 | 7.86 | YES | rs10069690 | 0.63 |
| 20q13.33 | GMEB2 | 3.05 × 10−16 | 2/2 | −8.26/−8.16 | −8.21 | 0.07 | YES | rs2297440 | 0.55 |
| 5p15.33 | NKD2 | 9.46 × 10−16 | 6/4 | −0.08/4.85 | 1.49 | 1.87 | YES | rs10069690 | 1.36 × 10−4 |
| 20q13.33 | LIME1 | 3.60 × 10−13 | 2/2 | −6.64/5.24 | −0.70 | 8.40 | YES | rs2297440 | 5.11 × 10−3 |
| 16q12.1 | HEATR3 | 3.48 × 10−10 | 13/1 | 4.86/6.73 | 6.03 | 0.52 | YES | rs10852606 | 0.82 |
| 22q13.1 | BAIAP2L2 | 8.61 × 10−9 | 1/1 | 5.76/5.76 | 5.76 | — | YES | rs2235573 | 0.27 |
| 7p11.2 | EGFR | 1.35 × 10−8 | 2/2 | −4.70/−4.44 | −4.57 | 0.18 | YES | rs723527,rs75061358 | 0.46 |
| 9p21.3 | CDKN2B | 3.11 × 10−8 | 1/1 | −5.53/−5.53 | −5.53 | — | YES | rs634537 | 0.38 |
| 22q13.1 | SLC16A8 | 4.88 × 10−8 | 3/3 | 5.45/5.54 | 5.48 | 0.05 | YES | rs2235573 | 0.84 |
| 20q13.33 | HAR1A | 2.33 × 10−7 | 11/5 | −1.48/4.08 | 0.32 | 1.48 | YES | rs2297440 | 0.90 |
| 20q13.33 | OPRL1 | 6.97 × 10−7 | 2/2 | −3.99/−2.02 | −3.00 | 1.39 | YES | rs2297440 | 0.03 |
| 1p31.3 | JAK1 | 9.29 × 10−7 | 4/3 | 4.11/5.36 | 4.87 | 0.56 | YES | rs12752552 | 0.16 |
| 20q13.33 | PCMTD2 | 1.07 × 10−6 | 5/5 | −1.85/3.17 | 0.91 | 2.34 | YES | rs2297440 | 0.02 |
| 22q13.1 | CTA-228A9.3 | 1.38 × 10−6 | 4/4 | 1.24/5.04 | 3.80 | 1.76 | YES | rs2235573 | 0.44 |
| 22q13.1 | PICK1 | 1.90 × 10−6 | 7/5 | 3.12/5.78 | 4.77 | 1.00 | YES | rs2235573 | 0.38 |
| 20q13.33 | DIDO1 | 2.11 × 10−6 | 4/3 | −2.10/3.16 | 0.38 | 2.16 | YES | rs2297440 | 0.92 |
| 5.08 × 10−6 | −5.25/−1.75 | −3.32 | 1.61 | YES | rs2562152 (GBM) | 0.36 | |||
| 16p13.3a | IL9RP3 | 4/4 | rs3751667 (non-GBM) | 9.42 × 10−6 | |||||
| 20q13.33 | TCEA2 | 5.45 × 10−6 | 3/3 | 1.68/5.10 | 3.55 | 1.73 | YES | rs2297440 | 0.42 |
| 12q13.13 | GALNT6 | 5.68 × 10−6 | 10/3 | 3.10/5.26 | 4.43 | 0.68 | NO | rs3782473 | 0.82 |
| Non-GBM | |||||||||
| 11q23.3 | PHLDB1 | 4.08 × 10−32 | 2/2 | 0.02/12.0 | 6.02 | 8.49 | YES | rs12803321 | 3.71 × 10−4 |
| 11q23.3 | TREH | 1.90 × 10−16 | 1/1 | 8.23/8.23 | 8.23 | — | YES | rs12803321 | 1.20 × 10−4 |
| 20q13.33 | ZGPAT | 1.18 × 10−11 | 3/3 | 0.04/6.20 | 3.90 | 3.36 | YES | rs2297440 | 2.17 × 10−4 |
| 20q13.33 | SLC2A4RG | 4.44 × 10−11 | 1/1 | 6.59/6.59 | 6.59 | — | YES | rs2297440 | 0.09 |
| 11q23.3 | RPL5P30 | 2.09 × 10−9 | 2/2 | −6.32/−3.99 | −5.16 | 1.65 | YES | rs12803321 | 0.55 |
| 5p15.33 | TERT | 5.09 × 10−7 | 2/2 | 0.50/5.38 | 2.94 | 3.45 | YES | rs10069690 | 0.96 |
| 20q13.33 | LIME1 | 3.78 × 10−6 | 2/2 | −3.66/4.24 | 0.29 | 5.58 | YES | rs2297440 | 3.32 × 10−3 |
| 11q23.3 | TMEM25 | 5.15 × 10−6 | 2/2 | 4.74/4.91 | 4.83 | 0.12 | YES | rs12803321 | 0.23 |
NOTE: Detailed are the S-MultiXcan P for association between gene expression and GBM/non-GBM risk and the corresponding Z-scores quantifying this relationship (e.g., a positive score indicates increased gene expression increases risk of GBM or non-GBM glioma). N and Nindep indicate the total number of single-tissue results used for S-MultiXcan analysis and the number of independent components after singular value decomposition, respectively.
aSpecific GBM and non-GBM signals have been reported at 16p13.33 (2).
Regional plot of association results, recombination rates, and chromatin state segmentation tracks at 12q13.33 in GBM glioma. Plot shows discovery association results of both genotyped (triangles) and imputed (circles) SNPs in the GWAS samples and recombination rates. −Log10 P values (y-axes) of the SNPs are shown according to their chromosomal positions (x-axes). The lead SNP rs3782473 is shown as a large circle. The color intensity of each symbol reflects the extent of LD with the top genotyped SNP, white (r2 = 0) through to dark red (r2 = 1.0). Genetic recombination rates, estimated using HapMap samples from Utah residents of western and northern European ancestry (CEU) are shown with a light blue line. Physical positions are based on NCBI build 37 of the human genome. Also shown are the relative positions of GENCODE v19 genes mapping to the region of association. Below the association plot, the location of GALNT6 eSNPs is indicated, as well as the relative positions of GENCODE v19 genes mapping to the region of association and the chromatin state segmentation tracks (ChromHMM) for H1 and H9 neural progenitor cells derived from the epigenome roadmap project, as per the legend. TSS, transcriptional start sites.
Regional plot of association results, recombination rates, and chromatin state segmentation tracks at 12q13.33 in GBM glioma. Plot shows discovery association results of both genotyped (triangles) and imputed (circles) SNPs in the GWAS samples and recombination rates. −Log10 P values (y-axes) of the SNPs are shown according to their chromosomal positions (x-axes). The lead SNP rs3782473 is shown as a large circle. The color intensity of each symbol reflects the extent of LD with the top genotyped SNP, white (r2 = 0) through to dark red (r2 = 1.0). Genetic recombination rates, estimated using HapMap samples from Utah residents of western and northern European ancestry (CEU) are shown with a light blue line. Physical positions are based on NCBI build 37 of the human genome. Also shown are the relative positions of GENCODE v19 genes mapping to the region of association. Below the association plot, the location of GALNT6 eSNPs is indicated, as well as the relative positions of GENCODE v19 genes mapping to the region of association and the chromatin state segmentation tracks (ChromHMM) for H1 and H9 neural progenitor cells derived from the epigenome roadmap project, as per the legend. TSS, transcriptional start sites.
For many loci, our TWAS findings broadly support the involvement of a number of genes that have previously been proposed to be implicated in defining glioma risk (3). Specifically, single-gene associations were identified at 1p31.3 (JAK1), 7p11.2 (EGFR), 9p21.3 (CDKN2B), and 16q12.1 (HEATR3). However, at a number of loci our analysis identified multiple significant genes, notably 5p15.33 (TERT and NKD2), 11q23.3 (PHLBD1, TREH, RPL5P30, and TMEM25), and 20q13.33 (ZGPAT, SLC2A4RG, ARFRP1, STMN3, GMEB2, LIME1, HAR1A, OPRL1, PFMTD2, DIDO1, and TCEA2). No significant genes were identified at nine previously reported glioma risk loci (3p14.1, 8q24.21, 10q24.33, 10q25.2, 11q23.2, 12q12.1, 14q12, 15q24.1, and 17p13.1).
To explore the possibility of generic eQTL effects, we considered S-PrediXcan analyses at the 31 identified genes using 922 whole-blood samples from the DGN study (Supplementary Table S8). Twelve genes were significantly associated at P < 0.05 and had a consistent direction of effect with S-MultiXcan analyses (GBM: JAK1 at 1p31.3, TERT at 5p15.33, GALNT6 at 12q13.13, HEATR3 at 16q12.1, ZGPAT, ARFRP1, GMEB2, LIME1, and PCMTD2 at 20q13.33; non-GBM: TERT at 5p15.33, TMEM25 at 11q23.3, and ZGPAT at 20q13.33), 6 genes were inconsistent (GBM: CDKN2B at 9p21.3, SLC2A4RG, STMN3, and OPRL1 and TCEA2 at 20q13.33; non-GBM: PHLDB1 at 11q23.3), 4 genes were not significantly associated (GBM: DIDO1 at 20q13.33 and BAIAP2L2, and PICK1 at 22q13.1; non-GBM: SLC2A4RG at 20q13.33), and 8 genes could not be assessed (GBM: NKD2 at 5p15.33, EGFR at 7p11.2, IL9RP3 at 16p13.3, HAR1A at 20q13.33 and SLC16A8, and CTA-228A9.3 at 22q13.1; non-GBM: TREH and RPL5P30 at 11q23.3).
Following on, we further investigated the relationship between the 31 genes identified as significantly associated with GBM or non-GBM by examining associations after adjusting for the top eSNP/s at each gene (Supplementary Fig. S8). For most loci, association signals were abrogated after adjusting for the top eSNP/s, consistent with variation in expression of the identified gene being functional. In contrast, the association signals at 11q23.3 and 20q13.33 were only really affected by adjusting for multiple rather than individual gene eSNPs, raising the possibility of combinatorial effects. Intriguingly at 7p11.2, which is characterized by two independent risk loci (marked by rs75061358 and rs723527, respectively), after adjustment for the EGFR eSNPs the rs75061358 signal disappears, while the rs723527 signal is unaffected, perhaps indicative of an additional distinct as yet unidentified functional mechanism.
Finally, we compared overlap of the 31 identified genes with presence in the COSMIC cancer gene census as an oncogene or tumor suppressor gene, as well as whether the given gene is subject to copy number gains and/or losses (Supplementary Table S9). Most TWAS directions of effect are consistent with the gene's probable role in tumorigenesis, such as the tumor suppressor gene CDKN2B, whereby decreased expression is associated with increased glioma risk. However, at 7p11.2 increased expression of the oncogene EGFR, which is commonly upregulated in gliomas, was found by S-MultiXcan analyses to be negatively associated with glioma risk, perhaps indicative of different mechanisms before and after tumor initiation.
Discussion
In this large TWAS involving 12,488 glioma cases of European ancestry, we identified genetically predicted expression levels in 23 genes associated with GBM, and 8 with non-GBM glioma risk. One of these genes, GALNT6, is located at least 55 Mb away from any previously identified GWAS glioma variant, consistent with it representing a potential novel risk locus. All other 30 genes identified were located within 1 Mb of known GWAS loci, including 14 genes at three loci that had not previously been associated with glioma risk.
Our findings provide further support study for a number of the genes previously implicated by GWAS whose expression influences the risk of developing glioma. These include JAK1 at 1p31.3, PHLDB1 at 11q23.3, EGFR at 7p11.2, and HEATR3 at 16q12.1. In addition, our TWAS implicates new genes at known glioma loci, including TMEM25 at 11q23.3 and NKD2 at 5p15.33 as playing a role in defining risk of non-GBM and GBM tumors, respectively. TMEM25 has been identified as a member of the immunoglobulin superfamily, whose members are implicated in immune responses, growth factor signaling, and cell adhesion (13). Intriguingly, NKD2 encodes a Wnt pathway inhibitor that is hypermethylated in a large proportion of GBM tumors (14). The functional consequence of rs10069690 at 5p15.33 has previously been reported to be due to the risk allele (A) creating an additional splice donor site in the fourth intron of TERT, resulting in expression of a dominant negative transcript inhibiting telomerase (15). Therefore, the TWAS association with TERT may not be directly due to cis-regulatory effects but as an indirect consequence of this dominant negative effect, with a possible, albeit currently undetermined, effect on expression of NKD2.
In addition to refining the genes underscoring previously reported GWAS associations, our TWAS study identified a new gene, GALNT6 at 12q13.33, a locus not previously identified as playing a role in GBM. The gene product of GALNT6 is polypeptide N-acetylgalactosaminyltransferase 6, which is a class of proteins frequently disrupted in cancers (16). Of note is that GALNT6 expression regulates EGFR activity (17). While requiring further investigation, GALNT6 and rs3782473 represent a promising new glioma risk locus.
A large number of genes associated with glioma risk were located at 20q13.33. These include DIDO1, PCMTD2, HAR1A, and TCEA2. HAR1A expression is reduced in GBM and has been shown to be a prognostic biomarker for diffuse glioma (18). While DIDO1, PCMTD2, and TCEA2 have not previously been shown to be associated with glioma, DIDO1 promotes cell fate differentiation in embryonic stem cells (19) and TCEA2 encodes transcription elongation factor A protein 2, which interacts with BRCA1 (20). Future work will be required to reveal the contribution of these genes to glioma development and determine whether any are acting as “passengers.”
A number of previously reported glioma risk loci were not implicated in our TWAS. The reason may be obvious for some loci where the demonstrated functional mechanism is not mediated through a cis-regulatory effect on gene expression and therefore is unlikely to be detected by TWAS (e.g., at 17p13.1 the SNP rs78378222 directly affects TP53 mRNA poly-adenylation; ref. 21). At other loci such as 8q24.21, it is less obvious why an association was not detected. It may be that adult brain tissues do not represent the best model for these loci, as many genes in this region were not retained for the TWAS [genes were only retained if the nested cross-validated correlation between predicted and actual levels >0.10 (R2 > 1%) and P of the correlation test < 0.05]. Indeed, we observed a far larger number of significant genes for GBM than non-GBM loci. Speculatively, models at earlier developmental stages may yield greater insights at these loci, especially if they are influencing differentiation down oligodendrocyte/astrocyte lineages. In addition, other mechanistic effects may explain the functional basis of such loci, including methylation and splicing.
Our ability to identify genes significantly associated with glioma risk in this TWAS has inevitably been affected by tissue specificity and the sample size of the dataset used in the genetic prediction model of gene expression. Because of the importance of tissue- or cell-specific regulators in governing development and function, we have sought to analyze the most appropriate tissue-specific model to best capture the transcriptional regulatory mechanisms relevant to deciphering glioma development. Here we have sought to analyze an appropriate tissue transcriptome to model gene expression. We acknowledge that brain tissue does, however, comprise both neurons and glial cells (which include oligodendrocytes, astrocytes, ependymal cells, Schwann cells, microglia, and satellite cells). However, in light of abundant shared cis-regulation of expression across multiple brain tissues (22), by combining data on multiple brain tissues, we would expect any model to yield greater power as the number of tissues in which a variant is functional increases. Hence we aimed to robustly capture genetically regulated genes expression using a large sample size.
In conclusion, this study identified new genes whose predicted expression is associated with glioma and serves to illustrate that the TWAS approach can be a useful method of utilizing preexisting GWAS to identify new susceptibility genes. On the basis of the power calculation, our TWAS analysis had only 80% power to detect an OR of around 1.1 or 1.2 for GBM or non-GBM glioma risk per one SD increase (or decrease) in the expression level of a gene whose cis-heritability is 60% and 20%, respectively. Hence, the application of TWAS based on larger eQTL and GWAS datasets is likely to provide further insights into the genetics of glioma.
Disclosure of Potential Conflicts of Interest
R.K. Lai is a senior medical director at Zymeworks Biopharmaceuticals. No potential conflicts of interest were disclosed by the other authors.
Authors' Contributions
Conception and design: I. Atkins, B. Kinnersley, M. Simon, J.L. Bernstein, K.-H. Jöckel, C. Johansen, B. Melin, M.L. Bondy, R.S. Houlston
Development of methodology: I. Atkins, K. Labreche, J.L. Bernstein, C. Johansen, M.L. Bondy
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): I. Atkins, D. Il'yasova, J.E. Eckel-Passow, M.J. Schoemaker, M.M. Nöthen, J.S. Barnholtz-Sloan, A.J. Swerdlow, M. Simon, P. Rajaraman, S.J. Chanock, J. Shildkraut, J.L. Bernstein, P. Hoffmann, K.-H. Jöckel, R.K. Lai, E.B. Claus, S.H. Olson, C. Johansen, M.R. Wrensch, B. Melin, R.B. Jenkins, M. Sanson, M.L. Bondy, R.S. Houlston
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): I. Atkins, B. Kinnersley, K. Labreche, J.S. Barnholtz-Sloan, M. Simon, P. Rajaraman, J.L. Bernstein, R.K. Lai, R.B. Jenkins, M.L. Bondy, R.S. Houlston
Writing, review, and/or revision of the manuscript: I. Atkins, B. Kinnersley, Q.T. Ostrom, K. Labreche, J.E. Eckel-Passow, M.J. Schoemaker, M.M. Nöthen, J.S. Barnholtz-Sloan, A.J. Swerdlow, M. Simon, P. Rajaraman, S.J. Chanock, P. Hoffmann, K.-H. Jöckel, R.K. Lai, E.B. Claus, S.H. Olson, C. Johansen, B. Melin, R.B. Jenkins, M. Sanson, M.L. Bondy, R.S. Houlston
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): B. Kinnersley, G.N. Armstrong, R.B. Jenkins, M.L. Bondy
Study supervision: B. Kinnersley, G.N. Armstrong, K.-H. Jöckel, C. Johansen, M.R. Wrensch, M.L. Bondy, R.S. Houlston
Acknowledgments
I. Atkins was supported by a Wellcome Trust Summer Student bursary. In the UK, additional funding was provided by Cancer Research UK (C1298/A8362). The Glioma International Case-Control Consortium Study was supported by grants from the NIH, Bethesda, MD (R01CA139020, R01CA52689, P50097257, and P30CA125123). The UK Interphone Study was supported by the European Commission Fifth Framework Program “Quality of Life and Management of Living Resources” and the UK Mobile Telecommunications and Health Programme. The Mobile Manufacturers Forum and the GSM Association provided funding for the study through the scientifically independent International Union against Cancer (UICC). R code for power calculations was kindly provided by Lang Wu.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
