Glioblastoma is the most prevalent primary malignant brain tumor in adults and is characterized by poor prognosis and universal tumor recurrence. Effective glioblastoma treatments are lacking, in part due to somatic mutations and epigenetic reprogramming that alter gene expression and confer drug resistance. To investigate recurrently dysregulated genes in glioblastoma, we interrogated allele-specific expression (ASE), the difference in expression between two alleles of a gene, in glioblastoma stem cells (GSC) derived from 43 patients. A total of 118 genes were found with recurrent ASE preferentially in GSCs compared with normal tissues. These genes were enriched for apoptotic regulators, including schlafen family member 11 (SLFN11). Loss of SLFN11 gene expression was associated with aberrant promoter methylation and conferred resistance to chemotherapy and PARP inhibition. Conversely, low SLFN11 expression rendered GSCs susceptible to the oncolytic flavivirus Zika. This discovery effort based upon ASE revealed novel points of vulnerability in GSCs, suggesting a potential alternative treatment strategy for chemotherapy-resistant glioblastoma.
Assessing allele-specific expression reveals genes with recurrent cis-regulatory changes that are enriched in glioblastoma stem cells, including SLFN11, which modulates chemotherapy resistance and susceptibility to the oncolytic Zika virus.
Glioblastoma ranks among the most lethal of human malignancies with current therapies only offering palliation (1). Reasons for treatment failure are myriad, with tumor heterogeneity at the genetic and transcriptional levels contributing to the malignancy of glioblastoma (2, 3). Glioblastoma displays a functional cellular hierarchy with stem-like, self-renewing glioblastoma stem cells (GSC) residing at the apex (4, 5). GSCs contribute to resistance to chemotherapy and radiotherapy, neoangiogenesis, invasion into normal brain, and escape from the immune system (6–8). Therefore, targeting GSCs may improve current glioblastoma management and extend the lives of patients.
Although glioblastoma is one of the most deeply characterized solid tumors, precision medicine has not significantly benefited most neuro-oncology patients. Most studies and targeted therapeutic strategies have thus far focused on protein-coding mutations. However, many important mutations lie in noncoding DNA where they function by perturbing gene regulation. Noncoding mutations likely help drive glioblastoma tumorigenesis and drug resistance but are more challenging to identify and impact gene regulation in many different ways. Gene dysregulation can be caused by copy-number alterations (CNA; refs. 9, 10), as well as by mutations that affect splice sequences (11), untranslated regions (UTR; ref. 12), insulators (13, 14), promoters (15, 16), and enhancers (12, 17). Moreover, genes are often regulated by multiple enhancers, which can be located hundreds of kilobases away from their targets (18). Regulatory mutations are, therefore, diverse and frequently spread over very large regions of the genome. As a result, standard recurrence analyses that identify driver mutations in coding sequences are unlikely detect many important regulatory mutations in cancer genomes. Alternate approaches to discover regulatory mutations can by stymied by the myriad of noncoding mutations whose function is difficult to predict from sequence alone. Rather than relying on interpretation of noncoding mutations, unbiased identification of genes with altered regulation can pinpoint functionally important genes that are unlikely to be discovered through other methods.
Here, we leveraged allele-specific expression (ASE) as a new approach to interrogate recurrently dysregulated genes in glioblastoma. ASE is the difference in expression between two alleles of a gene and can be estimated from mapped RNA sequencing (RNA-seq) reads that overlap heterozygous variants within exons (19, 20). Specifically, by counting the number of RNA-seq reads that contain the reference or alternate allele, the relative expression of each copy of the gene can be measured. Unlike standard RNA-seq expression levels, ASE is particularly sensitive to cis-acting mutations because it is generally unaffected by trans-acting or environmental effects that impact both alleles equally. Heterozygous sites within a gene's exons are required to obtain a readout of ASE; however, these heterozygous sites are not necessarily the cause of a gene's ASE. For example, the effect of a distal cis-regulatory mutation on gene expression can be measured using heterozygous variants within a gene, even when the distal mutation has not been directly observed by sequencing. Thus, a major advantage of ASE is that it can identify genes that are dysregulated by cis-acting regulatory mutations, even when the specific identities of the regulatory mutations are unknown. However, ASE is not solely caused by somatic mutations and can also result from common germline polymorphisms (21), imprinting (22), or random monoallelic expression (23). Therefore, to discover genes that are dysregulated by cis-acting factors such as regulatory mutations in cancer, the frequency of ASE in disease samples must be compared with a panel of normal samples. This approach has recently been used to identify new pathogenic genetic variants in muscle disease (20) and oncogenic mutations in T cell–acute lymphoblastic leukemia that would be difficult to identify using traditional techniques (14).
Based on this background, we hypothesized that a discovery effort based upon ASE could reveal novel points of fragility in the most resistant tumor cells, the GSCs. GSCs maintained in serum-free conditions maintain both genetic and transcriptional signatures found in the tumors from which they were derived, while removing the contaminating nontransformed cells that complicate genetic discovery. Here, we interrogated ASE in 43 patient-derived GSCs and compare the frequency of ASE to normal tissues to reveal novel dysregulated molecular targets that promote drug resistance and confer therapeutic vulnerabilities.
Materials and Methods
Derivation and maintenance of GSCs
Glioblastoma samples were obtained from surgical resection from patients at Duke University or Case Western Reserve University with informed consent in accordance with the Cleveland Clinic Institutional Review Board-approved protocol 090401. Prior to use, all samples were reviewed and verified by a neuropathologist. All patient studies were conducted in accordance with the Declaration of Helsinki. GSC23 was acquired via a material transfer agreement from The University of Texas MD Anderson Cancer Center (Houston, TX). GSCs were cultured in Neurobasal media (Invitrogen) supplemented with B27 without vitamin A (Invitrogen), EGF, and bFGF (20 ng/mL each; R&D Systems), sodium pyruvate, and glutamax. Short tandem repeat analyses were performed to authenticate the identity of each tumor model used in this article both annually and prior to high-throughput sequencing. Cells were stored at −160°C when not being cultured. To minimize cell culture–based artifacts, patient-derived xenografts were produced and propagated as a renewable source of tumor cells for study. Cultured cells underwent PCR testing for mycoplasma routinely every 6 months, and all RNA-seq data were evaluated for mycoplasma content during quality control.
Exome-seq reads were aligned to the GRCh37 (hg19) assembly of the reference genome using BWA-MEM with default parameters. Mapped reads were filtered for mapping quality score ≥ 30 and duplicate reads were removed using samtools (1.9, RRID:SCR_00215; ref. 24). Genotypes were generated for each individual using GATK's HaplotypeCaller and jointly processed using the GenotypeGVCFs function in GATK (4.1.1, RRID:SCR_001876). Following genotyping, single-nucleotide variants (SNV) were extracted and filtered using the variant quality score recalibration (VQSR) method in GATK.
RNA-seq alignment and processing
RNA-seq reads were aligned end-to-end to the GRCh37 (or hg19) assembly of the reference genome using STAR (2.5.3a, RRID:SCR_004463). Mapped reads were filtered using mapping quality score ≥ 20 and duplicate reads were removed using samtools (1.9). Read counts for GENCODE (v28) genes were computed using FeatureCounts (1.6.3, RRID:SCR_012919; ref. 25) and fragments per kilobase per million (FPKM) values were estimated using DESeq2 (1.14.1; ref. 26). For downstream analysis, the FPKM values were quantile normalized, and then converted to z-scores by mean-centering and standardizing across samples.
Estimating allele-specific expression
RNA-seq read alignments were corrected for mapping bias and allele-specific read counts at heterozygous positions were collected using WASP (27). Heterozygous sites covered by at least 10 RNA-seq reads were retained for allele-specific expression analysis.
We estimate a for each gene by maximum likelihood under the alternative model of allele imbalance. Then we use a likelihood ratio test to compare the alternative model to the null model of no allele imbalance (i.e., with a fixed to a = 0). We correct the P values from the likelihood ratio test for multiple testing using the Benjamini–Hochberg method. To make it clear we are referring to allele imbalance in RNA-seq reads, we refer to a as aRNA in the rest of the manuscript.
Enrichment compared with a normal whole-blood and brain tissues from GTEx
To generate a reference profile of ASE in normal samples, we obtained RNA-seq data for 369 whole blood and 216 brain samples distributed across 13 brain regions from the GTEx consortium (19) and analyzed these data using our ASE algorithm described in the previous section. To discover genes that were enriched for ASE in GSCs, we compared the frequency of samples with significant ASE for each gene between GSCs and whole blood using Fisher exact test (FET) and adjusted the resulting P values using the Benjamini–Hochberg procedure. We only analyzed genes that were testable for ASE (i.e., had ≥ 1 heterozygous site with ≥10 reads) in both GSCs and whole-blood tissues. For further analysis of ASE frequency for individual genes, we also compared estimated allele imbalance from our model (aRNA) among whole blood, 13 brain tissues, and GSCs. Manhattan plots for enriched genes were generated using ggbio (1.30). Gene ontology enrichment analysis for genes showing recurrent ASE in GSCs was carried out using topGO (2.34, RRID:SCR_014798; ref. 28).
Association between DNA methylation and gene expression
We downloaded precomputed genome-wide methylation data for 43 GSCs from Gene Expression Omnibus (GSE119774, RRID:SCR_005012; ref. 29). These methylation data were generated using the Illumina Infinium Epic Methylation Array. In this assay, DNA methylation levels at CpG sites are represented by β, which is the ratio of the methylated (C) to unmethylated (T) signal. We annotated the CpG probe positions based on GENCODE (v28, RRID:SCR_014966) genes and computed the mean β values for promoter regions (i.e., 1 kb upstream to 500 bp downstream of annotated transcription start sites; βpromoter). To discover ASE genes that may be dysregulated by aberrant DNA methylation, we computed Spearman rank correlation between βpromoter and normalized gene expression. We corrected correlation P values for multiple testing using the Benjamini–Hochberg procedure. For this analysis, we only considered genes with ≥ 3 CpG probes mapping to their promoter regions.
Analysis of H3K27ac chromatin immunoprecipitation and sequencing data
We downloaded H3K27ac chromatin immunoprecipitation and sequencing (ChIP-seq) data for 35 GSCs from Gene Expression Omnibus (GSE119755, RRID:SCR_005012; ref. 29). ChIP-seq reads were aligned to the GRCh37 (hg19) assembly of the reference genome using BWA-MEM with default parameters. The mapped reads were filtered using mapping quality score of ≥20, and duplicate reads were removed using samtools (1.9). H3K27ac peaks were called using MACS2 in paired-end mode with custom parameters (–nomodel –extsize 200 –qvalue 0.05; ref. 30). To generate a unified set of test regions, we divided the genome in 1 kb nonoverlapping genomic bins and kept the bins that overlapped MACS2.0 peaks in at least one GSC. We recounted the reads mapping to these genomic bins using exomeCopy (1.28.0, RRID:SCR_001276) and calculated FPKM using DESeq2 (1.22.2, RRID:SCR_000154; ref. 26). FPKM measurements were quantile normalized and mean-centered for downstream analysis.
The 1 kb genomic bins generated from H3K27ac ChIP-seq peaks were annotated using the ChIPSeeker (1.18.0) package. To discover cis-regulatory elements (CRE), we selected all distal intergenic and intronic genomic bins located within 100 kb of promoters (i.e., 1 kb upstream and 500 bp downstream of transcription start sites) of ASE genes. We performed a Spearman correlation analysis between normalized coverage for bins and normalized expression for genes to identify CREs. We corrected the P values for multiple testing using the Benjamini–Hochberg procedure.
TRIzol reagent (Sigma-Aldrich) was used to isolate total cellular RNA from cell pellets, and the qScript cDNA Synthesis Kit (Quanta BioSciences) was used for reverse transcription. Quantitative real-time PCR was performed using SYBR-Green PCR Master Mix (Thermo Fisher Scientific) on an Applied Biosystems 7900HT cycler.
Cells were collected and lysed in RIPA buffer (50 mmol/L Tris-HCl, pH 7.5; 150 mmol/L NaCl; 0.5% NP-40; 50 mmol/L NaF with protease inhibitors) and incubated on ice for 30 minutes. Lysates were centrifuged at 4°C for 15 minutes at 14,000 rpm, supernatant was collected, and protein concentration was confirmed using a Pierce BCA protein assay kit (Thermo Scientific, cat #23225). Equal amounts of protein samples were mixed with SDS Laemmli loading buffer, boiled for 10 minutes, electrophoresed using NuPAGE Bis-Tris gels and transferred onto PVDF membranes. Membranes were blocked for 1 hour with TBS-T plus 5% nonfat dry milk, then incubated in primary antibodies overnight at 4°C. Blots were washed 3 times for 5 minutes each with TBS-T and then incubated in TBS-T plus 5% milk for 1 hour with appropriate secondary antibodies. Blots were imaged using Bio-Rad Image Lab software and processed using Adobe Illustrator (RRID:SCR_010279) to create the figures. The following antibodies were used for Western blot: SLFN11 (Santa Cruz Biotechnology; cat #SC-374339, RRID:AB_10989536) and HRP-conjugated GAPDH (Proteintech; cat #HRP-60004, RRID:AB_2737588).
Lentiviral constructs expressing shRNAs directed against SLFN11 (Sigma TRCN, TRCN0000152057) or a nontargeting control shRNA (TRCN0000231489) with no targets in the human genome were obtained from Sigma-Aldrich. The SLFN11 expression vector was obtained from VectorBuilder along with an empty vector control with the same lentiviral backbone. 293T cells (ECACC; cat #12022001, RRID:CVCL_0063) were used to generate lentiviral particles by cotransfection of packaging vectors psPAX2 (RRID:Addgene_12260) and pMD2.G (RRID:Addgene_12259) using a standard PEI transfection method in DMEM plus 1% penicillin/streptomycin. GSCs were transduced with the lentiviral constructs, and selection was started 48 hours later using 1 μg/mL of puromycin for 72 hours, at which times, cells were assayed for SLFN11 expression.
In vitro treatment and cell viability
For in vitro cell viability assays, 2,000 cells/well for individual or 5,000 cells/well for combinatorial drug studies were plated in a 96-well plate. Cells were then treated 24 hours later with temozolomide (Selleck Chem; cat #S1237), olaparib (Selleck Chem; #S1060), both drugs, or DMSO at an equivalent percent volume to the highest drug concentration. Cell viability was assayed 4 days later following a 12-hour incubation with alamarBlue (Thermo Fisher; cat #DAL1025) and detected using a fluorescence-based plate reader. For Zika virus studies, 5,000 cells/well were infected with Dakar 41519 strain ZIKV at a multiplicity of infection of 5 (MOI = 5; ref. 31). Viability was assayed at 3 days of post infection using CellTiter-Glo according to the manufacturer's instructions.
Tests for ASE were performed using likelihood ratio tests as described above. Frequencies of samples with significant ASE were compared using FET, as described above. For analysis of all in vitro data, comparisons between more than two samples were performed using ANOVA followed by Tukey HSD or using Student t test. The number of replicates and specific statistical tests used are indicated in the figure legends.
Data availability statement
RNA-seq (GSE119834), H3K27ac ChIP-seq (GSE119755), and human methylation 450K array data (GSE119774) for GSCs used in the paper were previously published and can be downloaded from GEO (GSE119776; ref. 29).
Glioblastoma target gene discovery leveraging recurrent allele-specific expression
To discover genes with ASE, we implemented a statistical model to estimate allele imbalance for every gene (Fig. 1A). Our systematic methodology utilizes RNA-seq reads overlapping heterozygous sites and reduces false positives by accounting for technical sources of variation, including sequencing errors, genotyping errors, and overdispersion of RNA-seq read counts. We quantified ASE with RNA allelic imbalance (aRNA), which is the difference between the reference allele proportion and the expected value of 0.5. We applied our model to RNA-seq data from 43 patient-derived GSCs representing a diverse population of patients (age, biological sex, etc.) and tumor features (genetics, transcriptional subgroups, etc.) that we previously reported (29). We identified 5,808 genes with significant ASE in at least one GSC under a false discovery rate (FDR) of 10% (Supplementary Table S1). ASE was restricted to a single sample for most genes (3,948 out of 5,808). However, 1,860 genes showed ASE in 2 or more GSCs, and 298 genes showed highly recurrent ASE that was present in 5 or more GSCs (Fig. 1B).
To discover genes that showed recurrent ASE specific to GSCs relative to normal tissues, we compared the frequency of ASE for each gene in GSCs to both normal whole-blood and normal brain samples from the Genotype Tissue Expression (GTEx) project using the FET. Whole blood has, by far, the largest number of samples in GTEx, so we used it as the initial reference tissue; however, we performed subsequent comparisons with ASE in 13 different brain tissues to account for tissue-specific imprinting, as described below. Under an FDR of 10%, 118 genes displayed significantly enriched ASE in GSCs (Supplementary Table S2). To illustrate the power of this approach, we examined the ASE patterns of the noncoding RNA gene H19, which is maternally imprinted (32), and RHOB, which is important for glioblastoma tumorigenesis (33, 34). Imprinting causes loss of expression of the allele inherited from one parent, and as expected for an imprinted gene, H19 showed ASE in almost all normal and tumor samples (Fig. 1C, top). In contrast, ASE of RHOB was restricted to GSCs (Fig. 1C, bottom).
Recurrent ASE in cancer genomes can be caused by frequent CNAs or loss-of-heterozygosity. In the presence of large CNAs, many genes are expected to exhibit ASE, but these genes would be clustered into the genomic regions that undergo frequent CNAs. However, in our data set, recurrent ASE genes were distributed throughout the genome, suggesting that their dysregulation was not due to large-scale CNAs, but was instead caused by localized cis-regulatory mutations or epigenetic changes (Supplementary Fig. S1A).
To examine the biological function of the 118 genes with recurrent ASE in GSCs, we performed Gene Ontology (GO) enrichment analysis, revealing overrepresentation of recurrent ASE genes involved in regulation of cell cycle and apoptosis (FDR ≤ 5%). Twenty-eight genes with recurrent ASE were associated with the Biological Process GO category “programmed cell death,” whereas only 13 genes were expected by chance (Fig. 1D). These genes included the kinase IP6K2 (35, 36), which exhibited a marked enrichment of ASE in GSCs compared with both whole-blood (FET P = 0.001; FDR-adjusted P = 0.06) and brain tissues from GTEx (Supplementary Fig. S1B and Supplementary Table S2).
Allele-specific gene expression is associated with H3K27ac marks at distal regulatory elements
Cis-acting pathogenic variants impact gene expression by disrupting regulatory sequences, such as enhancers (17, 37). To illuminate the cis-acting mechanisms that underlie ASE in GSCs, we tested whether the expression of the 118 ASE genes was associated with the activity of nearby CREs. We identified putative CREs within 100 kb of the ASE gene promoters using histone H3 lysine 27 acetyl (H3K27ac) chromatin immunoprecipitation followed by deep sequencing (ChIP-seq) data that we previously generated for the GSCs (29). We divided the genome into 1 kb genomic bins and labeled bins that overlapped H3K27ac peaks in at least one GSC as putative CREs.
To connect CREs with genes, we correlated normalized H3K27ac levels with gene expression, focusing on distal intergenic and intronic elements. We used gene expression rather than ASE for this purpose, because ASE can only be estimated in samples that carry at least one heterozygous SNP within exons. Furthermore, ASE does not detect regulatory changes that affect both alleles equally. Of the 118 ASE genes, 56 had a significant Spearman rank correlation with the activity of a distal CRE (FDR ≤ 5%). In many cases, a single gene was associated with the activity of multiple CREs, such that 227 CRE bins were associated with the expression of these 56 genes (Supplementary Table S3). Thus, gene expression of many ASE genes was associated with activity levels of distal regulatory elements, as measured by H3K27ac, suggesting that cis-acting mutations within these CREs are a plausible and potentially common mechanism for the dysregulation of these genes in glioblastoma.
One example of an ASE gene that is associated with multiple CREs is NOTCH1. ASE of NOTCH1 was enriched in GSCs compared with both whole-blood (FET P = 5.09e−5; FDR-adjusted P = 0.0094) and brain samples, and we observed extreme biases in reference and alternate allele proportions at all heterozygous sites (Fig. 2A and B). NOTCH1 gene expression correlated with H3K27ac levels of 29 nearby CREs (Fig. 2C and D; Supplementary Table S3). To confirm that these CREs were more strongly associated with NOTCH1 expression than expected by chance, we created an empirical null distribution for each CRE by correlating H3K27ac levels with the expression of 1,000 randomly selected genes. Using these null distributions, 75% (22 out of 29) CREs remained significantly correlated with NOTCH1 gene expression (empirical P ≤ 5%). These regions may contain cis-regulatory noncoding mutations that lead to NOTCH1 gene dysregulation and may be excellent targets for future regulatory screens dissecting the regulation of NOTCH1 expression in GSCs.
Promoter methylation is associated with the expression of ASE genes
Aberrant DNA methylation of gene promoters is associated with widespread gene expression changes and chemotherapy resistance in brain tumors. In gliomas, MGMT CpG-rich promoter methylation is associated with decreased expression and improved response to temozolomide (TMZ) treatment (38). To determine whether genes with recurrent ASE were associated with aberrant DNA methylation, we analyzed the CpG methylome of the GSCs (29). We estimated the promoter methylation of the 118 ASE genes (βpromoter) and computed Spearman rank correlations with normalized gene expression. We identified 30 genes that displayed correlation between promoter methylation and gene expression at FDR ≤ 10% (Supplementary Table S4), of which, 16 overlapped with the 56 genes identified above with correlated gene expression and H3K27ac levels at CREs. Thus, 70 of the 118 ASE genes were associated with promoter methylation, CRE activity or both, suggesting possible mechanisms for their dysregulation.
SLFN11 promoter methylation associates with its gene expression
One of the 30 genes with correlated gene expression and promoter methylation was SLFN11 (rho = −0.68, FDR corrected P = 0.0001; Fig. 3A–C). SLFN11 is notable because it inhibits DNA replication and promotes cell death in response to DNA damage (39, 40). Loss of SLFN11 causes resistance to poly ADP ribose polymerase (PARP) inhibitors in small cell lung cancer, suggesting that it may be an important marker of chemotherapy resistance (41). ASE of SLFN11 was highly enriched in GSCs with 4 out of 10 testable GSCs exhibiting ASE, compared with 0 out of 159 testable normal whole blood tissues from GTEx (FET P = 6.4e−6; FDR-adjusted P = 2.2e−3; Supplementary Table S2). Although this gene showed enrichment of ASE in GSCs compared with whole blood samples, we detected ASE in a small number of normal brain tissue samples (4 out of 224; Fig. 3D). This suggests that rare germline variants or somatic events, such as mutations or DNA methylation, may affect SLFN11 expression in some phenotypically normal individuals. In the four GSCs with significant SLFN11 ASE the RNA-seq read counts show strong allelic bias at multiple SNVs at different locations within the gene, which rules out the possibility that the observed ASE is due to genotyping errors or read mapping artifacts that would be more likely to affect a single site (Fig. 3E).
Based on the promoter methylation and gene expression of SLFN11, GSCs can be divided in three distinct classes: (i) GSCs with high methylation and low expression; (ii) GSCs with intermediate methylation and intermediate expression; and (iii) GSCs with low methylation and high expression (Fig. 3A). Four GSC samples with ASE of SLFN11 had intermediate methylation, which is consistent with methylation and reduced expression of one allele (Fig. 3A and B). However, samples like GSCs 2907 and 007B had high promoter methylation levels (>60%) and very low expression of SLFN11 (Fig. 3A). Under these circumstances, genes would not be found by ASE because the expression of both alleles is reduced. These results demonstrate that low expression of SLFN11 in GSCs is associated with increased promoter methylation and that the samples with detectable ASE likely have methylation of one allele, but not the other.
SLFN11 augments chemotherapy resistance in GSCs
SLFN11 regulates cellular responses to DNA damaging agents (39, 40). Therefore, we hypothesized that SLFN11 expression in GSCs would be associated with chemotherapeutic resistance to an alkylating agent, TMZ, and a PARP inhibitor, olaparib. To test this hypothesis, we utilized four patient-derived GSCs: two with high SLFN11 expression and no evidence for ASE (839 and MNK1) and two with low SLFN11 expression and strong ASE (2012 and 1552). We confirmed both SLFN11 mRNA and SLFN11 protein abundance by RT-PCR (Fig. 4A) and immunoblot (Fig. 4B). We treated cells with drug concentrations ranging from 0 to 1,000 μmol/L for TMZ and 0 to 50 μmol/L for olaparib, then measured drug effects on cell survival to generate a concentration–response matrix where an effect of 100% corresponds to complete killing of all cells and 0% corresponds to no difference in cell survival (Fig. 4C and D; Supplementary Fig. S2). We estimated synergy between the two drugs using SynergyFinder 2.0 (42). Two GSCs with low expression and ASE of SLFN11 had reduced responses to drug treatment compared the GSCs with high expression of SLFN11. The maximum response for both ASE GSCs ranged from 40% to 50%, whereas the GSCs with high expression of SLFN11 had responses ranging from ∼70% to 80% (Fig. 4D). The two ASE GSCs also had lower drug synergy scores for the combination of the drugs (Fig. 4E).
To directly test whether SLFN11 expression affects chemotherapeutic drug sensitivity in GSCs, we also performed knockdown (KD) and overexpression (OE) experiments. Specifically, we performed KD of SLFN11 in MNK1 GSCs, which have high baseline expression of SLFN11 (Fig. 5A), using two nonoverlapping shRNAs. After confirming shRNA-mediated knockdown of SLFN11 in MNK1 GSCs, we measured cell survival in response to treatment with TMZ and olaparib. SLFN11 KD in MNK1 cells reduced cell death in response to drug treatment (Fig. 5A and B). Conversely, overexpression of SLFN11 in 2012 GSCs, which have low baseline expression of SLFN11, increased cell death in response to combinatorial drug treatment (Figs. 5C and D). Knockdown and overexpression of SLFN11 similarly changed sensitivity to single-drug treatments, although not significantly so (0.05 < P < 0.1; Supplementary Fig. S3). Synergistic drug response to TMZ and olaparib treatment decreased upon SLFN11 KD and increased with OE (Fig. 5E). Thus, SLFN11 expression is causally linked to sensitivity to the chemotherapeutic drugs TMZ and olaparib.
GSCs with low expression of SLFN11 are sensitive to Zika virus
SLFN11 may be involved in cellular response to viral infection. SLFN11 is upregulated following virus-induced type I interferon response and restricts flavivirus replication in human tumor cell lines (43–45). Oncolytic viruses that infect the central nervous system can be leveraged to treat brain tumors (46). We recently demonstrated that Zika virus, which is a member of the flavivirus genus of RNA viruses, preferentially infects and kills GSCs compared with differentiated tumor cells and normal neuronal cells (47). Based on this background, we hypothesized that tumor cells with promoter methylation of SLFN11 would be more susceptible to oncolytic destruction by Zika because they are unable to increase SLFN11 expression in response to interferon stimulation (47). Under this hypothesis, interferon exposure would induce upregulation of most interferon-responsive genes but would fail to induce SLFN11 expression in cells with SLFN11 promoter methylation.
To test the above hypothesis, we examined the association between SLFN11 gene expression and type I interferon response by analyzing gene expression data from 669 glioblastoma tumors obtained from The Cancer Genome Atlas (48). To quantify interferon response in each sample, we computed a single-sample gene set enrichment analysis (ssGSEA) score using the interferon alpha response hallmark gene set (49, 50). Glioblastoma had higher gene expression of SLFN11 compared with grade II and III gliomas (Fig. 6A), and SLFN11 gene expression correlated with the interferon alpha (IFNα) ssGSEA score (R = 0.4, P < 2.2e−16), consistent with induction of SLFN11 expression in response to interferon signaling (Fig. 6B). Glioblastomas with higher expression of SLFN11 had gene set enrichment for activation of immune response, complement activation, defense response to virus, and response to interferon alpha compared with tumors with lower SLFN11 expression (Fig. 6C). Enrichment of these immune pathways likely reflects the constitutive activation of autocrine interferon signaling, which might facilitate immune escape for tumors (51).
To more directly test whether interferon signaling induces SLFN11 gene expression, we treated GSCs with IFNα and measured the expression of SLFN11 and type 1 immune response genes OAS1, ISG20, and IFITM 8 hours after treatment. All four GSCs increased expression of OAS1, ISG20, and IFITM following IFNα treatment, although the level of induction in the GSC 839 was modest (2–3-fold) compared with the other GSCs (Fig. 6D). IFNα treatment also increased SLFN11 gene expression in MNK1 and 839 GSCs, which have high baseline expression of SLFN11 and no ASE. However, IFNα treatment did not increase expression of SLFN11 in 2012 and 1552 GSCs, both of which showed ASE of SLFN11 and high promoter methylation (>25% βpromoter). These results suggest that DNA methylation of the SLFN11 promoter blocks upregulation of SLFN11 by interferon signaling (Fig. 6D).
To determine if GSCs with promoter methylation and low expression of SLFN11 were more susceptible to killing by Zika, we treated the same four GSCs with Zika or saline control (Fig. 6E; Supplementary Fig. S4). GSCs with promoter methylation of SLFN11 (2012 and 1552) displayed decreased viability following Zika infection compared with GSCs without promoter methylation (839 and MNK1). To directly test whether SLFN11 expression affects susceptibility to Zika, we performed SLFN11 OE in 2012 GSCs (which have low baseline expression of SLFN11) and SLFN11 KD in the 839 and MNK1 GSCs (which have high baseline expression of SLFN11). SLFN11 OE (Fig. 5C) decreased susceptibility to Zika (Fig. 6F), whereas KD (Fig. 5A) increased susceptibility to Zika (Fig. 6G).
Although each Zika experiment was performed with only a single biological replicate per GSC and more biological replicates would be desirable, we note that experiments were performed in multiple GSCs, and each experiment consisted of 3–10 technical replicates. Furthermore, our conclusions are supported by multiple lines of evidence. First, Zika infection and control experiments were performed on two GSCs with high SLFN11 expression (MNK1 and 839) and two GSCs with low SLFN11 expression (2012 and 1552). As expected, both of the low expression cell lines had lower viability in the presence of Zika (Fig. 6E). Second, overexpression of SLFN11 in the 2012 GSC (with low baseline SLFN11 expression) increased resistance to Zika (Fig. 6F). Third, knockdown of SLFN11 in the MNK1 GSC (with high baseline SLFN11 expression) increased susceptibility to Zika (Fig. 6G). In combination, these results provide compelling support for the hypothesis that SLFN11 expression is important for cell viability following Zika infection.
ASE analysis of tumor genomes is a new approach that can discover biomarkers and therapeutic targets by illuminating genes that are recurrently dysregulated. ASE analysis makes it possible to prioritize genes that are affected by cis-regulatory factors such as regulatory mutations, even when the precise identity of the mutational events that drive the transcriptional changes is unknown (52, 53). This differs from standard differential gene expression analysis, which identifies thousands of genes when comparing normal and tumor tissues. As gene expression is affected by the environment and by the activity of upstream factors, most of these differentially expressed genes are “passengers” that are downstream of the primary events that drive transcriptional changes. In contrast, by analyzing ASE we discover a relatively small number (118) of candidate disease genes that are recurrently dysregulated by cis-regulatory factors in GSCs, but not in normal tissues.
Many of the candidate genes identified by our ASE analysis have an established role in tumor biology. For example, IP6K2, a proapoptotic protein kinase showed ASE almost exclusively in GSCs. IP6K2 selectively binds to HSP90, which decreases its catalytic activity and inhibits apoptosis (36). Disruption of this interaction by cisplatin and novobiocin, chemotherapeutic compounds that bind to the C-terminus of HSP90, restores its catalytic function and promotes apoptosis (54). Furthermore, knockdown of IP6K2 in colorectal cancer cells has been demonstrated to selectively impair p53-mediated apoptosis, instead favoring cell-cycle arrest (55). These observations from previous studies suggest that IP6K2 may be an important tumor suppressor in glioblastoma.
NOTCH1 also exhibited ASE that was specific to GSCs. NOTCH1 regulates neural stem cell fate during neurogenesis and high expression of NOTCH1 has been reported in many high-grade gliomas (56–58). Notch1 signaling promotes invasion, self-renewal, and growth of GSCs (59, 60); NOTCH1-KD suppresses cell proliferation and induces apoptosis (61). Furthermore, inhibition of the Notch1 signaling pathway sensitized tumor cells to apoptosis induced by ionizing radiation, the death ligand TRAIL (tumor necrosis factor-related apoptosis-inducing ligand), or the Bcl-2/Bcl-XL inhibitor ABT-737 (62). These studies suggest that NOTCH1 may help maintain the stem cell–like behavior of GSCs and promote tumor progression. The multiple CREs correlated with NOTCH1 expression are potentially excellent targets for subsequent studies and cis-regulatory screens.
Recurrent ASE of SLFN11 is an important finding because this gene has recently emerged as a biomarker of drug sensitivity in cancer (63). We demonstrate that in GSCs, SLFN11 gene expression is associated with DNA methylation of its promoter and its expression is required for the antitumor activities of the DNA alkylating agent TMZ and the replication inhibitor olaparib. The current standard of care for glioblastoma patients is maximum safe surgical resection followed by concurrent TMZ and radiotherapy (64). Similar to MGMT promoter DNA methylation, SLFN11 promoter methylation may be a biomarker that predicts the efficacy of DNA damaging agents, such as TMZ and olaparib. In addition, promoter methylation and reduced expression of SLFN11 may reflect evolution of resistance to TMZ within tumor cells. In GSCs, SLFN11 gene expression was regulated by type 1 interferons and frequently upregulated in high-grade tumors, alongside constitutive activation of autocrine interferon signaling that facilitates immune evasion of GBM cancer cells (51). However, in the presence of promoter CpG methylation, SLFN11 is unresponsive to IFNα cytokine treatment, rendering GSCs vulnerable to killing by oncolytic viruses, such as Zika (47). Thus, tumors refractory to DNA damaging agents may be more amenable to treatment with genetically modified viruses. As GSCs with low SLFN11 expression are susceptible to Zika, but resistant to chemotherapy and vice versa, the combination of oncolytic viruses and chemotherapy may be a powerful treatment approach (Fig. 7).
One limitation of our study is that it does not incorporate the role of the tumor microenvironment in the model systems. Future efforts could investigate interactions between immune cells and tumor cells to elucidate how SLFN11 expression in tumor cells affects tumor maintenance and the antitumor immune responses.
J.N. Rich reports grants and personal fees from Synchronicity Pharma outside the submitted work. G. McVicker reports grants from NIH/NCI, Padres Pedal the Cause/RADY, and NIH/NHGRI during the conduct of the study. No disclosures were reported by the other authors.
A. Sen: Conceptualization, software, investigation, visualization, methodology, writing–original draft. B.C. Prager: Conceptualization, investigation, visualization, methodology, writing–original draft. C. Zhong: Validation, investigation. D. Park: Investigation. Z. Zhu: Investigation. R.C. Gimple: Investigation. Q. Wu: Investigation. J.A. Bernatchez: Investigation. S. Beck: Investigation. A.E. Clark: Investigation. J.L. Siqueira-Neto: Resources, supervision. J.N. Rich: Conceptualization, supervision, funding acquisition, project administration, writing–review and editing. G. McVicker: Conceptualization, supervision, funding acquisition, project administration, writing–review and editing.
This research was supported by the NCI funded Salk Institute Cancer Center (NIH/NCI CCSG: P30 014195) by a grant from Padres Pedal the Cause/RADY #PTC2019 (G. McVicker), by a Pioneer Fund Postdoctoral Scholar Award (A. Sen), and the following NIH grants and fellowships: CA217066 (B.C. Prager), CA217065 (R.C. Gimple), HG011315 (G. McVicker), CA197718, CA238662, and NS103434 (J.N. Rich). G. McVicker was supported by the Frederick B. Rentschler Developmental Chair. Figures 1A and 7 were created with BioRender.com.
The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.
Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).