A Functional Variant on 9p21.3 Related to Glioma Risk Affects Enhancer Activity and Modulates Expression of CDKN2B-AS1.

Genome-wide association studies have identified SNPs associated with glioma risk on 9p21.3, but biological mechanisms underlying this association are unknown. We tested the hypothesis that a functional SNP on 9p21.3 affects activity of an enhancer, causing altered expression of nearby genes. We considered all SNPs in linkage disequilibrium with the 9p21.3 sentinel SNP rs634537 that mapped to putative enhancers. An enhancer containing rs1537372 exhibited allele-specific effects on luciferase activity. Deletion of this enhancer in GBM cell lines correlated with decreased expression of CDKN2B-AS1. Expression quantitative trait loci analysis using non-diseased brain samples showed rs1537372 to be a consistently significant eQTL for CDKN2B-AS1. Additionally, our analysis of Hi-C data generated in neural progenitor cells showed that the bait region containing rs1537372 interacted with the CDKN2B-AS1 promoter. These data suggest rs1537372, a SNP at the 9p21.3 risk locus, is a functional variant that modulates expression of CDKN2B-AS1.

INTRODUCTION

Gliomas, the most common form of brain cancer, account for around 13,000 deaths in the USA each year (Bondy et al., 2008; Ostrom et al., 2015). Gliomas are comprised of a heterogeneous group of tumors and are typically associated with a poor prognosis. Glioblastoma multiforme (GBM), the most common type of glioma, has a median overall survival of only 10–15 months (Bondy et al., 2008).

To date, at least 25 inherited risk variants have been identified by genome‐wide association studies (GWAS) of glioma (Kinnersley et al., 2018), accounting for approximately 30% of heritable risk (Melin et al., 2017). Discovering the biological mechanism underlying these risk variants has the potential to reveal novel insights into glioma development. However, deciphering the molecular underpinnings of risk alleles has proven to be challenging because few index GWAS single nucleotide polymorphisms (SNPs) are themselves functional. There is growing evidence that many functional SNPs that affect expression levels of target genes are in linkage disequilibrium (LD) with GWAS sentinel SNPs and map to enhancers or promoters (Ali et al., 2020; Biancolella et al., 2014; Fortini et al., 2014).

Several GWAS have identified 9p21.3 as a risk locus associated with glioma, and in particular with GBM (Kinnersley et al., 2015; Rajaraman et al., 2012; Shete et al., 2009; Wrensch et al., 2009). The largest GWAS meta‐analysis (cases: 12,496, controls: 18,190) confirmed the association between 9p21.3 and glioma and identified rs634537 as the lead SNP (NC_000009.11:g.22032152T>G) (P = 3.21 E−45) (Melin et al., 2017). These studies reveal that 9p21.3 is among the most consistently validated GWAS locus for GBM; however, a functional variant(s) within the region has yet to be identified.

The top GWAS SNP on 9p21.3, rs634537, maps to intron 1 of CDKN2B Antisense RNA 1 (CDKN2B‐AS1). It is unlikely that this SNP is functional since it does not map to any regulatory elements in astrocytes or neural stem cells (Figure 1a). We therefore tested the hypothesis that the functional SNP(s) in this region was in linkage disequilibrium (LD) with this top GWAS SNP. To identify candidate functional variants, we systemically screened SNPs in LD with rs634537 that intersected with regulatory elements/enhancers as catalogued in publicly available annotations. Each candidate enhancer/SNP region identified was tested in luciferase reporter assays for SNP‐dependent allele‐specific effects on enhancer activity. Using this approach, we identified a functional SNP, rs1537372 (NC_000009.11:g.22103183G>T), that mapped to an enhancer and had allele‐specific effects on enhancer activity. CDKN2B‐AS1 was identified as a target gene of this enhancer following CRISPR‐Cas9 deletion. We further demonstrated that this variant correlated with the expression of CDKN2B‐AS1 using expression quantitative trait loci (eQTL) analysis of non‐diseased brain tissues. Additionally, our Hi‐C analysis of data generated in neural progenitor cells showed that the bait region containing rs1537372 interacted significantly with the promoter of CDKN2B‐AS1. Taken together, our results support that rs1537372 is a functional SNP on 9p21.3 that modulates activity of an a local enhancer and that CDKN2B‐AS1, a gene previously implicated in glioma risk, is the most consistent target gene.

