Distinctive epigenomes characterize glioma stem cells and their response to differentiation cues
© The Author(s). 2018
Received: 4 December 2017
Accepted: 12 March 2018
Published: 27 March 2018
Glioma stem cells (GSCs) are a subpopulation of stem-like cells that contribute to glioblastoma (GBM) aggressiveness, recurrence, and resistance to radiation and chemotherapy. Therapeutically targeting the GSC population may improve patient survival, but unique vulnerabilities need to be identified.
We isolate GSCs from well-characterized GBM patient-derived xenografts (PDX), characterize their stemness properties using immunofluorescence staining, profile their epigenome including 5mC, 5hmC, 5fC/5caC, and two enhancer marks, and define their transcriptome. Fetal brain-derived neural stem/progenitor cells are used as a comparison to define potential unique and common molecular features between these different brain-derived cells with stem properties. Our integrative study reveals that abnormal expression of ten-eleven-translocation (TET) family members correlates with global levels of 5mC and 5fC/5caC and may be responsible for the distinct levels of these marks between glioma and neural stem cells. Heterogenous transcriptome and epigenome signatures among GSCs converge on several genes and pathways, including DNA damage response and cell proliferation, which are highly correlated with TET expression. Distinct enhancer landscapes are also strongly associated with differential gene regulation between glioma and neural stem cells; they exhibit unique co-localization patterns with DNA epigenetic mark switching events. Upon differentiation, glioma and neural stem cells exhibit distinct responses with regard to TET expression and DNA mark changes in the genome and GSCs fail to properly remodel their epigenome.
Our integrative epigenomic and transcriptomic characterization reveals fundamentally distinct yet potentially targetable biologic features of GSCs that result from their distinct epigenomic landscapes.
Glioblastoma (GBM) is the most common and malignant glioma subtype in adults . Average survival is only about 15 months, making GBM one of the most aggressive cancers. It has been hypothesized that distinct subpopulations with a survival advantage and enhanced resistance to therapy are contained within the bulk tumor, yet this property is not readily unveiled by gross histopathologic examination. Genetic and transcriptomic studies have subtyped GBM based on expression profiling into classical, mesenchymal, proneural, and neural, with each subtype characterized by distinct amplification/loss/mutation of genes such as EGFR, NF1, CDKN2A, and PTEN . Single-cell RNA sequencing (RNA-seq) revealed that an individual tumor contains a spectrum of GBM subtypes, suggesting that intratumoral heterogeneity is extensive [3, 4].
Theories underlying tumor evolution support the marked heterogeneity observed within individual gliomas. The cancer stem cell (CSC) theory postulates existence of a subpopulation of tumor cells residing at the apex of the hierarchy, propagating tumor formation in a hierarchical manner. CSCs are characterized by an ability to self-renew and differentiate, contributing to the heterogeneity and complexity of tumors. CSCs resemble normal stem cells in a number of properties, including the ability to form spheres on non-adherent culture surfaces in serum-free media . Relative quiescence coupled with low levels of apoptosis and slow cell cycling contribute to CSC resistance to chemotherapy, while their asymmetric division gives rise to poorly differentiated daughter cells that facilitate tumor recurrence [6, 7]. Oncogenic mutations occurring in normal stem cells could contribute to their malignant transformation into cancer stem cells. Early studies showed that manipulating the ARF/p53 pathway in neural stem/progenitor cells resulted in high-grade glioma [8, 9]. Glioma stem cells (GSCs) identified within bulk GBM tumors might therefore share biologic similarities with normal neural stem cells, but also possess distinct genetic and epigenetic alterations that underpin their malignant growth potential. Elucidating such differences is key to improving therapeutic targeting, efficiency, and specificity; however, such targetable epigenetic and transcriptomic differences between NSCs and GSCs remain largely unknown.
GSCs acquire both genetic and epigenetic mutations . Epigenetic changes, like genetic changes, act as driver events in transformation or collude with genetic events to drive transformation. In contrast to genetic alterations, epigenetic changes are, in principle, reversible and therefore represent attractive therapeutic targets. DNA methylation (5mC, mediated by the DNA methyltransferases DNMT1, 3A, and 3B) and DNA hydroxymethylation (5hmC, mediated by the ten-eleven translocation TET1, 2, 3 family) are extensively disrupted in GBM. The TET protein family is responsible for producing 5hmC, 5-formylcytosine (5fC), and 5-carboxylcytosine (5caC); and like the DNMTs, their expression is tightly regulated during development. Tumor cells, including GBM, are in general depleted for both 5mC and 5hmC, accompanied by reduced TET expression [11, 12]. GBM patients with G-CIMP (glioma-CpG island hypermethylator phenotype resulting from IDH1/2 mutation), or overall elevated 5mC and/or 5hmC levels, exhibit better clinical outcome [13, 14]. Although DNA methylation and transcriptional alterations have been studied extensively in GBM, details of the epigenetic reprogramming events that contribute to and/or define glioma stem cells, especially in terms of the TET-regulated DNA marks and their association with enhancer function, remain poorly characterized. In the current study, we aimed to define epigenetic abnormalities linked to GSC stemness and their regulation in response to differentiation cues relative to those characteristics of neural stem cells, to identify features unique to GSCs that may serve as novel therapeutic targets.
Characterization and global assessment of DNA epigenetic modifications in glioma stem cells
Expression of the DNA modification machinery in GSCs is variable and distinct from that in NSCs. TET1 transcription is overall downregulated, whereas TET2 and TDG are upregulated in the majority of GSCs; TET3 is the most variably expressed (Fig. 1d). Global quantification of DNA marks using mass spectrometry (MS)  show that NSCs overall possess higher amounts of 5mC and 5hmC, but lower levels of 5fC, compared to GSCs (Additional file 2: Figure S1B). Their levels were further compared to those of established serum cultured glioma cell lines T98G and A172, and normal brain tissue. Notably, the cancer cell lines exhibit lower levels of 5mC compared to normal tissue (Additional file 2: Figure S1C top), but much reduced levels of 5hmC and 5fC (Additional file 2: Figure S1C middle and bottom). 5caC levels are not shown because our LC-MS/MS method failed to capture the 5caC-containing fractions from all samples analyzed. In keeping with these findings, PDX tissue derived from GBM6, GBM64, and GBM84 exhibit reduced immunohistochemical (IHC) staining intensity in the nucleus for both 5mC and 5hmC, relative to surrounding normal tissue (Fig. 1e). We correlated transcription level of the TETs to the level of each DNA mark (Additional file 2: Figure S1D) and observed that TET2 and TDG expression are significantly negatively correlated with 5mC level, while TET2 is positively correlated with 5fC (Fig. 1f, Additional file 2: Figure S1D), suggesting that TET2 expression contributes to the increased level of 5fC in GSCs.
Transcriptome profiling reveals novel gene networks and GSC-specific targets
We identified several highly conserved downregulation events for microRNAs (miRNAs) in GSCs, including hsa-mir-31  and hsa-mir-210  (conserved in 6–9 GSCs), that have been reported to be downregulated in GBM, and seven poorly characterized miRNAs deregulated in GBM (conserved in > 9 GSCs) (Additional file 1: Table S5, validation of select miRNAs in Additional file 2: Figure S2C). GSC-downregulated miRNAs are potentially involved in processes modulating nucleosome assembly, chromatin silencing, epigenetic regulation, and telomere organization through their target genes (Additional file 1: Table S5, Additional file 2: Figure S2E). Coordination of these deregulated miRNAs may together lead to transcriptional deregulation in cancer via a large network (Fig. 2e), in which KAT6A and TET3 may also participate. Because we could not directly derive prognostic value of these miRNAs, we examined expression of genes predicted to be targets of the deregulated miRNAs for links to patient survival. A total of 193 and 543 genes are significantly associated with better and worse survival, respectively, leveraging clinical information and RNA-seq data of TCGA GBM samples (Additional file 1: Table S6 for complete list). We identified 20 GSC-upregulated miRNA-targeted genes associated with significantly reduced 5-year survival in GBM patients, including previously identified glioma-promoting genes HMGA2 , MYO1C , RGS4 , and several novel targets, such as HPCAL1, IMPDH1, and YWHAG (Additional file 1: Table S5, orange-colored genes, Additional file 2: Figure S2F). Therefore, deregulation of miRNAs in GSCs may contribute to poor outcome in part through modulating expression of these miRNA-targeted loci.
Ectopic gain and loss of enhancer marks in GSCs contribute to deregulated gene expression
Gene clusters surrounding (< 50 kb) GSC-specific active enhancers (common to at least six GSCs) are primarily linked to DNA damage responses, p53 signaling, and angiogenesis (Fig. 3d top); genes near NSC-specific active enhancers (lost in at least six GSCs) are enriched for stem-cell differentiation, apoptosis, and epigenetic gene regulation (Fig. 3d bottom, Additional file 1: Table S7). We next investigated relationships between ectopic enhancers and gene transcription. Many more RNA-seq reads map to enhancer-gained regions in GSCs (common to at least six GSCs); fewer reads map to enhancer-lost regions (Fig. 3e top and bottom, cluster 1). Meanwhile, downregulation events are proportionally greater than upregulation events for genes linked to regions that lose enhancer marks in GSCs (common to at least six GSCs) and vice versa (Additional file 2: Figure S3E). These findings suggest that ectopic gain/loss of active enhancers is strongly linked to gene deregulation events in GSCs. To confirm our findings, several genomic regions positive for H3K4me1 and/or H3K27ac modifications were validated using locus-specific ChIP-qPCR. Results are consistent with ChIP-seq signals for both enhancer marks at these loci (Additional file 2: Figure S3F and G).
HOX loci are upregulated in primary tumor-isolated GSCs despite the presence of promoter hypermethylation . Our data suggest that the HOX clusters acquire both H3K4me1 and H3K27ac in GSCs, which may explain the elevated HOX gene expression despite the presence of promoter methylation (Fig. 3f, Additional file 2: Figure S3H). Loss of flanking active enhancers is associated with gene repression of the PLEKHA5, PLEKHH2, and ST3GAL5 loci (Additional file 2: Figure S3H). PLEKHA5 facilitates cell migration and invasion and is linked to brain tumor metastasis via the PI3K-AKT pathway [24, 25], whereas ST3GAL5 encodes a sialyltransferase that catalyzes the formation of ganglioside GM3, a suppressor of EGFR/PI3K-AKT signaling and cell proliferation, and a inducer of apoptosis [26, 27]. Both gene clusters have known functions related to cancer cell properties, although their specific role in glioma stem cells has not been previously reported.
Genome-wide single-base resolution mapping of DNA epigenetic marks in GSCs and NSCs
We mapped DNA epigenetic modifications in greater detail by specifically examining 5mC, 5hmC, and 5fC/5caC, using modified-reduced representation bisulfite sequencing (RRBS) methods, including TET-assist bisulfite sequencing (TAB-RRBS) for 5hmC  and M.SssI-assisted bisulfite sequencing (MAB-RRBS) for 5fC + 5caC collectively . Modification level for each mark is stratified into low, medium, and high tiers. Migration from the high to the medium tier, rather than complete loss, is observed for 5mC in GSCs compared with NSCs (Fig. 4b, left). In contrast, 5hmC undergoes substantial loss in GSCs, as evidenced by the dominant tier shift from high/medium to low (Fig. 4b, middle). Interestingly, GSCs accumulate medium tier 5fC/5caC sites at the expense of lowly modified CpGs (Fig. 4b, right). We validated our RRBS results at several genomic regions in GSCs and NSCs using DNA immunoprecipitation with antibodies specific for each of the four DNA marks, coupled with real-time PCR (DIP-qPCR). DIP-qPCR results at these loci correlate well with RRBS as shown in the browser views for each DNA mark (Additional file 2: Figure S4C–E). More importantly, DIP-qPCR differentiates 5fC from 5caC at these loci (which our MAB-RRBS protocol does not), unveiling region- and cell line-specific preferences for 5fC and/or 5caC occupancy (Additional file 2: Figure S4C).
The genomic localization of highly-modified CpG sites is modification- and cell type-dependent. In the NSC profile, 5fC/5caC exhibits greater occupancy at promoter TSS1500 regions (5fC/5caC vs 5mC/5hmC: 24% vs 18%) and CpG islands (CGI) (9% vs 4%) (Additional file 2: Figure S4F). However, localization of highly modified sites differs in GSCs, with CGIs becoming more enriched for all three DNA marks at the expense of open sea for 5mC (GSC vs NSC: 19% vs 13%) and 5fC/5caC (28% vs 21%), and shore regions for 5hmC (CGI: 10% vs 12%, shore: 15% vs 13%) (Fig. 4c). We next compared the entire dataset for each mark between cell types and observed global loss of 5mC and 5hmC accompanied by an overall gain of 5fC/5caC in GSC. Specifically, loss of 5mC and 5hmC occurs mainly at intergenic and CpG-poor regions; gain of 5mC occurs more frequently at promoters, gene bodies, and CGIs (Fig. 4d left and middle). In contrast, gain/loss of 5fC/5caC lacks feature-specificity, evidenced by the even distribution across the genome (Fig. 4d right).
We then examined the relationship between each DNA mark, and between NSCs and GSCs. Substantial loss of both 5mC and 5hmC is observed simultaneously in GSC compared to NSC (n = 4658), with a small population of sites losing 5mC to 5hmC (n = 651) or losing 5hmC to 5mC (n = 348) (Fig. 4e, left). While GSCs have globally lower 5hmC levels, a considerable number of these sites meanwhile gain 5fC/5caC (n = 1392), followed by 884 sites that lose both 5hmC and 5fC/5caC (Fig. 4e, middle). Interestingly, we observed that hypo-5mC events are largely accompanied by increases in 5fC/5caC (n = 3962), followed by 2190 sites that lose both 5mC and 5fC/5caC (Fig. 4e, right). These results taken together suggest that GSCs replace 5mC and 5hmC extensively with 5fC/5caC in the background of a global loss of all DNA marks.
Distribution of 5hmC and 5fC/5caC reveals novel motif binding factors
Binding factors with cell-type specificity could lead to differential functional outcomes between NSC and GSC. We identified 16 GSC-specific and 130 NSC-specific motifs for 5mC, and 23 GSC-specific and 123 NSC-specific motifs for 5hmC (Additional file 1: Table S8). While binding of many of these factors has not previously been shown to be modulated by DNA marks, the expression of several of these factors correlates with levels of their associated DNA mark, and with GBM patient survival. For example, the NKX3–1 binding motif (CG-free) is unique at sites of GSC-specific 5mC and 5hmC, its expression ten times higher in GSCs than in NSCs (FPKM: 1.17 vs 0.12), and it is linked inversely to patient survival (Additional file 2: Figure S5B, top). The binding motif for KLF14 (CG-containing) is enriched only at 5mC-containing sites in NSC, and its transcription is higher in NSC relative to GSCs (Additional file 1: Table S8) and inversely related to patient survival (Additional file 2: Figure S5B, bottom).
Integration of DNA epigenetic modifications with the transcriptome and enhancer histone marks
Active enhancers in general are depleted of 5mC but enriched with 5hmC ; much less is known about 5fC/5caC at these regions. We examined co-localization of each type of highly modified DNA mark region with enhancer histone marks in NSC, compared with an expected pattern based on chance by referencing global distributions of all CpG sites derived from RRBS datasets. We observe that both 5mC and 5hmC are most enriched at H3K4me1-marked regions (1.9-fold expected for 5mC, 2.2-fold expected for 5hmC), with 5mC preferentially enriched at H3K4me1-marked poised enhancers (5mC vs 5hmC: 29% vs 17%) and 5hmC enriched at H3K4me1 and H3K27ac co-marked active enhancers (5mC vs 5hmC: 45% vs 60%). In contrast, 5fC/5caC, while being consistently enriched at H3K4me1-marked regions, shows a onefold increased enrichment at H3K27ac marked regions relative to that of 5mC and 5hmC (5fC/5caC vs 5mC vs 5hmC: 50% vs 26% vs 23%) (Fig. 6b).
We defined 16 chromatin states by integrating the histone marks from two NSCs and 12 GSCs into a ChromHMM model and then grouped them into five scenarios: GSC loss, NSC-GSC conserved, NSC-GSC transition, GSC gain, and NSC-GSC switch (Fig. 6c, left column). When we examined the relationship between these groups and tissue-specific enhancers defined by H3K4me1, we observe that GSC-lost active enhancers co-localize with a subset of fetal brain-specific enhancers that also exhibit enrichment for other tissue-specific enhancers (states 16, 15, and 2). Meanwhile, GSC-gained enhancers overlap with other tissue-specific enhancers at multiple states, including transitioning state 5 and 12, and switching state 9, in which only NSC possess H3K27ac or H3K4me1, but there is a switch to poised or active enhancer status in GSC. GSC-unique gained enhancer groups (states 8, 6, and 13) show minimal co-localization with other tissue-specific enhancers (Fig. 6c middle). Such enhancer mark-switching patterns suggest that GSCs lack the tissue-specificity of NSCs in part due to loss of tissue-specific footprints established by brain/neural-specific enhancers and gain of other tissue-specific and/or random enhancer marks in the genome.
We demarcated distinct DNA mark switching events between GSCs and NSCs and investigated their co-localization with the most prevalent ChromHMM state changes (Additional file 2: Figure S6A, Additional file 1: Table S9). Several states having a unique combination of enhancer and DNA mark changes were analyzed for their potential impact on nearby genes. A collection of hypermethylation events, regardless of the change in 5hmC and 5fC/5caC, are a major driving force clustering a subgroup of states with exclusive loss of H3K27ac in GSC relative to NSC (state 10). A total of 283 genes functionally relevant to nervous system development and cell adhesion are adjacent to such regions (< 10 kb) (Additional file 1: Table S10, column C), with downregulated expression in GSC (n = 41) being the dominant effect (Fig. 6d, top). We compared this result to TCGA GBM data and found 24 genes consistently repressed in primary tumors, including RANBP17, an activator of p21  and an indicator of poor survival (Additional file 2: Figure S6B). An additional 17 downregulated genes are unique to GSC lines, including FOXE1, a developmental transcription factor whose promoter and adjacent enhancer are specifically targeted for DNA hypermethylation in GSC lines (but not bulk tumor based on comparisons with TCGA data). Low FOXE1 expression also correlates with reduced survival of GBM patients (Additional file 2: Figure S6C). Consistent with their hypermethylation-mediated downregulation in GSC lines, ectopic re-expression of RANBP17 and FOXE1 leads to attenuated cell viability in the GSC6 and GSC14 lines that do not express these genes from the endogenous loci (Additional file 2: Figure S6D, E), suggesting they represent growth suppressors targeted via epigenetic mechanisms to promote GSC growth or survival.
Exclusive loss of 5mC and gain of 5hmC in GSC is enriched at states 12 and 8 and correlates with a set of GSC-specific H3K4me1 peaks. Gene deregulation in GSCs is observed near state 12 (loss of 5mC), including 37 upregulation and 105 downregulation events, among which 14 and 36 genes exhibit the same form of deregulation in TCGA GBM samples, respectively (Additional file 1: Table S10, column G). Upregulated genes, such as EGFR and MYO1C, are associated with glioma cell proliferation and migration, whereas downregulated genes such as MMP11 and SREBF1 are related to cell adhesion and lipid metabolism (Fig. 6d, bottom). Regions exclusively gaining 5hmC in state 12 are not linked to altered expression events, whereas 5hmC loss exhibits slightly higher enrichment at states 2, 3, and 16. This lack of enrichment specificity compared to altered 5mC events is consistent with the ubiquitous loss of 5hmC observed in GSCs.
Focusing on switching events among 5mC and 5fC/5caC marks between NSC to GSC reveals a novel correlation; loss of 5mC and gain of 5fC/5caC strongly occupies state 9, a collection of territories that lose H3K27ac to H3K4me1 in GSC (Fig. 6c). A total of 43 genes reside within 10 kb of these regions (Additional file 1: Table S10, column B), including those involved in cancer pathways (e.g. MMP2), angiogenesis (e.g. VEGFA), hypoxia response (e.g. SCAP), and neural development (e.g. KCNQ2), with downregulation of these loci being the dominant transcriptional effect (Additional file 1: Table S3). On the other hand, 5mC gain with 5fC/5caC loss events are enriched at state 14, regions that lose H3K4me1 to H3K27ac (Fig. 6c). Few genes are located near these regions and only a small subset are transcriptionally deregulated in GSC including ANK1, a gene known to maintain cell morphology and promote cancer cell growth . Interestingly, while ANK1 is upregulated in GSCs, it is downregulated overall in primary GBM based on TCGA data (Additional file 2: Figure S6F, top panel), suggesting that GBM intratumoral heterogeneity masks molecular characteristics in GSCs. Using survival information from these 151 patients, we observe that high expression of ANK1 in GBM is associated with poor survival (Additional file 2: Figure S6F, bottom panel).
Differentiation of neural and glioma stem cells is characterized by distinctive gene regulation events
In contrast to the phenotypic changes in NSCs, expression programming of GSC stem regulators  SOX-2, OLIG2, POU3F2, and SALL2, and lineage marks is heterogeneous under differentiating conditions. While OLIG2 is consistently repressed, responses of other loci vary among GSC lines in response to differentiation cues (Additional file 2: Figure S7C). Induction of GFAP occurs in a subset of GSCs, including GSC12, 39, 84, and 102; TUBB3 activation is observed in GSC6, 10, 12. Protein levels of these markers were quantified by IF with data summarized as the percent positive cells in Additional file 1: Table S1, bottom. DNMT1, DNMT3A, and DNMT3B are downregulated in several GSCs upon differentiation (Additional file 2: Figure S7D). TET2 expression is largely unchanged, whereas TET1, TET3, and TDG respond variably (Additional file 2: Figure S7E). We examined relationships between expression of stem/lineage markers and DNMT/TET expression during GSC differentiation, and observe three clusters of positive correlation, including DNMT3A correlating with OLIG2, SALL2, NESTIN, and TUBB3, and TET1/TET2 correlating with stem markers POU3F2 and SOX2, and lineage marker GFAP. Of note, several inversely correlated gene pairs were observed, including DNMT1, TET3, and TDG, whose expression is anti-correlated with DNMT3B, DNMT3L, and GFAP (Fig. 7c).
In addition to altered expression, dynamic changes in TET subcellular localization occur during GSC differentiation. IF staining for TET expression is consistent with our previous qRT-PCR data (Fig. 1d) and shows that TET1 is highly expressed in NSCs, while TET2 is highly expressed in GSCs (Additional file 2: Figure S7F). Moreover, TET1 and TET3 are more prominently localized to the nucleus in GSCs compared to NSCs. We co-stained for each TET and GFAP in GSC84, a GSC with robust astrocytic lineage differentiation (Fig. 7d). While all three TETs translocate from the nucleus to the cytoplasm and many cells display co-expression of cytoplasmic TET3 and GFAP, ~ 20% of the cells fail to induce GFAP expression under FBS-differentiation conditions (Fig. 7e). These cells express mainly nuclear TET3 and are devoid of GFAP (Fig. 7f, white boxed areas), suggesting that TET3 subcellular redistribution is important for lineage differentiation and the ongoing presence of nuclear TET3 inhibits GFAP expression.
Atypical reprogramming of the transcriptome and epigenome during glioma stem-cell differentiation
Differentiation induction led to extensive DNA epigenetic modification changes in a cell lineage-specific manner. Astrocytic and neuronal induction triggers an almost equal number of hyper- and hypomethylation events in NSC, whereas GSCs under FBS conditions trigger a more variable set of 5mC alterations (Additional file 2: Figure S8C, top). Because 5mC changes are heavily influenced by cell type, samples clustered by cell line rather than differentiation status (Additional file 2: Figure S8D). In contrast, differentiation induced a much greater number of 5hmC gains in both NSCs and GSCs (Additional file 2: Figure S8C, bottom), which in turn led to two distinct clusters formed by differentiation status rather than cell type (Additional file 2: Figure S8E). In addition, GSCs share very few common modification-change events with NSCs during differentiation, again showing that differential regulation of 5mC and 5hmC is highly specific to cell type (Additional file 2: Figure S8F).
Genomic distribution of differentiation-induced 5mC and 5hmC changes did not differ greatly by gene feature or CpG-density. 5mC and 5hmC changes occur more frequently in CpG-sparse and intergenic regions and less commonly in promoters and CpG islands (Fig. 8c, d). We next mapped these sites to the enhancer landscape in the undifferentiated state described earlier and identify distinct localization of DNA marks among GSC and NSC (Fig. 8e). For example, NSCs lose 5hmC at states 9 and 10, while GSCs lose 5hmC at state 16 where NSCs display minimal 5hmC changes. We examined whether differentiation-induced loss of 5hmC in NSC at H3K27ac-marked regions (states 9 and 10) affects expression of nearby genes (Fig. 8e box 1). Downregulation events are indeed identified for such genes including GFRA1, GRIN2D, NKX2–2, and PDGFB, whereas upregulation of CLDN11 occurs in differentiated NSC. This finding suggests that 5hmC loss during differentiation at previously active regions is associated with decreased expression. NKX2–2 expression decreases upon astrocyte induction , whereas CLDN11 is upregulated during neural induction . Our data, in keeping with a previous study, show that NKX2–2 is suppressed under astrocytic conditions, while CLDN11 is activated in neuronal conditions. GSCs under differentiation conditions, however, do not show comparable changes, indicating a differentiation defect likely resulting from failure to remodel their epigenome. Hypermethylation occurs specifically during GSC differentiation at state 16 (Fig. 8e, box 2), a collection of regions co-occupied by active enhancers in NSC but poised enhancers in GSC. Genes adjacent to such regions are involved in development and differentiation-related processes (Fig. 8f), thus an increase in 5mC at flanking enhancers may exert an inhibitory effect on differentiation. Indeed, we identified a subset of genes, including MBP (a glial marker) , COL18A1 (a neural retina gene) , SREBF1 (nerve fatty acid synthesis) , and WNT1 (neurogenesis gene)  that are upregulated in differentiation-induced NSC27, but not in any GSCs. Taken together, these data suggest that NSCs acquire 5mC/5hmC reprogramming events at regulatory regions of developmental genes during differentiation; this reprogramming process, however, is less specific in GSCs, which may account for their more varied response to differentiation cues and the ability of “differentiated” derivatives to continue to proliferate and contribute to the bulk tumor mass.
Co-expression analysis and differential promoter DNA modification levels link TET2 and TET3 to the DNA damage response
A co-expression analysis approach was taken to identify genes potentially regulated by each TET in both our GSC dataset and the TCGA GBM dataset. Using our RNA-seq data, a total of 878, 2229, and 771 genes are linked to TET1, TET2, and TET3 expression, respectively (Additional file 1: Table S14). Based on this analysis each TET relates to a unique set of genes, but TET2 and TET3 share more common genes while TET1 is linked to a distinct group (Fig. 9b). Genes co-expressed with TET1 or TET3 are not significantly enriched for specific pathways (although several DNA damage repair and cell-cycle-regulating genes are in the TET3 list but do not reach the cutoff for statistical significance [not shown]). Genes positively correlating with TET2 expression are, however, strongly enriched in cellular processes including cell-cycle regulation, DNA damage repair, and telomere maintenance; TET2 negatively correlates with genes associated with the lysosome and cell adhesion processes (Fig. 9c, Additional file 1: Table S15). Next, we examined whether similar co-expression patterns (downloaded from cBioPortal) exist in bulk tumor and from a larger independent dataset using TCGA GBM data. Unlike GSCs, the TETs share a greater number of commonly co-expressed genes in primary GBM (Fig. 9d). Again, genes involved in lysosome regulation, metabolic pathways, and oxidative phosphorylation are negatively correlated with TET expression in general, while biological processes including DNA damage repair, cell proliferation, and pluripotency of stem cells are enriched for positively correlated gene sets for all three TETs (Additional file 1: Table S16). These co-expression patterns based on primary tumor data are well-preserved in GSCs for TET2, but not for TET1 and TET3, strongly suggesting that TET2 is involved in regulating cell-cycle control and DNA damage repair, functions that have been previously linked to TET2 and 5hmC .
The links we identified between TET2 and expression of genes regulating aspects of DNA repair motivated us to examine this link functionally. We monitored the viability of a representative GSC line with high TET2 expression (GSC6), a low TET2-expressing line (GSC84), and NSC30 (TET2 levels comparable to GSC84) after exposure to bleomycin, a chemotherapeutic agent that induces DNA double-strand breaks and mimics the effects of ionizing radiation . NSC30 exhibits reduced proliferation at low bleomycin concentration. GSC84, whose TET2 expression is close to that of NSC30, also shows significantly attenuated proliferation at low bleomycin concentration. Proliferation of GSC6 (highest TET2 levels), however, is least affected by bleomycin (Fig. 9e). To directly investigate the contribution of TET2 to the response of GSC6 toward bleomycin-induced DNA damage, we performed the same assay under TET2- and TET3-small interfering RNA (siRNA) knockdown conditions (Additional file 2: Figure S9D, TET3 is included because it too exhibits a modest expression correlation with DNA repair genes). Results show that proliferation of TET2-depleted GSC6 is significantly reduced by bleomycin treatment relative to a no-target control siRNA transfection (P value = 0.009), while TET3 depletion also decreases cell proliferation, but to a lesser extent (P value = 0.015) (Fig. 9f). That suppression of TET2 and TET3 in GSC6 leads to attenuated proliferation under conditions of DNA damage suggests that they play a direct role in modulating the DNA damage response process and may contribute to the chemo−/radio-resistance of glioma stem cells.
Non-CpG methylation exhibits unique patterns in GSCs and in response to differentiation
Upon differentiation there is an overall increase in non-CpG methylation (delta > 0.1), with neuronal induction in NSCs resulting in the greatest number of hypermethylated non-CpG sites (Fig. 10c). Meanwhile, hypermethylation occurs more frequently within CpA, while the majority of hypomethylation events occur in CpC and CpA context (Fig. 10d). Correlation analysis reveals a similar pattern to that of CpG methylation, in which differentiation induced non-CpG methylation changes are highly specific to cell type and vary by genomic feature (Fig. 10e, Additional file 2: Figure S10E). Differentiation-induced non-CpG hypermethylation is also lineage-specific as neuronal differentiation triggers a substantial and unique set of gains (n = 3910) in non-CpG methylation in NSC (Fig. 10f left), while hypomethylation events are more common between two lineages (Fig. 10f right). In contrast, differential non-CpG methylation is largely unique to each GSC, especially for hypomethylation events (Fig. 10g), suggesting a distinct differentiation path between NSC and GSC.
Glioma stem cells exist as a minor subpopulation within the bulk tumor, but actively contribute to overall radiation and chemo-resistance, and tumor recurrence . We referenced fetal brain-derived neural stem cells, since they are a possible cell of origin for GSCs, and profiled genome-wide DNA epigenetic modifications, enhancer marks, and transcriptomes in multiple PDX-derived GSCs. We performed an integrated analysis and compared epigenetic signatures before and after differentiation in each cell type to pinpoint cancer stem-cell-specific epigenetic abnormalities. Our data reveal several novel findings: (1) deregulation of TET expression correlates with global loss of 5mC and 5hmC, and gain of 5fC/5caC in GSCs, switching patterns of which display a distinct co-localization with enhancer mark changes; (2) despite heterogeneity within GSCs at both the transcriptional and enhancer levels, these cells exhibit substantial overlap in the deregulated genes/pathways that are highly associated with the differentially modified enhancer landscape between GSC and NSC; (3) TET expression responds differently during differentiation in GSCs compared to NSCs, with subcellular localization of TET3 associated with differentiation proficiency; (4) differentiation-induced reprogramming of 5mC and 5hmC localizes to distinct enhancer regions and contributes to the regulation of developmental genes; (5) integrative analysis identifies novel epigenetically deregulated loci relevant to GBM patient survival and GSC growth; and (6) high-level TET2 expression in GSCs contributes to their chemo-resistance, likely through modulation of the DNA damage response and repair system.
A primary goal of our study is to identify unique epigenomic features and associated gene networks in glioma stem cells that could be targeted clinically without harming neighboring normal cells. Global hypo-5mC/5hmC and promoter-specific hyper-5mC events are frequent in cancer. A potential division of labor exists within the TET protein family, with TET2 more efficient at generating 5fC and 5caC  and playing a greater role at enhancers . For the first time, we show a global elevation in 5fC/5caC levels in GSCs, which is potentially attributable to elevated TET2 expression. Further support for this notion comes from our finding that TET2 expression strongly correlates with 5fC/5caC quantity and higher levels of TET2 transcription and 5fC/5caC in GSCs. Through global correlation analysis, we also discovered that TET2 expression is associated with a more robust response to DNA damage in GSCs. Previous studies demonstrated that GSCs exhibit enhanced ability to repair DNA damage through activation of checkpoint proteins such as ATM, RAD17, and CHK1/2 after exposure to radiation . This survival advantage is eliminated by inhibition of CHK1/2; however, normal stem cells in this environment may suffer by acquiring undesired mutations due to the low-selectivity of CHK1/2 inhibition. Our results reveal that GSCs express DNA repair genes at a higher level than NSCs, which may contribute to their resistance to radiation- or chemotherapy-induced DNA damage. Meanwhile, because TET2 expression is higher in GSCs than NSCs, and TET2 expression linearly correlates with expression of repair genes, TET2 may modulate expression of this group of genes and thus serve as a potential therapeutic target. Given that TET2 is lowly expressed in neural stem cells and is further repressed during NSC differentiation, TET2 inhibition might have less impact on NSCs than on GSCs. In support of this, depletion of TET2 using siRNA knockdown sensitizes a GSC line to bleomycin treatment, suggesting that TET2 indeed functions to promote DNA repair or survival after DNA damage. A recent study showed that 5hmC localizes to DNA damage foci and facilitates DNA repair, a process requiring TET2 expression . Our study suggests that TET2 and 5fC/5caC might be specifically involved in this process in GSCs. In addition, TET2 may also be responsible for generating promoter 5fC/5caC at cell cycle and DNA damage response-related genes whose expression is higher in GSC than NSC and positively correlated with TET2 levels, suggesting a possible mechanism through which TET2 modulates expression of these loci.
Another potential approach to eliminate or repress GSCs is to trigger an effective and lasting differentiation process. Previous studies revealed that in vitro induction of differentiation by growth factor withdrawal and addition of serum/bone morphogenetic proteins (BMPs) only triggers a process of pseudo-differentiation. These “differentiated” cells are vulnerable to cell-cycle re-entry and do not properly remodel their DNA methylation patterns . Based on our study, we observe that the DNA methylation/demethylation machineries are not regulated in the same way in GSCs compared to NSCs, even though stem-cell markers, cell-cycle genes, and DNA repair pathways are inhibited upon differentiation induction as expected. This phenomenon might be directly attributable to the distinctive patterns of methylation/hydroxymethylation reprogramming in NSC and GSC, in which gain and loss of region-specific epigenetic marks is required for proper and unidirectional differentiation. More importantly, in addition to the well-established role of the TETs in brain development, we observed a specific abnormality in TET3 expression in GSC upon differentiation. Although TET expression in GSCs did not decrease in response to differentiation induction, the localization pattern of TET shifted from the nucleus to the cytoplasm and was accompanied by increased GFAP expression as shown in GSC84, one of the most efficiently differentiated GSCs. An unbiased correlation analysis revealed that TET3 expression is inversely correlated with GFAP expression. IF analysis indeed showed that when TET3 fails to translocate to the cytoplasm, GFAP expression is also attenuated. This irregular TET expression change and/or subcellular distribution under differentiation conditions might contribute to the inability of GSCs to effectively remodel their epigenome from a stem to a lineage committed and differentiated status, and thus account in part for the pseudo-differentiation phenotype.
One possible origin for GSCs is through transformation of neural stem/progenitor cells. As such, GSCs may possess attributes that contribute to normal stem-cell maintenance, while also possessing unique epigenetic characteristics associated with or driving their transformation. A previous study showed that neural progenitors committed to a specific lineage give rise to distinct glioma subtypes, even in response to the same driver mutation . Given that malignant transformation of neural progenitors may occur at different developmental stages (lineage restricted or not) and originate from different regions of the brain, it is possible that GSCs derived from this process preserve distinct epigenetic signatures of their ancestry, which in part contributes to heterogeneity within the GSC population. An earlier study reported that reprogramming of the DNA methylome in GSCs to a non-neural lineage reactivates expression of tumor suppressors and inhibits tumorigenesis . This finding suggests that even though GSCs originate through different mechanisms and/or cell types reflected by their genetic and epigenetic heterogeneity, these events collectively play a critical role in maintaining their cancer stem-cell identity and likely converge on common biological characteristics, as shown in our results and those of others . Meanwhile, we observe heterogeneity among GSCs at multiple levels, including transcription, DNA modifications, and enhancer marks. This heterogeneity is also partially reflected at the level of TET expression among the GSCs. Recently, Jin et al. reported that heterogeneity among GSCs gives rise to distinct sensitivities to different epigenetic targeting mechanism . Our observations related to TET2 and its relationship to DNA damage repair also support this conclusion, in that high TET2-expressing GSCs might be clinically more resistant to drug therapy and that a treatment targeting TET2 is likely to be effective in GSCs expressing TET2 at high level, instead of all GSCs, when combined with other cancer therapies.
The present study has identified several novel epigenetically and transcriptionally convergent events, pinpointed key relationships between the TET machinery and unique epigenomic features of GSCs, and revealed distinctive GSC-specific patterns of TET expression during differentiation. These GSC-specific features likely bestow this cell population with unique molecular profiles and epigenetic properties, which represent potential therapeutic targets for inhibiting or eliminating the GSC population.
Glioma and neural stem-cell culture
Tumor tissue from a well-characterized panel of GBM patient-derived xenografts (PDX) developed as part of the Mayo Clinic Brain SPORE, was mechanically dissociated and selected for the stem-cell population in stem-cell culture medium (ThermoFisher, StemPro NSC SFM A1050901) containing L-glutamine and penicillin-streptomycin using laminin (#L2020) coated flasks to reduce culture condition-introduced heterogeneity [15, 58, 59]. Fetal brain-derived neural stem cells, NSC23, NSC27, and NSC30, were obtained from Dr. Philip H. Schwartz at the Children’s Hospital at Orange County (CHOC) . Neural stem cells were expanded on Matrigel (Corning 356,230) or fibronectin-coated 6-well plates in stem-cell growth medium (GM) containing DMEM/F-12 plus Glutamax (Life Technologies 10,565–018), BIT9500 serum substitute (1×, Stem Cell Technologies 09500), heparin (20 μg/mL Sigma H3149), primocin (1×, InvivoGen ant-pm-1), human bFGF (20 ng/mL, Stemgent 03–0002), and human EGF (20 ng/mL, R&D Biosystems 236-EG) . Low passage cultures of cells (< 5 and typically < 3 passages) were used for all experiments and deep sequencing.
Differentiation of glioma and neural stem cells
Cells were seeded at low density (~ 0.2–0.5 × 10^6 per well of 6-well plate). Culture vessels were coated with poly-l-lysine followed by laminin to promote cell adherence. For neuronal lineage differentiation, growth factors bFGF and EGF were withdrawn for GSCs; neural differentiation medium containing Neurobasal medium (ThermoFisher, 21,103), B-27 supplement (ThermoFisher, 17,504), GlutaMax-I supplement (ThermoFisher, 35,050), dibutyryl cAMP (0.5 mM added on day 7 for three days) and primocin was used for NSCs. For astrocytic lineage differentiation, 10% FBS was added to the culture medium for GSCs; astrocyte differentiation medium containing DMEM (ThermoFisher, 11,995), N-2 supplement (ThermoFisher, 17,502), GlutaMax-I supplement, 10% FBS, and primocin was used for NSCs. Cells were incubated in differentiation media for 20–30 days.
Immunofluorescence, immunohistochemistry, and imaging
For IF, cells were seeded onto a fibronectin- or laminin-coated surface and fixed with 4% paraformaldehyde (Electron Microscopy Sciences, 15,710) and permeabilized with 0.1% Triton X-100, followed by 5% goat serum blocking before primary antibody incubation overnight at 4 °C. IHC staining for FFPE tissues from PDX (5-μm thickness) was performed at the Mayo Clinic Pathology Research Core. Slides were retrieved using Epitope Retrieval 2 (EDTA; Leica), incubated in Rodent Block M (Biocare) for 30 min, followed by 15 min incubation of the primary antibody diluted in Background Reducing Diluent (Dako). DAB and DAB buffer (1:19 mixture) and Schmidt hematoxylin were applied for visualization (Bond Polymer Refine Detection System). The dilution factor for each antibody is as follows, anti-Nestin (AB92391,1:200), anti-SOX2 (MAB2018,1:50), anti-GFAP (AB4648, 1:100), anti-CD44 (HPA005785, 1:100), anti-Tuj1 (801,201, 1:100), anti-Galactocerebroside (MAB342, 1:100), anti-5mC (Calbiochem; Clone 16233D3, 1: 200), and anti-5hmC (Active Motif, 39,769, 1:1500). Imaging was performed using an EVOS FL cell imaging system (Life Technologies). Positive cells were counted from at least three fields (each containing > 50 cells). The “stemness” was ranked based on the proportion of cells staining positive for stem and lineage markers.
Ectopic expression and knockdown, bleomycin treatment, and cell viability measurements
An infection complex composed of 20 μg of FOXE1 (VB170518-1071veg) or RANBP17 (VB170518-1058gjg) expression vector (Cyagen Biosciences Inc.), 10 μg packaging plasmid psPAX2, and 5 μg envelope plasmid pMD2.G were combined with 70 μL X-tremeGENE HP DNA Transfection Reagent in 2 mL opti-MEM media, which was then applied to HEK293T cells. Virus was collected at day 4 and day 5 from the supernatant and filtered through a 0.45-μm filter. Virus was concentrated at 26,000 rpm for 2 h at 4 °C then resuspended in 1 mL serum-free DMEM. For transduction, 2–3 × 105 cells per well of GSC6 and GSC14 were seeded into a 6-well plate, followed by addition of 50–100 μL concentrated virus and 6.5 μL polybrene (4 μg/mL) for 8–12 h. Fresh GSC media was then added and cells were harvested at the indicated timepoints for viability assay (Promega CellTiter 96 Cell Proliferation Assay), RNA isolation, and real-time PCR. For knockdown experiments, GSC6 cells were seeded at a density of 0.5 × 10^6 per well in 6-well plates. Twelve hours after seeding, cells were transfected with siRNA no-target control (siNTC) and siTET2/TET3 (GE Dharmacon, SMARTpool) at 50 nM and 10 uL PepMute™ siRNA Transfection Reagent (SignaGen Laboratories) per well. After 24 h, cells were transfected again with the same reagents at the same concentration to increase knockdown efficiency. For assessing the effect of siRNA knockdown, cells were collected for RNA isolation and real-time PCR after three days. For bleomycin treatment, transfected cells were re-seeded into 96-well plates at a density of 5000 cells per well and supplied with different concentrations of bleomycin for three days. An MTS assay (Promega, Madison, WI, USA) was then performed to assess viability.
RNA isolation, real-time PCR, and transcriptomic profiling by RNA-seq
RNA was extracted using Trizol reagent following the manufacturer’s protocol (Life Technologies). Reverse transcription and real-time PCR were performed as described . MiRNA was extracted with the mirVana miRNA Isolation Kit (Life Technologies #AM1561), and reverse transcribed and detected by real-time PCR using the miScript PCR system (Qiagen). Primer sequences are in Additional file 1: Table S17. Quality control of mRNA, library preparation, and paired-end sequencing were conducted at the University of Minnesota Genomics Center using an Illumina HiSeq 2500. Approximately 30 million reads were obtained for each sample (Additional file 1: Table S18). Raw reads were mapped to the human genome hg19 and transcripts assembled using TopHat. Quantification and pair-wise comparison of gene expression was conducted through Cufflinks and Cuffdiff packages, respectively. Fragments per kilobase of transcript per million mapped reads (FPKM) was generated to represent transcriptional level and for generating figures. A total of 57,883 unique transcripts were mapped, including 1416 miRNAs. Differentially expressed genes were defined using both fold-change and statistical q-value generated by the Cuffdiff package (false discovery rate [FDR]) (cutoff > 1 or < − 1).
Differential and co-expression analysis using RNA-seq data
To identify gene sets exhibiting co-expression with each TET we required a Pearson correlation P < 0.05 and log2-transformed expression fold change between the highest and lowest sample > 0.5. For TCGA GBM data, the co-expression gene list was acquired from cBioPortal using a Pearson correlation > 0.3 or < − 0.3. MiRNA targeting and network analysis was performed with miRDB (http://www.mirdb.org/miRDB/) and miRNet (http://www.mirnet.ca/) databases. Relationships between gene expression and patient survival were determined by the survival P value, calculated with a R package “Survival” under the “survdiff” function, whereas the odds ratio and its P value were calculated under the “coxph” function. Only genes with both P < 0.05 and expressed in more than half of all patient samples (n > 75) were considered to have prognostic value.
Enhancer mapping using chromatin immunoprecipitation sequencing (ChIP-seq)
In brief, 10 million cells were fixed with 1% formaldehyde. Fixed cells were pelleted and lysed with lysis buffer. Cells were subjected to MNase digestion (NEB, M0247S), followed by sonication using a Diagenode Twin sonicator. Chromatin input was taken from the supernatant and quantified with a Qubit 3.0 fluorimeter. Antibodies recognizing H3K4me1 (AB8895) and H3K27ac (AB4729) were added to immunoprecipitate enhancer-associated DNA fragments using GE magnetic sepharose protein G beads (#28–9670-70). The beads were washed with high-salt, low-salt, and Tris/LiCl buffer. The binding complex was eluted from the beads and reverse cross-linked. DNA was then treated with RNase and proteinase K and purified using a Qiagen mini-elute kit. Paired-end sequencing was performed on an Illumina HiSeq 2500 at the Medical Genome Facility at Mayo Clinic using the ThruPLEX DNA-seq kit (Rubicon Genomics #R400428). The coverage report can be found in Additional file 1: Table S19.
Stratifying cell type and tissue specific enhancers with ChIP-seq data
ChIP-seq peak calling was performed with the SICER package at the default settings and corrected for input. Differential enrichment was identified using SICER-diff at FDR < 0.01. Unique and common differentially enriched peaks across all glioma stem-cell lines vs neural stem-cell lines were identified with BEDOPS v2.4.20. The R package “ChIPseeker” was used to annotate ChIP peaks with general genomic features and compare feature distribution across different subsets of peaks. Tag density plots of transcript abundance centered at enhancer marks was plotted with deepTools2 software. Genes near enhancer peaks were mapped using Genomic Regions Enrichment of Annotations Tool (GREAT). Human tissue-specific enhancer regions marked with H3K4me1 were extracted from ten different tissue types using BEDOPS v2.4.20, including H1ES (GSM537679), fetal brain (GSM806934), gastric (GSM910574), pancreas (GSM910576), small intestine (GSM956019), lung (GSM1059443), thymus (GSM1059446), bladder (GSM1059450), liver (GSM1059451), and spleen (GSM1120351). Peaks called by SICER with FDR < 0.01 and enrichment score ≥ 30 were used for selecting regions uniquely marked in each tissue type. On average, 60,000 peaks from each tissue were assessed for overlap. When at least 60% of a single peak from one tissue overlapped with another, they were considered shared between those tissues. A total of 32,640 unique H3K4me1 peak regions were identified as brain-specific enhancers.
Genomic DNA isolation, quantification by mass spectrometry, and DNA mark profiling by (modified) reduced representation bisulfite sequencing (RRBS)
Genomic DNA was isolated using a serial phenol: chloroform extraction and was quantified for total cellular levels of 5mC, 5hmC, 5fC, and 5caC by LC-MS/MS as described . DNA from human ES, neural stem cells, human fetal brain, normal astrocytes, oligodendrocyte precursor cells, normal white/Gy matter-regions, and PDX-derived GSCs was assayed for methylation using the Infinium HumanMethylation450K BeadChip. DNA methylation of the primary PDX tumors was assayed on the Infinium MethylationEPIC BeadChip (only probes overlapping with the 450 K beadchip were used). For RRBS, genomic DNA was subjected to MspI digestion, size selection, and library preparation following a standard RRBS protocol. Paired-end sequencing was conducted on an Illumina HiSeq 2000 at the Medical Genome Facility at Mayo Clinic Rochester. Raw sequencing reads were aligned to human genome hg19 and analyzed through a Streamlined Analysis and Annotation Pipeline (SAAP-RRBS) . A methylation score was assigned based on the ratio of total methylated reads/ total sequenced reads. Only CpGs with > 10X coverage were used for downstream analysis. For TET-assist bisulfite sequencing (TAB-RRBS), MspI digested DNA was first treated with β-glucosyltransferase to protect all 5hmC sites from the subsequent Tet1-based oxidation reaction following the company’s standard protocol (WiseGene kit #K001). Spike-in 5hmC and 5mC controls were used for evaluating the efficiency of 5hmC protection and Tet1 oxidation. For the TAB kit we used, the 5hmC protection rate was nearly 100% and the conversion efficiency of recombinant Tet1 was 97.04%. For M.SssI assisted bisulfite sequencing (MAB-RRBS, to assay 5fC + 5caC levels), genomic DNA was methylated by M.SssI according to a previously published protocol , before bisulfite conversion. We assessed methylation efficiency by deep-sequencing MAB-treated genomic DNA containing unmethylated lambda phage DNA (1:1000 ratio) and by deep-sequencing MAB-treated lambda phage DNA alone. The estimated methylation efficiency is 96.20%; therefore, the expected error rate is 3.8%. We applied a binomial test to each called CpG site, accounting for the M.SssI methylation failure rate and sequencing coverage of that base. A P value for each CpG is calculated as 1-binom (number of times sequenced as T, total reads, fail rate). Because of the relatively low levels of 5fC and 5caC in the genome, we considered 5fC and 5caC only at sites with at least 20X coverage with a P value < 0.1. The coverage report can be found in Additional file 1: Table S20.
DNA epigenetic modification analysis at CpG and non-CpG sites
The phyloepigenetic evolutionary tree was generated using the top 150,000 most variably methylated CpGs by incorporating 154 TCGA human Infinium 450 K datasets from GBM patients with our sample set. We applied a cutoff of 0.15/0.02/0.02 for average 5mC/5hmC/5fC + 5caC changes at promoter regions (TSS1500) or a particular feature (enhancers) to assess patterns of differential methylation. A cutoff of 0.2/0.05/0.05 was used for individual site comparisons. CpGs were annotated with Bedtools to genomic features. DNA modification of a promoter CpG is associated with TET expression when the Pearson correlation P value is < 0.01 and the modification exhibits a difference of at least 0.1 for 5mC, 5hmC, and 5fC between the maximum and minimum values. To characterize the genome-wide non-CpG methylation pattern in NSC and GSC, and the effect of differentiation, we used non-CpG sites captured by RRBS with at least 10× coverage and analyzed them in a strand-specific manner .
Validation of DNA and histone epigenetic marks by chromatin and DNA immunoprecipitation coupled with real-time PCR (ChIP-qPCR and DIP-qPCR)
For ChIP-qPCR, samples were prepared as described above for ChIP-seq. For DIP-qPCR, DNA was extracted and incubated with 2 μL RNAse (GenDNA RNAse) at 37 °C for 10 min, followed by phenol-chloroform extraction. DNA samples were subject to shearing with a Covaris S220 to achieve 150-bp fragment size. Fragmented DNA was incubated with 2.5 μg 5mC (Epigentek), or 5 μg 5hmC (in-house) , 5fC (in-house), or 5caC (Cell Signaling) antibody, respectively, at room temperature for 3 h, followed by addition of 3 μL bridge antibody (Active Motif) and pre-washed magnetic protein beads at 4 °C overnight. The complex was washed with MeDIP (10 mM sodium phosphate, pH 7.0, 140 mM NaCl, 0.05% Triton X-100) and TE buffer, and then eluted at 65 °C for 15 min. Eluate was digested with proteinase K and purified using a Qiagen purification column. Eluted DNA was measured for enrichment at selected loci using primers in Additional file 1: Table S17. The specificity of each antibody was tested using spike-in controls; validation results are shown in Additional file 1: Table S21.
Motif search and gene ontology analysis
Motifs were identified using MEME-suite/CentriMO/HOMER software at the default settings. Matching frequencies to the input sequences (± 300 bp of each highly modified 5hmC or 5fC/5caC site) and significance assessment (E-value) were considered for selecting the best hits. Ontology analysis was performed using the DAVID bioinformatics database with Benjamini correction for multiple testing and Ingenuity Pathway Analysis (IPA, Qiagen) with standard program parameters. GREAT (version 3.0.0) was used to predict genes near enhancer histone marks within 50-kb distance and to generate functional annotation and basic ontology analysis for these genes (e.g. biological process, gene cluster).
Results of real-time PCR are expressed as the mean ± SEM. Significant differences were tested using one-way ANOVA on ranks. Post-hoc analysis was performed with pair-wise comparison using Dunn’s method. Differences in overall modification level among groups was performed with one-way ANOVA on ranks, with post-hoc analysis between each cell line performed with the Mann–Whitney U-test using Sigma Plot software 12.0 (Systat Software, Inc). Significance level was set at P < 0.05 level for all statistical tests.
We kindly thank Yue Wang for developing phython programs for data analysis, Ann Tuma and Alissa Caron for GBM tissue samples, and Leonard Collins for conducting LC/MS-MS.
This work was supported by a Pilot Project Grant from the Mayo Clinic SPORE in Brain Cancer (P50CA108961, KDR), the National Institutes of Health (R01CA114229, KDR), the Mayo Clinic Center for Individualized Medicine Epigenomics Program, and the Mayo Clinic Cancer Center.
Availability of data and materials
All authors had access to the study data and reviewed and approved the final manuscript. All data generated in this manuscript are included in NCBI GEO under accession number GSE92469 Superseries . Additional datasets used for analysis are previously published and deposited in GEO under accession numbers GSE17312  and GSE16368 . Level 3 DNA methylation Infinium 450 K array data, mRNA expression (RNA-seq data V2 RSEM), and clinical information for glioblastoma primary tumor are acquired from the Genomic Data Commons (https://portal.gdc.cancer.gov/) via the GDC client tool (https://gdc.cancer.gov/access-data/gdc-data-transfer-tool).
DZ designed and conducted laboratory experiments, performed data and statistical analyses, interpreted the results, wrote and revised the manuscript. BMA, MAS, SL, and JNS performed laboratory experiments and contributed to writing. JJT and RAH assisted with data analysis and revised the manuscript. JQ provided tissues. JHL conducted library preparation for ChIP-seq. KDR conceived and designed the study, interpreted the results, and wrote the manuscript. All authors have read and approved the final manuscript.
Ethics approval and consent to participate
All experimental animal procedures described in this work were reviewed and approved by the Mayo Clinic Institutional Animal Care and Use Committee (#A56814–14, “Creation, maintenance and characterization of brain tumor xenograft panel”). All research involving human tumor material was performed in accordance with the Declaration of Helsinki. Patient brain tumor samples for mouse xenograft generation were obtained with informed consent under a Mayo Clinic Institutional Review Board approved protocol (#07–007623, “Primary xenograft model for studying brain tumor biology and therapy”).
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Ostrom QT, Bauchet L, Davis FG, Deltour I, Fisher JL, Langer CE, et al. The epidemiology of glioma in adults: a “state of the science” review. Neuro-Oncology. 2014;16:896–913.View ArticlePubMedPubMed CentralGoogle Scholar
- Verhaak RG, Hoadley KA, Purdom E, Wang V, Qi Y, Wilkerson MD, et al. Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell. 2010;17:98–110.View ArticlePubMedPubMed CentralGoogle Scholar
- Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. 2014;344:1396–401.View ArticlePubMedPubMed CentralGoogle Scholar
- Sottoriva A, Spiteri I, Piccirillo SG, Touloumis A, Collins VP, Marioni JC, et al. Intratumor heterogeneity in human glioblastoma reflects cancer evolutionary dynamics. Proc Natl Acad Sci U S A. 2013;110:4009–14.View ArticlePubMedPubMed CentralGoogle Scholar
- Azari H, Millette S, Ansari S, Rahman M, Deleyrolle LP, Reynolds BA. Isolation and expansion of human glioblastoma multiforme tumor cells using the neurosphere assay. J Vis Exp. 2011;56:e3633.Google Scholar
- Moore N, Lyle S. Quiescent, slow-cycling stem cell populations in cancer: a review of the evidence and discussion of significance. J Oncol. 2011;2011:396076.View ArticlePubMedGoogle Scholar
- Campos B, Gal Z, Baader A, Schneider T, Sliwinski C, Gassel K, et al. Aberrant self-renewal and quiescence contribute to the aggressiveness of glioblastoma. J Pathol. 2014;234:23–33.View ArticlePubMedGoogle Scholar
- Wang Y, Yang J, Zheng H, Tomasek GJ, Zhang P, McKeever PE, et al. Expression of mutant p53 proteins implicates a lineage relationship between neural stem cells and malignant astrocytic glioma in a murine model. Cancer Cell. 2009;15:514–26.View ArticlePubMedPubMed CentralGoogle Scholar
- Zhu Y, Guignard F, Zhao D, Liu L, Burns DK, Mason RP, et al. Early inactivation of p53 tumor suppressor gene cooperating with NF1 loss induces malignant astrocytoma. Cancer Cell. 2005;8:119–30.View ArticlePubMedPubMed CentralGoogle Scholar
- Shukla S, Meeran SM. Epigenetics of cancer stem cells: Pathways and therapeutics. Biochim Biophys Acta. 2014;1840:3494–502.View ArticlePubMedGoogle Scholar
- Kraus TF, Globisch D, Wagner M, Eigenbrod S, Widmann D, Munzel M, et al. Low values of 5-hydroxymethylcytosine (5hmC), the “sixth base,” are associated with anaplasia in human brain tumors. Int J Cancer. 2012;131:1577–90.View ArticlePubMedGoogle Scholar
- Ficz G, Gribben JG. Loss of 5-hydroxymethylcytosine in cancer: cause or consequence? Genomics. 2014;104:352–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Noushmehr H, Weisenberger DJ, Diefes K, Phillips HS, Pujara K, Berman BP, et al. Identification of a CpG island methylator phenotype that defines a distinct subgroup of glioma. Cancer Cell. 2010;17:510–22.View ArticlePubMedPubMed CentralGoogle Scholar
- Johnson KC, Houseman EA, King JE, von Herrmann KM, Fadul CE, Christensen BC. 5-Hydroxymethylcytosine localizes to enhancer elements and is associated with survival in glioblastoma patients. Nat Commun. 2016;7:13177.View ArticlePubMedPubMed CentralGoogle Scholar
- Carlson BL, Pokorny JL, Schroeder MA, Sarkaria JN. Establishment, maintenance and in vitro and in vivo applications of primary human glioblastoma multiforme (GBM) xenograft models for translational biology studies and drug discovery. Curr Protoc Pharmacol. 2011;Chapter 14:Unit 14–6.Google Scholar
- Ito S, Shen L, Dai Q, Wu SC, Collins LB, Swenberg JA, et al. Tet proteins can convert 5-methylcytosine to 5-formylcytosine and 5-carboxylcytosine. Science. 2011;333:1300–3.View ArticlePubMedPubMed CentralGoogle Scholar
- Lee EJ, Rath P, Liu J, Ryu D, Pei L, Noonepalle SK, et al. Identification of global DNA methylation signatures in glioblastoma-derived cancer stem cells. J Genet Genomics. 2015;42:355–71.View ArticlePubMedPubMed CentralGoogle Scholar
- Rajbhandari R, McFarland BC, Patel A, Gerigk M, Gray GK, Fehling SC, et al. Loss of tumor suppressive microRNA-31 enhances TRADD/NF-kappaB signaling in glioblastoma. Oncotarget. 2015;6:17805–16.View ArticlePubMedPubMed CentralGoogle Scholar
- Shang C, Hong Y, Guo Y, Liu YH, Xue YX. MiR-210 up-regulation inhibits proliferation and induces apoptosis in glioma cells by targeting SIN3A. Med Sci Monit. 2014;20:2571–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Kaur H, Ali SZ, Huey L, Hutt-Cabezas M, Taylor I, Mao XG, et al. The transcriptional modulator HMGA2 promotes stemness and tumorigenicity in glioblastoma. Cancer Lett. 2016;377:55–64.View ArticlePubMedPubMed CentralGoogle Scholar
- Edimo WE, Ramos AR, Ghosh S, Vanderwinden JM, Erneux C. The SHIP2 interactor Myo1c is required for cell migration in 1321 N1 glioblastoma cells. Biochem Biophys Res Commun. 2016;476:508–14.View ArticlePubMedGoogle Scholar
- Weiler M, Pfenning PN, Thiepold AL, Blaes J, Jestaedt L, Gronych J, et al. Suppression of proinvasive RGS4 by mTOR inhibition optimizes glioma treatment. Oncogene. 2013;32:1099–109.View ArticlePubMedGoogle Scholar
- Creyghton MP, Cheng AW, Welstead GG, Kooistra T, Carey BW, Steine EJ, et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc Natl Acad Sci U S A. 2010;107:21931–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Jilaveanu LB, Parisi F, Barr ML, Zito CR, Cruz-Munoz W, Kerbel RS, et al. PLEKHA5 as a biomarker and potential mediator of melanoma brain metastasis. Clin Cancer Res. 2015;21:2138–47.View ArticlePubMedGoogle Scholar
- Eisele SC, Gill CM, Shankar GM, Brastianos PK. PLEKHA5: A key to unlock the blood-brain barrier? Clin Cancer Res. 2015;21:1978–80.View ArticlePubMedGoogle Scholar
- Park HJ, Chae SK, Kim JW, Yang SG, Jung JM, Kim MJ, et al. Ganglioside GM3 induces cumulus cell apoptosis through inhibition of epidermal growth factor receptor-mediated PI3K/AKT signaling pathways during in vitro maturation of pig oocytes. Mol Reprod Dev. 2017;84:702–11.View ArticlePubMedGoogle Scholar
- Kawashima N, Nishimiya Y, Takahata S, Nakayama KI. Induction of glycosphingolipid GM3 expression by valproic acid suppresses cancer cell growth. J Biol Chem. 2016;291:21424–33.View ArticlePubMedPubMed CentralGoogle Scholar
- Lee J, Kotliarova S, Kotliarov Y, Li A, Su Q, Donin NM, et al. Tumor stem cells derived from glioblastomas cultured in bFGF and EGF more closely mirror the phenotype and genotype of primary tumors than do serum-cultured cell lines. Cancer Cell. 2006;9:391–403.View ArticlePubMedGoogle Scholar
- Hahn MA, Li AX, Wu X, Pfeifer GP. Single base resolution analysis of 5-methylcytosine and 5-hydroxymethylcytosine by RRBS and TAB-RRBS. Methods Mol Biol. 2015;1238:273–87.View ArticlePubMedPubMed CentralGoogle Scholar
- Neri F, Incarnato D, Krepelova A, Parlato C, Oliviero S. Methylation-assisted bisulfite sequencing to simultaneously map 5fC and 5caC on a genome-wide scale for DNA demethylation analysis. Nat Protoc. 2016;11:1191–205.View ArticlePubMedGoogle Scholar
- Chang F, Xing P, Song F, Du X, Wang G, Chen K, et al. The role of T-box genes in the tumorigenesis and progression of cancer. Oncol Lett. 2016;12:4305–11.View ArticlePubMedPubMed CentralGoogle Scholar
- Thomas J, Morle L, Soulavie F, Laurencon A, Sagnol S, Durand B. Transcriptional control of genes involved in ciliogenesis: a first step in making cilia. Biol Cell. 2010;102:499–513.View ArticlePubMedGoogle Scholar
- Wu H, D’Alessio AC, Ito S, Wang Z, Cui K, Zhao K, et al. Genome-wide analysis of 5-hydroxymethylcytosine distribution reveals its dual function in transcriptional regulation in mouse embryonic stem cells. Genes Dev. 2011;25:679–84.View ArticlePubMedPubMed CentralGoogle Scholar
- Neri F, Incarnato D, Krepelova A, Rapelli S, Anselmi F, Parlato C, et al. Single-base resolution analysis of 5-formyl and 5-carboxyl cytosine reveals promoter DNA methylation dynamics. Cell Rep. 2015; https://doi.org/10.1016/j.celrep.2015.01.008.
- Stroud H, Feng S, Morey Kinney S, Pradhan S, Jacobsen SE. 5-Hydroxymethylcytosine is associated with enhancers and gene bodies in human embryonic stem cells. Genome Biol. 2011;12:R54.View ArticlePubMedPubMed CentralGoogle Scholar
- Sangel P, Oka M, Yoneda Y. The role of Importin-betas in the maintenance and lineage commitment of mouse embryonic stem cells. FEBS Open Bio. 2014;4:112–20.View ArticlePubMedPubMed CentralGoogle Scholar
- Omura N, Mizuma M, MacGregor A, Hong SM, Ayars M, Almario JA, et al. Overexpression of ankyrin1 promotes pancreatic cancer cell growth. Oncotarget. 2016;7:34977–87.View ArticlePubMedPubMed CentralGoogle Scholar
- Caren H, Stricker SH, Bulstrode H, Gagrica S, Johnstone E, Bartlett TE, et al. Glioblastoma stem cells respond to differentiation cues but fail to undergo commitment and terminal cell-cycle arrest. Stem Cell Reports. 2015;5:829–42.View ArticlePubMedPubMed CentralGoogle Scholar
- Labeed FH, Lu J, Mulhall HJ, Marchenko SA, Hoettges KF, Estrada LC, et al. Biophysical characteristics reveal neural stem cell differentiation potential. PLoS One. 2011;6:e25458.View ArticlePubMedPubMed CentralGoogle Scholar
- Suva ML, Rheinbay E, Gillespie SM, Patel AP, Wakimoto H, Rabkin SD, et al. Reconstructing and reprogramming the tumor-propagating potential of glioblastoma stem-like cells. Cell. 2014;157:580–94.View ArticlePubMedPubMed CentralGoogle Scholar
- Wu Y, Liu Y, Levine EM, Rao MS. Hes1 but not Hes5 regulates an astrocyte versus oligodendrocyte fate choice in glial restricted precursors. Dev Dyn. 2003;226:675–89.View ArticlePubMedGoogle Scholar
- Poloni A, Maurizi G, Foia F, Mondini E, Mattiucci D, Ambrogini P, et al. Glial-like differentiation potential of human mature adipocytes. J Mol Neurosci. 2015;55:91–8.View ArticlePubMedGoogle Scholar
- Yoshikawa H. Myelin-associated oligodendrocytic basic protein modulates the arrangement of radial growth of the axon and the radial component of myelin. Med Electron Microsc. 2001;34:160–4.View ArticlePubMedGoogle Scholar
- Marneros AG, Keene DR, Hansen U, Fukai N, Moulton K, Goletz PL, et al. Collagen XVIII/endostatin is essential for vision and retinal pigment epithelial function. EMBO J. 2004;23:89–99.View ArticlePubMedGoogle Scholar
- Cermenati G, Audano M, Giatti S, Carozzi V, Porretta-Serapiglia C, Pettinato E, et al. Lack of sterol regulatory element binding factor-1c imposes glial Fatty Acid utilization leading to peripheral neuropathy. Cell Metab. 2015;21:571–83.View ArticlePubMedGoogle Scholar
- Fathi E, Farahzadi R, Charoudeh HN. L-carnitine contributes to enhancement of neurogenesis from mesenchymal stem cells through Wnt/beta-catenin and PKA pathway. Exp Biol Med (Maywood). 2017;242:482–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Kafer GR, Li X, Horii T, Suetake I, Tajima S, Hatada I, et al. 5-Hydroxymethylcytosine marks sites of DNA damage and promotes genome stability. Cell Rep. 2016;14:1283–92.View ArticlePubMedGoogle Scholar
- Hecht SM. Bleomycin: new perspectives on the mechanism of action. J Nat Prod. 2000;63:158–68.View ArticlePubMedGoogle Scholar
- Jang HS, Shin WJ, Lee JE, Do JT. CpG and non-CpG methylation in epigenetic gene regulation and brain function. Genes (Basel). 2017;8:E148.View ArticleGoogle Scholar
- Ghanem N, Andrusiak MG, Svoboda D, Al Lafi SM, Julian LM, McClellan KA, et al. The Rb/E2F pathway modulates neurogenesis through direct regulation of the Dlx1/Dlx2 bigene cluster. J Neurosci. 2012;32:8219–30.View ArticlePubMedGoogle Scholar
- Boonen M, Staudt C, Gilis F, Oorschot V, Klumperman J, Jadot M. Cathepsin D and its newly identified transport receptor SEZ6L2 can modulate neurite outgrowth. J Cell Sci. 2016;129:557–68.View ArticlePubMedGoogle Scholar
- Putiri EL, Tiedemann RL, Thompson JJ, Liu C, Ho T, Choi JH, et al. Distinct and overlapping control of 5-methylcytosine and 5-hydroxymethylcytosine by the TET proteins in human cancer cells. Genome Biol. 2014;15:R81.View ArticlePubMedPubMed CentralGoogle Scholar
- Bao S, Wu Q, McLendon RE, Hao Y, Shi Q, Hjelmeland AB, et al. Glioma stem cells promote radioresistance by preferential activation of the DNA damage response. Nature. 2006;444:756–60.View ArticlePubMedGoogle Scholar
- Alcantara Llaguno SR, Wang Z, Sun D, Chen J, Xu J, Kim E, et al. Adult lineage-restricted CNS progenitors specify distinct glioblastoma subtypes. Cancer Cell. 2015;28:429–40.View ArticlePubMedPubMed CentralGoogle Scholar
- Stricker SH, Feber A, Engstrom PG, Caren H, Kurian KM, Takashima Y, et al. Widespread resetting of DNA methylation in glioblastoma-initiating cells suppresses malignant cellular behavior in a lineage-dependent manner. Genes Dev. 2013;27:654–69.View ArticlePubMedPubMed CentralGoogle Scholar
- Mazor T, Pankov A, Johnson BE, Hong C, Hamilton EG, Bell RJ, et al. DNA methylation and somatic mutations converge on the cell cycle and define similar evolutionary histories in brain tumors. Cancer Cell. 2015;28:307–17.View ArticlePubMedPubMed CentralGoogle Scholar
- Jin X, Kim LJY, Wu Q, Wallace LC, Prager BC, Sanvoranart T, et al. Targeting glioma stem cells through combined BMI1 and EZH2 inhibition. Nat Med. 2017;23:1352–61.PubMedGoogle Scholar
- Pollard SM, Yoshikawa K, Clarke ID, Danovi D, Stricker S, Russell R, et al. Glioma stem cell lines expanded in adherent culture have tumor-specific phenotypes and are suitable for chemical and genetic screens. Cell Stem Cell. 2009;4:568–80.View ArticlePubMedGoogle Scholar
- Woolard K, Fine HA. Glioma stem cells: better flat than round. Cell Stem Cell. 2009;4:466–7.View ArticlePubMedGoogle Scholar
- Schwartz PH, Bryant PJ, Fuja TJ, Su H, O’Dowd DK, Klassen H. Isolation and characterization of neural progenitor cells from post-mortem human cortex. J Neurosci Res. 2003;74:838–51.View ArticlePubMedGoogle Scholar
- Zhou D, Hlady RA, Schafer MJ, White TA, Liu C, Choi JH, et al. High fat diet and exercise lead to a disrupted and pathogenic DNA methylome in mouse liver. Epigenetics. 2017;12:55–69.View ArticlePubMedGoogle Scholar
- Sun Z, Baheti S, Middha S, Kanwar R, Zhang Y, Li X, et al. SAAP-RRBS: streamlined analysis and annotation pipeline for reduced representation bisulfite sequencing. Bioinformatics. 2012;28:2180–1.View ArticlePubMedPubMed CentralGoogle Scholar
- Guo JU, Su Y, Shin JH, Shin J, Li H, Xie B, et al. Distribution, recognition and regulation of non-CpG methylation in the adult mammalian brain. Nat Neurosci. 2014;17:215–22.View ArticlePubMedGoogle Scholar
- Tiedemann RL, Putiri EL, Lee JH, Hlady RA, Kashiwagi K, Ordog T, et al. Acute depletion redefines the division of labor among DNA methyltransferases in methylating the human genome. Cell Rep. 2014;9:1554–66.View ArticlePubMedPubMed CentralGoogle Scholar
- Zhou D, Robertson KD. GEO data series, accession GSE92469: Characterizing the epigenome of glioma stem cells. 2018. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE92469.Google Scholar
- Bernstein BE, Meissner A. GEO data series, accession GSE17312: BI Human Reference Epigenome Mapping Project. 2009. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE17312.Google Scholar
- BCCA Genome Sciences Centre. GEO data series, accession GSE16368: UCSF-UBC Human Reference Epigenome Mapping Project. 2009. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE16368.Google Scholar