Figure 1

Region of chromosome 9 associated with glioblastoma risk and Allele‐specific enhancer activity of enhancer region 6. (A) Detailed view of the region defined by LD block with lead SNP rs634537, r > 0.2 in CEU population using the UCSC Genome Browser showing putative enhancer elements containing SNPs in LD with rs634537. SNPs in LD are noted below genes in the region. Histone ChIP‐Seq tracks for H3K27ac, H3K4me1 from normal human astrocytes (NHA), and H3K4me1 from MGG4 glioblastoma stem cells aligned below SNPs indicate potential enhancer elements. Region 2 and 3 denote regions with no enhancer activity in luciferase assays, Regions 1, 4 and 5 denote region that exhibited enhancer activity but were unaffected by haplotype. Region 6 (marked with an asterisk) denotes allele‐specific enhancer region, which includes rs15237372. Note that the size of the regions tested for enhancer activity are not to scale. (B) All enhancer regions seen in Figure 1a were cloned into a luciferase enhancer assay construct and tested for enhancer activity. Here we show data for enhancer region 6 that includes SNPs rs1537372, rs1537373, and rs1333042 haplotypes TGG, GGG, and GTA (four clones for each allele, three independent experiments for each cell line: (+) represent experiment 1, (×) represent experiment 2, and (*) represent experiment 3). The construct with the GGG and GTA alleles demonstrated statistically significantly higher activity than the TGG allele, as shown in box plots from U‐87 MG and U‐118 MG cells

Using our systematic approach of enhancer region analysis (detailed in the Supp. Materials and Methods section), we identified six putative enhancers (Figure 1a). All six putative enhancers were cloned separately into luciferase enhancer activity vectors (Biancolella et al., 2014; Fortini et al., 2014), in both directions, and enhancer activity of different alleles/haplotypes were independently tested following transfection of constructs into two GBM cell lines, U‐87 MG and U‐118 MG. Two of the six candidate enhancer regions (Region 2 including rs2069418 (NC_000009.11:g.22009698G>C) and Region 3 including rs1333037 (NC_000009.11:g.22040765C>T) and rs1360590 (NC_000009.11:g.22041443T>C) in Figure 1a) had no enhancer activity in luciferase assays (Supp. Figure S1). Four regions (Regions 1, 4, 5 and 6 in Figure 1a) demonstrated enhancer activity in both cell lines in at least one orientation (Supp. Figure S2). Summary statistics for luciferase assays included in Supp. Figure S1 and Supp. Figure S2 are listed in Supp. Table S1. There were a total of nine SNPs (r of ≥ 0.2 with rs634537) within these four enhancer regions: rs2518723 (NC_000009.11:g.21995882C>T) in Region 1; rs6475604 (NC_000009.11:g.22052734T>C), rs10757267 (NC_000009.11:g.22052810G>C), rs10965219 (NC_000009.11:g.22053687A>G), and rs7027048 (NC_000009.11:g.22053709A>G) in Region 4; rs1537370 (NC_000009.11:g.22084310C>T) in Region 5; and rs1537372 ((NC_000009.11:g.22103183G>T), rs1537373 (NC_000009.11:g.22103341T>G), and rs1333042 (NC_000009.11:g.22103813A>G) in Region 6. Only one SNP, rs1537372 in Region 6 (Figure 1a), demonstrated allele‐specific effects. rs157372 (minor allele frequency (MAF), T=0.4) is 71.1 kb away from the top GWAS SNP rs634537 (MAF, G=0.44) (pairwise LD r = 0.25, D'=0.69 in the CEU population). While the candidate enhancer in Region 6 demonstrated activity in luciferase assays in both directions, rs1537372 showed allele‐specific effects on enhancer activity in the forward orientation only (Figure 1b), and not the reverse orientation in either cell line (Supp. Figure S2). The fragment containing the rs1537372 G allele (the reference allele) correlated with higher activity than the fragment containing the T allele (the alternative allele) in both cell lines and all replicates tested.

To provide evidence that the candidate functional SNP rs1537372 identified in our cell‐based luciferase assays correlated with altered expression of genes mapping in cis, we used CRISPR‐Cas9 genome editing to delete the region containing rs1537372. The GBM cell lines U‐87 MG and U‐118 MG were chosen for these experiments because they demonstrated enhancer activity of the fragment containing rs1537372. We designed guide RNAs (gRNA1 and gRNA2) to delete an approximately 1 kb fragment containing the enhancer and SNP rs1537372 (Figure 2a). Each target sequence was cloned into a gRNA expression plasmid and used together with Cas9 expression plasmids to induce targeted deletion. U‐87 MG and U‐118 MG cells were transfected with either Cas9‐only expressing vector, Cas9 vector and gRNA empty vector, or Cas9 vector and guide RNA target vectors. Transfected cells were placed under geneticin selection for at least five days, after which time DNA was harvested and used to assess the deletion efficiency. Deletion efficiency was measured by PCR using primers designed to amplify an approximately 2 kb region across the putative enhancer containing the candidate functional SNP rs1537372 (PCR Forward and PCR Reverse in Figure 2a). The cell population with the unedited genome revealed a 2 kb band, while the population with edited cells revealed two bands: the unedited 2 kb fragment and the edited 1 kb fragment. Our results show that a large proportion of cells (though not all) in both cell lines were successfully edited (Figure 2b). RNA was isolated from the same cellular pools, and the expression of cis genes for which a TaqMan assay was available was tested by qPCR.

Figure 2

CRISPR‐Cas9 genome editing of enhancer region 6 on 9p21.3 and the associated gene expression changes. (A) Chromosome view of section of putative enhancer region 6 targeted by CRISPR‐Cas9 genome editing technique in the UCSC Genome Browser. Total region (~2 kb) represents part of putative enhancer region 6, containing the candidate functional SNP rs1537372, amplified by PCR (using Forward and Reverse primers). Region highlighted in grey represents region (~1 kb) targeted by CRISPR gRNAs (gRNA1 and gRNA2). (B) 1 kb band in DNA gel electrophoresis demonstrates the targeted deletion of putative enhancer region 6, containing SNP rs1537372, in U‐87 MG and U‐118 MG cells. NC: mock transfected parental cells; Cas9: cells transfected with Cas9 only (no guide RNAs); Cas9+Control: cells transfected with Cas9 vector and gRNA empty vector; Cas9+Target: cells transfected with Cas9 vector and guide RNA target vectors. (C) Genomic region corresponding to enhancer region 6 was targeted for deletion in U‐87 MG and U‐118 MG (Figure 2b) cells using CRISPR‐Cas9 technology. Pools of transfected cells were analyzed using qPCR and TaqMan gene expression assays for CDKN2B‐AS1 and TBP (control) custom assays in triplicate, in three independent experiments. NC: mock transfected parental cells; cas9: cells transfected with Cas9 only (no guide RNAs); Cas9+Control: cells transfected with Cas9 vector and gRNA empty vector; Cas9+Target: cells transfected with Cas9 vector and guide RNA target vectors. Targeting of region resulted in a decrease in CDKN2B‐AS1 expression levels. *P < 0.05; **P < 0.01; ***P < 0.001; and ****P < 0.0001 indicate the levels of significance

We quantified gene expression within 0.5 Mb either side of rs1537372 using TaqMan qPCR expression assays in CRISPR genome‐edited GBM cell lines and compared expression to mock CRISPR‐edited cells. We observed a significant reduction in the expression of CDKN2B‐AS1 (Figure 2c). No detectable levels of the following genes were observed in the GBM cell lines used in our study: CDKN2A, CDKN2B, MTAP, and DMRTA1. The two GBM cell lines used in our CRISPR deletion experiments, U‐87 MG and U‐118 MG, have been reported to harbor homozygous deletion in the 9p21.3 locus affecting the CDKN2A, CDKN2B, and MTAP genes (Garcia‐Claver et al., 2010), which is consistent with our results.

To control for off target effects and to further confirm the relationship between functional SNP rs1537372 and CDKN2B‐AS1 gene expression, we designed guide two new RNAs (gRNA3 and gRNA4) to delete an approximately 0.5 kb fragment containing the enhancer and SNP rs1537372 (Supp. Figure S3). U‐87 MG and U‐118 MG cells were electroporated with either Cas9 only, gRNAs only, or Cas9 and gRNA target vectors. DNA was harvested and used to assess the deletion efficiency. As before, deletion efficiency was measured by PCR using primers designed to amplify an approximately 2 kb region across the putative enhancer containing the candidate functional SNP rs1537372 (PCR 2 Forward and PCR Reverse 2 in Supp. Figure S3). The cell population with the unedited genome revealed a 2 kb band, while the population with edited cells revealed two bands: the unedited 2 kb fragment and the edited 1.5 kb fragment. Our results showed that a large proportion of cells (though not all) in both cell lines were successfully edited (Supp. Figures S3B). RNA was isolated from the same cellular pools, and the expression of CDKN2B‐AS1 was determined. We observed a significant reduction in the expression of CDKN2B‐AS1 (Supp. Figure S4), confirming our prior results.

eQTL mapping of the genes within the a 0.5 Mb window in either direction of the functional SNP rs1537372 using non‐diseased brain samples revealed CDKN2B‐AS1 as the only gene significantly associated with allelic dosage at rs1537372 (β = 0.11 (standard error (SE) 0.03), P = 4.45 E−04, FDR Q value 0.01, Supp. Table S2). Further eQTL analysis of genes outside of the 0.5 Mb window of rs1537372 did not reveal any significant associations (data not shown). Direct Imputation of summary STatistics (DIST) imputed rs1537372 with an imputation INFO score of 0.93 in GBM. The GWAS Z score was −7.84 (P = 4.06 E−09). COLOC showed PP3 was 0 and PP4 was 0.81 for the target gene CDKN2B‐AS1. This suggests rs1537372 has high probability of being a causal variant for both the eQTL and glioma GWAS signals in GBM. Additionally, Hi‐C analysis of data generated in neural progenitor cells showed that SNP rs1537372 interacted significantly with the promoter of CDKN2B‐AS1 (FitHiC Q value = 4.57 E−04) and 23 other genes (Supp. Table S3). Supp. Figure S5 illustrates the interaction of glioma functional SNP rs1537372 with target genes in the 9p21.3 locus.

To summarize, we provide evidence that rs1537372 is a functional SNP on 9p21.3 mapping to a risk enhancer using cell‐based enhancer activity assays. We show that rs1537372 had an allele‐specific effect on enhancer activity in the forward direction only. We previously noted unidirectional, rather than bi‐directional, allele‐specific effects on enhancer activity in colorectal cancer and glioma functional studies (Ali et al., 2020; Biancolella et al., 2014; Fortini et al., 2014). CRISPR deletion of the risk enhancer demonstrated that this SNP affects glioma risk potentially through altered expression of CDKN2B‐AS1. We applied complementary methods of eQTL mapping and Hi‐C interaction data that also support CDKN2B‐AS1 being the most consistent target gene across orthogonal methodologies.

The functional SNP identified in our study, rs1537372, maps to intron 3 of the CDKN2B‐AS1 gene. CDKN2B‐AS1, also known as ANRIL (Antisense Noncoding RNA In The INK4 Locus), encodes a functional RNA molecule that epigenetically represses tumor suppressors CDKN2A and CDKN2B through interactions with polycomb repressive complex‐1 (PRC1) and −2 (PRC2) (Aguilo et al., 2011). CDKN2B‐AS1 has been shown to be implicated in several diseases, including cardiovascular diseases (Chan et al., 2013; Holdt & Teupser, 2012), diabetes (Cheng et al., 2011), and cancers such as endometrial cancer, colorectal cancer, renal cell carcinoma, and osteosarcoma (Gu et al., 2013). A link between polymorphisms in the CDKN2B‐AS1 gene and glioma risk has been reported by several groups. Previously, an association between risk variants mapping to the CDKN2B‐AS1 gene and the risk of glioma has been found (Wibom et al., 2015). CDKN2B‐AS1 has been shown to promote glioma cell growth, migration, and invasion by Dong et al. (Dong et al., 2018). Further, CDKN2B‐AS1 knockdown resulted in decreased proliferation and migration in glioma cells (Gao et al., 2020). These studies, as well as our own finding, implicate CDKN2B‐AS1 as an important gene involved in glioma risk, pathogenesis and progression. Importantly, the functional SNP rs1537372 maps to an enhancer region within the CDKN2B‐AS1 gene and regulates its expression. These data are consistent with our eQTL and Hi‐C analyses, and together provide strong evidence that CDKN2B‐AS1 is a key gene that mediates risk of GBM on 9p21.3.

We cannot discount the possibility that there are additional functional SNPs on 9p21.3 that may regulate gene expression, and this includes alternate mechanisms such as microRNA binding or alternative splicing. While we used a number of publicly available datasets to identify candidate enhancers, these may not have captured all relevant enhancers across this region. Our in vitro assessment of enhancer activity was conducted in only two GBM cell lines, and we cannot exclude the possibility that a relevant enhancer(s) in this region is inactive in the cell lines tested but active in other cell lines. Additionally, the GBM cell lines that were used for the CRISPR deletion experiments may have limitations, because additional target genes of this enhancer may be seen in other glioma cell lines, neural stem cells or normal astrocytes. U‐87 MG and U‐118 MG cell lines have been reported to have deletions in the 9p21.3 locus, with the loss of genes CDKN2A, CDKN2B, and MTAP. However, U‐87 MG and U‐118 MG cell lines were selected for our CRISPR deletion experiments as the candidate enhancers tested in our study showed activity in these cell lines only. Activity was not detected in the other GBM cell lines that we tested: LN‐229, LN‐19, and A172. As a result of our using these cell lines, we cannot assess the relationship between the functional SNP and expression of the CDKN2A, CDKN2B, and MTAP genes. Another limitation of these CRISPR experiments is that we have not determined the impact of SNP rs1537372 directly, but of the enhancer containing that SNP. The functional SNP rs1537372 identified in this study is in relatively low LD with sentinel SNP rs634537 (r = 0.25 in the CEU population). This is not unprecedented. It has been suggested that a correlation r threshold value of 0.2 or even less may be needed when determining LD structure, as the strength of the correlation between the sentinel SNP and the causal SNP may be low (Freedman et al., 2011). Baskin et. al., identified a candidate functional SNP in LD with the GWAS sentinel SNP on 11q23.3 using r >= 0.2. The study identified rs73001406 as a candidate functional variant with PHLDB1 and DDX6 as potential target genes. SNP rs73001406 was in a relatively low LD with the glioma risk SNP rs498872 on 11q23.3 (r = 0.38) (Baskin et al., 2015). Finally, our eQTL analyses were based upon bulk expression estimates that represent the composite results of a number of CNS cells. Single cell data generated in astrocytes is still lacking, and single cell eQTL datasets generated to date have less power for discovery, given the high cost to generate data (van der Wijst et al., 2018). Despite this we believe we provide strong supportive evidence for the identification of rs1537372 as a functional SNP relevant to gliomagenesis on 9p21.3.

In summary, we report identification and characterization of a functional SNP, rs1537372, that affects activity of an enhancer on 9p21.3 that leads to modulated expression of CDKN2B‐AS1, a gene implicated in glioma risk. Further studies to evaluate the molecular mechanism of this gene on glioma growth and development may lead to novel therapeutic strategies in the future.

WEB RESOURCES

Enhancer region analysis was performed using University of California Santa Cruz (UCSC) Genome Browser: http://genome.ucsc.edu/.

ChIP‐Seq peaks for histone modifications were obtained from the Encyclopedia of DNA Elements (ENCODE): https://www.encodeproject.org/.

CRISPR gRNAs (guide RNAs) were designed using https://www.crisprscan.org/.

eQTL mapping was done using CommonMind Consortium (CMC) release 1: http://commonmind.org/.

The Glioma International Case–Control study (GICC) GWAS summary statistics is available at European Genome‐phenome Archive (http://www.ebi.ac.uk/ega/; accession number: EGAS00001003372).

CONFLICT OF INTEREST

The authors declare that there are no conflict of interests.

AUTHOR CONTRIBUTION

G.C. and R.K.L.: conception and design. M.W.A., S.P. and C.H.D. performed and analyzed luciferase assays. M.W.A., C.K., and M.A. conceived and performed CRISPR‐Cas9 experiments. C.P.K.P. performed eQTL mapping. M.W.A and M.D. performed and analyzed qPCR experiments. M.W.A., G.C., C.P.K.P and R.K.L. wrote the manuscript with input from all authors. All authors read and approved the manuscript. All authors consent to publication.

Supplementary information.

Click here for additional data file.