Analysis of SOX2-Regulated Transcriptome in Glioma Stem Cells.

INTRODUCTION: Glioblastoma is the most malignant brain tumor in adults and is associated with poor survival despite multimodal treatments. Glioma stem-like cells (GSCs) are cells functionally defined by their self-renewal potential and the ability to reconstitute the original tumor upon orthotopic implantation. They have been postulated to be the culprit of glioma chemo- and radio-resistance ultimately leading to relapse. Understanding the molecular circuits governing the GSC compartment is essential. SOX2, a critical transcription regulator of embryonic and neural stem cell function, is deregulated in GSCs however; the precise molecular pathways regulated by this gene in GSCs remain poorly understood.
RESULTS: We performed a genome-wide analysis of SOX2-regulated transcripts in GSCs, using a microarray. We identified a total of 2048 differentially expressed coding transcripts and 261 non-coding transcripts. Cell adhesion and cell-cell signaling are among the most enriched terms using Gene Ontology (GO) classification. The pathways altered after SOX2 down-modulation includes multiple cellular processes such as amino-acid metabolism and intercellular signaling cascades. We also defined and classified the set of non-coding transcripts differentially expressed regulated by SOX2 in GSCs, and validated two of them.
CONCLUSIONS: We present a comprehensive analysis of the transcriptome controlled by SOX2 in GSCs, gaining insights in the understanding of the potential roles of SOX2 in glioblastoma.

Introduction

Glioblastoma is the most common and deadly primary brain tumors in adults and despite multiple treatments, the survival of glioma patients remains poor, with a median time between 12–15 months [1-3]

Glioblastoma contains a subpopulation of tumor propagating stem-like cells, known as glioma stem-like cells (GSCs) [4], which display the ability to self-renew, to differentiate into distinct lineages and to efficiently initiate and propagate tumors in xenografts models that recapitulate the phenotypic characteristics of the initial tumor from which they were derived [5,6]. Moreover, GSCs have been shown to increase resistance to radio-and chemotherapy [7,8], explaining in part the poor overall survival despite multiple treatments.

SOX2, a member of the SRY gene family, is a key transcription factor in the regulation of stemness properties, and it is essential in early embryonic development [9]. SOX2 has been reported to be deregulated in several human cancers [10-12] including glioblastoma where is over expressed due to several mechanisms such as amplification and promoter hypomethylation [13]. SOX2 is enriched in human-derived GSCs where it sustains stemness, migration, invasion and maintenance of tumorigenicity [13,14]. Although SOX2 response program in a glioblastoma cell line has been analyzed [15], to the best of our knowledge, an exhaustive analysis of SOX2-regulated molecular circuitries in GSCs has not been performed. Deciphering the molecular circuitries controlled by SOX2 in GSCs will provide insights about glioma development, biology and possible novel molecular therapies.

Besides coding genes, long non-coding RNAs (lncRNAs) are an emerging class of RNAs with no functional protein-coding ability that consists of more than 200 nucleotides [16]. Recent discoveries have proven that they play important roles regulating gene expression and function. These non-coding RNAs actively participate in many pathological processes in human malignancies [17-19] including cancer where a number of lncRNAs have been shown to act as oncogenes or tumor suppressors [20]. Recently, different groups published a signature of lncRNAs with aberrant expression in glioblastoma [21] and a set of prognostic lncRNAs that could have clinical implications in the sub-classification of this disease [22]; though the functional effect of lncRNAs in glioblastoma is not well understood.

Given that SOX2 is predominantly expressed in the GSCs compartment, which plays prominent roles in driving the growth, treatment resistance and recurrence of glioblastoma, the elucidation of the transcriptome and the molecular pathways involved in the generation and maintenance of this recalcitrant population is critical to understand the molecular underpinnings of glioblastoma malignancy. The aim of this work was to characterize the transcriptome regulated by SOX2 in GSCs. We set out to describe not only the coding genes but also the lncRNAs, which have been shown to play predominant roles in cancer. In this study we present a comprehensive analysis of the transcriptome controlled by SOX2 in GSCs, gaining insights in the understanding of the potential roles of SOX2 in glioblastoma.

Materials and Methods

Cell Lines, culture and transfection

Neurosphere cultures (GSC11 and GSC23), a kind gift of Dr. Lang at UT MD Anderson Cancer Center, were established from acute cell dissociation of human glioblastoma surgical specimens and maintained in Dulbecco's modified Eagle's medium/nutrient mixture F12 supplemented with B27 (Invitrogen, Carlsbad, CA), epidermal growth factor, and basic fibroblast growth factor (20 ng/mL each; Sigma-Aldrich, St Louis, MO).

To inhibit SOX2 expression, transient transfection assays were performed using two commercially available, specific siRNA against human SOX2 (si-SOX2, s13295 and s13294, Ambion) and a non-targeting control siRNA (si-scramble) (Ambion) in four independent experiments. The siRNA transfections were performed according to the manufacturer's instructions using Lipofectamine 2000 (Invitrogen). The cells were then cultured for 72 h after transfections and subjected to different analysis.

RNA extraction and Real Time PCR analysis

Total cellular RNA was isolated from the cultured cells using a Trizol reagent (Ambion) according to the manufacturers' protocols. For lncRNAs analysis, total RNA was subjected to DNase I treatment to digest the DNA. RNA quantity and quality were measured by NanoDrop ND-1000. The RNA samples were then reversely transcribed into cDNA using the Taqman miRNA Reverse Transcription Kit (Applied Biosystems) according to the manufacturer's instructions. Real Time PCR (RT-PCR) was performed using the Sybr Green Fast Master Mix (Applied Biosystems) in the ABI 7700 sequence detection system (Applied Biosystems, Foster City, CA). The quality of the products was controlled by the melting curve. Transcript levels were normalized against human GAPDH. Transcripts expression levels relative to GAPDH were calculated using the ddCt method. Primers for lncRNA detection and quantification were designed at Universal ProbeLibrary Assay Design Center (http://www.roche-applied-science.com/). All primer sequences are listed below (Tables 1 and 2):

Table 1

Primers of lncRNAs for qRT-PCR.

lncRNA position	TCONs	Forward Primer	Reverse Primer
chr19:28281401–28284848	TCONS_00027256	GCCCAAAGTTTGATTTCTCG	CGAGGTCTAACCCAGGTGTG
chr11:121899032–121899389	TCONS_00020142	GCTGAGCCTTCCATGAAAAT	GTGCAAATCACTCCAGTCACA

Table 2

Primers of genes for qRT-PCR.

Gene	Forward Primer	Reverse Primer
GAPDH	AGCCACATCGCTCAGACAC	GCCCAATACGACCAAATCC
SOX2	AGCTCGCAGACCTACATGAA	CCGGGGAGATACATGCTGAT
PLP1	ACCTATGCCCTGACCGTTG	TGCTGGGGAAGGCAATAGACT
COL2A1	TGGACGCCATGAAGGTTTTCT	TGGGAGCCAGATTGTCATCTC
ATP8B1	ACGACATTTGACGAGGATTCTC	GGTTTTGTTCTGGTTCAACAGC
PPP1R1B	CAAGTCGAAGAGACCCAACCC	GCCTGGTTCTCATTCAAATTGCT
CMTM5	GGAGGACCACATCCGCTAGAT	CCAGGGAGTGGAAGCAGAT
GALNT14	CACTGCTGGTGTATTGCACG	CGGATCAGATGCGTAGGGG
F11R	GTGCCTACTCGGGCTTTTCTT	GTCACCCGGTCCTCATAGGAA
SYT4	ATGGGATACCCTACACCCAAAT	TCCCGAGAGAGGAATTAGAACTT
SLC18A1	GTGGTGGTATTCGTCGCTTTG	CCGAGGTGCAGAGAAGAGT
ITLN2	GCAGGGCAACAAAGCAGACTA	CAGGGCGCTGTTTCTCCAA

Immunoblotting Assay

For the western blot assay, cells were lysed in RIPA buffer (Triton and PBS) for 30 min on ice. Samples containing identical amounts of protein (30 μg) were resolved in a 12% polyacrylamide gel, transferred to polyvinylidene membranes, and blocked in 5% nonfat milk in phosphate-buffered saline/Tween-20. Membranes were incubated with the following antibodies: SOX2 (Cell Signaling, Danvers, MA) and α-Tubulin (Sigma-Aldrich) using 1:1000 dilution. The membranes were developed according to the protocol for enhanced chemiluminiscence from Perkin Elmer.

Microarray expression analysis

Total RNA was isolated from scrambled and SOX2-siRNA GSC11 cells using Trizol extraction. RNA was purified by the QIAGEN RNAeasy mini kit (QIAGEN) according to the manufacturer´s protocol. One-color Cy3 RNA labeling, array hybridization to Agilent SurePrint G3 8 × 60 K Human Gene Expression Arrays (Agilent Technologies), data collection, and analysis were performed at the Bioinformatics Unit (Fundación para la Investigación Médica Aplicada, CIMA, Pamplona, Spain). Normalization of microarray data was performed using quantile algorithm. After quality assessment a filtering process was carried out to eliminate low expression probe sets. Applying the criterion of an expression value greater than 64 in 2 samples of at least one of the experimental conditions, 40986 probe sets were selected for statistical analysis. LIMMA (Linear Models for Microarray Data) [23] was used to find out the probe sets that showed significant differential expression between experimental conditions. Genes were selected as significant using a B statistic cut off B>0. Data processing and statistical analyses were performed with R and Bioconductor [24]. The microarray data from this study are publicly available at Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo) under accession number GSE79302.

Functional group analysis

In this study we applied Gene Ontology (GO) analysis to find the primary function of the differential expression of mRNAs regulated by SOX2, using online software DAVID (Database for Annotation, Visualization and Integrated Discovery, http://david.abcc.ncifcrf.gov/). GO analysis can organized genes into hierarchical categories (Gene Ontology Consortium). To find out the significant pathway of the differential genes participating we performed gene regulatory network analysis using Ingenuity Pathway Analysis (IPA) software (http://www.ingenuity.com), which can integrate gene-expression data with other molecular databases to facilitate the development of new and more complete pathway maps. Fisher's exact test was used to select the significant GO categories. The threshold of significance was defined by P value with a cut-off set in 0.05.

Statistical analysis

Experimental data are represented as the mean ± SD of three biologic replicates and were compared using Student's t-test. Significant P-values are indicated with asterisks as follows: *P < 0.05, **P < 0.01.

Results

Transcripts regulated by SOX2

Since SOX2 is a key driver in the maintenance of the GSCs phenotype and therefore in the perpetuation of this devastating tumor we down-regulated the expression of this gene in the GSC11 cell line in four independent experiments, using a SOX2 specific siRNA (Fig 1A). The efficiency of SOX2 knockdown was assessed by real-time PCR and western blot (Fig 1B) and was also confirmed in our array results. Microarray data identified a total of 2048 differentially expressed coding transcripts and 261 non-coding transcripts (B value >0) (Fig 1A and S1 Table).

Fig 1

Transcripts regulated by SOX2.

(A) Schematic representation of the research design employed to uncover the SOX2 transcriptome in GSC11 cells. (B) qRT-PCR and western blot confirmation of SOX2 inhibition in GSC11 cells after 72h of si-SOX2 or si-Scramble (si-Scrbl) transfection. SOX2 relative mRNA levels are presented as 2-ΔΔCt standardized with their constitutive gene GAPDH. Each bar represents the mean ± SD. For western blot tubulin was used as housekeeping control and shown as a representative blot of four independent experiments.

SOX2 controls a wide spectrum of protein-coding genes and pathways in GSCs

To further narrow the coding transcripts data a cut-off 1 logarithmic fold difference between SOX2 knockdown and scrambled GSC11 cells was set, identifying 35 up-regulated and 100 down-regulated genes, which suggest that SOX2 act primarily as a transcriptional activator. In Table 3 we showed the top-10 up or down-regulated protein coding-genes, and select the top 5 candidates of each group for further validation by qRT-PCR using the GSC11 and GSC23 cell lines. We confirmed the observed microarray expression changes in 5 out of 5 down-regulated coding-genes (Fig 2B) and in one out of 5 up-regulated coding genes in GSC11 cells (Fig 2A). Regarding GSC-23 cells, we down-modulate SOX2 expression using two different siRNAs against human SOX2, and we validate the expression of 5 out of 5 up-regulated coding-genes (Fig 2A) and of 4 out of 5 down-regulated coding genes (Fig 2B), partially validating the microarray results.

Table 3

The top 10 up- and down-regulated genes in Sox2-downmodulated GSC-11 cells, organized by logFC.

GeneName	logFC	P.Value	B
PLP1	2,443	4,73E-05	2,488
COL2A1	2,122	2,18E-05	3,281
ATP8B1	1,954	1,50E-06	5,943
PPP1R1B	1,925	7,27E-06	4,393
CMTM5	1,702	1,12E-04	1,586
ELMO1	1,66	5,97E-05	2,246
ITIH5L	1,555	6,67E-05	2,133
IGFBP5	1,522	1,82E-06	5,759
SCARNA9	1,498	7,10E-05	2,067
SCARNA17	1,435	4,46E-06	4,879
GALNT14	-3,006	0,000299	0,569
F11R	-2,824	1,13E-07	8,308
SYT4	-2,664	2,61E-06	5,405
SLC18A1	-2,335	4,01E-06	4,984
ITLN2	-2,329	3,57E-07	7,288
RASEF	-2,294	9,38E-07	6,39
GADD45G	-2,231	3,55E-08	9,277
CYP26A1	-2,056	7,66E-10	11,938
KRTAP21-1	-2,028	0,00021	0,94
PNLIPRP2	-1,992	0,000495	0,042

Fig 2

Analysis by qRT-PCR of the top 5 (A) up- and (B) down-regulated coding transcripts in SOX2 downmodulated GSC11 and GSC23 cells. Total RNA was extracted after 72h of si-Sc or si-SOX2 transfection in GSC11 and GSC23 cells. In GSC23 cells we used two different siRNAs against human SOX2, si-SOX2 (1) is referred to s13295 and si-SOX2 (2) is referred to s13294 from Ambion. Values are normalized to GAPDH and each bar represents the mean ± SD.

To understand the significance of differential gene expression, bioinformatics analysis related to Gene Ontology (GO) Classification and pathway analysis were performed. GO classifications using the DAVID web tool and pathway analysis using Ingenuity Pathway Analysis (IPA) was performed. For these analysis, gene lists were classified based upon decreased (logFC < -1) or increased (logFC > 1) expression relative to control and analyzed altogether as a single list.

Enrichment analysis of GO categories including biological process (BP), molecular function (MF), and cellular component (CC) were obtained using DAVID web tool (Fig 3). We observed the highest enrichment in the categories related to “cell adhesion”, “biological adhesion”, “cell-cell signaling”, “extracellular region” and “calcium ion binding”.

Fig 3

Top-10 GO Biological Processes analysis of protein-coding genes regulated by SOX2 in GSC11 cells.

Bar chart represents classification of GO Biological Processes, Cellular Component or Molecular Function as determined by DAVID web tool. Bars represent the number of genes in the specified category, organized by p-value.

We used IPA analysis to undercover the canonical pathways regulated by SOX2 in GSCs. Our results showed 13 pathways significantly altered (Table 4). Most of them related with amino-acid metabolism pathways, such as “histamine biosynthesis”, where histamine is an important regulator of numerous physiological processes including neurotransmission in the central nervous system (CNS) [25]; “L-cysteine degradation process” where cystathionine γ-lyase (CTH) activity has been related with glioblastoma treatment [26] and “serotonin receptor signaling pathway” being serotonin an important neurotransmitter in the CNS during neuronal development [27]. Other enriched pathways were “hematopoiesis from multipotent stem cells”, where KITLG has been reported to regulate neoplastic processes such as growth and invasion [28]; apoptosis [29] and cell adhesion [30]. Role of JAK2 in Hormone-like Cytokine Signaling stood out because GHR (growth hormone receptor) and IRS1 (insulin receptor substrate 1) has been linked with glioma progression [31,32]. A well-characterized pathway frequently altered in tumors is the NOTCH signaling cascade, which was also enriched in our analysis. The NOTCH pathway is a conserved intercellular signaling route that has been implicated in different developmental processes. Interestingly, NOTCH pathway is deregulated in human glioblastoma and plays a key role in maintaining the growth, the undifferentiated state of glioma cells and tumorigenesis [33-35]. The integrated analysis of SOX2 enriched canonical pathways revealed the link between this transcription factor and multiple cellular processes such as amino-acid metabolism and intercellular signaling cascades, like NOTCH pathway.

Table 4

List of top-13 canonical pathways identified by IPA software.

Pathway	-log(p-value)	Ratio	Molecules
Glycine Betaine Degradation	2,63E+00	2,50E-01	DMGDH,PIPOX
Hepatic Stellate Cell Activation	2,38E+00	4,42E-02	LY96,COL2A1,COL22A1,IGFBP5,COL28A1
Histamine Biosynthesis	2,03E+00	1,00E+00	HDC
L-cysteine Degradation II	2,03E+00	1,00E+00	CTH
Triacylglycerol Degradation	2,02E+00	1,25E-01	PNLIPRP2,CES1
Retinol Biosynthesis	2,02E+00	1,25E-01	PNLIPRP2,CES1
Hematopoiesis from Multipotent Stem Cells	1,73E+00	5,00E-01	KITLG
Cysteine Biosynthesis/Homocysteine Degradation	1,73E+00	5,00E-01	CTH
Role of JAK2 in Hormone-like Cytokine Signaling	1,61E+00	7,69E-02	GHR,IRS1
Serotonin Receptor Signaling	1,58E+00	7,41E-02	SLC18A1,HTR1D
Phenylethylamine Degradation I	1,56E+00	3,33E-01	AOC3
Notch Signaling	1,35E+00	5,56E-02	HES5,HEY1
Lysine Degradation V	1,34E+00	2,00E-01	PIPOX

This selection is organized by the negative logarithm of p-values (Fisher Test), calculated by IPA ([-Log (0.05) = 1.3]).

The IPA analysis also showed the most relevant biological functions and diseases in our data set. The most significant bio-functions altered following SOX2 down-modulation are showed in Table 5. The set of SOX2-associated genes were assigned mainly to the following networks: “cancer”, “organismal injury and abnormalities”, “cellular movement”, “tissue morphology”, “cellular development” and “hematopoiesis”. Interestingly, most of these networks involved very well-known functions of SOX2 such as morphology determination [36], development [37] and cellular proliferation and migration in glioma [13,13]. Fig 4 shows the most relevant selection of bio-function categories: disease and disorders, molecular and cellular functions and physiological system development and function, obtained by using IPA software organized by p-value.

Table 5

The top ten significant Bio-Functions altered following Sox2 down-modulation in the GSC11 cell line.

Category	p-value	Number of Targets
Dermatological Diseases and Conditions	1,63 x 10−08–9,35 x 10−03	74
Cancer	2,93 x 10−08–9,35 x 10−03	93
Organismal Injury and Abnormalities	2,93 x 10−08–9,35 x 10−03	100
Cellular Movement	3,98 x 10−06–9,35 x 10−03	41
Connective Tissue Development and Function	1,71 x 10−05–9,35 x 10−03	23
Tissue Morphology	1,71 x 10−05–9,35 x 10−03	39
Reproductive System Disease	2,52 x 10−05–9,35 x 10−03	16
Cellular Development	2,72 x 10−05–9,35 x 10−03	51
Hematological System Development and Function	2,72 x 10−05–9,35 x 10−03	31
Hematopoiesis	2,72 x 10−05–9,35 x 10−03	9

The p-value range indicates the p-values of the various pathways and processes belonging to that category. The number of targets indicates the total number of genes associated with the functional category.

Fig 4

10-Top Bio Functions categories altered following SOX2 inhibition.

The categories listed are Physiological System Development and Function, Molecular and cellular Functions and Disease and Disorders, identified using IPA software. Bars represent the number of genes in the specified category, organized by p-value.

These results established a signature of protein coding-genes regulated by SOX2 in GSCs with biological functions relevant to glioblastoma growth and maintenance of its malignant phenotype. The tight overlap between the existing literature and our enrichment analysis highlights the robustness of our results and predicts that this approach will be an excellent discovery platform to identify novel SOX2 targets.

SOX2–regulated non-coding RNAs in GSC

Reprogramming transcription factors, including SOX2, have been shown to regulate both coding and non-coding RNAs [38]. LncRNAs are emerging as key regulators of biological processes and disease [39] therefore, seems reasonable to hypothesize that SOX2 will regulate this class of genes as well. The strength of our data-sets allowed us to identify potential non-coding transcripts differentially expressed (B value > 0) regulated by SOX2 in GSCs. After biotype distribution analysis we identify protein coding RNAs (44% for up-regulated and 41% for down-regulated), while the rest were classified as different types of non-coding transcripts. Out of the total number of transcripts differentially expressed we identify 80 upregulated and 181 down-regulated and we classify them as intergenic RNAs, antisense, processed transcripts, transcripts derived from pseudogenes and unassigned transcripts (Fig 5). The transcripts classified as “others” correspond to transcripts derived from miRNAs, rRNAs, sense-overlaping and sense intronic transcripts. The lncRNA annotation was performed with the Bioconductor package ChIPpeakAnno [40] and using Gencode v19 as reference [41]. The gene type corresponding to the gene that overlaps with the lincRNA locus was assigned to each lincRNA. Table 6 shows the top 25 non-coding transcripts regulated by SOX2 in GSCs.

Fig 5

SOX2 regulated non-coding transcripts.

(A) A total of 261 transcripts were found differentially expressed (B > 0), which were distributed in 80 upregulated and 181 downregulated transcripts. (B) Biotype distribution of the differentially expressed transcripts following SOX2 down-modulation in GSC11 cells.

Table 6

List of the top 25 non-coding transcripts regulated by Sox2, organized by B value.

Probe	GeneName	Classification	logFC	B
A_19_P00320471	lincRNA:chr9:2535671–2536375_R	antisense	-1,28	11,75
A_19_P00315804	lincRNA:chr9:2530903–2539456_R	antisense	-1,15	10,77
A_19_P00320469	lincRNA:chr9:2535671–2536375_R	antisense	-1,27	9,70
A_19_P00811613	lincRNA:chr9:2452800–2552025_R	antisense	-1,17	7,94
A_33_P3397743	LOC100128088	pseudogene	-1,83	7,42
A_19_P00321203	lincRNA:chr6:72126155–72129954_R	lincRNA	-0,92	7,33
A_19_P00322118	lincRNA:chr2:39745746–39826668_F	antisense	-0,86	6,90
A_23_P3552	LOC730092	pseudogene	-0,65	6,50
A_19_P00322220	lincRNA:chr20:37055062–37063887_R	processed_transcript	-0,69	6,28
A_33_P3392460	LOC100128077	processed_transcript	-1,47	6,18
A_19_P00322149	lincRNA:chr6:72126142–72129923_R	lincRNA	-0,91	6,16
A_19_P00317793	lincRNA:chr20:37055062–37063916_R	processed_transcript	-0,68	6,10
A_19_P00808846	lincRNA:chr21:17992729–18010729_F	lincRNA	-0,64	6,09
A_19_P00318304	lincRNA:chr20:37050986–37063998_R	processed_transcript	-0,67	5,94
A_19_P00316341	lincRNA:chr7:130600800–130606702_F	lincRNA	-0,86	5,93
A_19_P00316985	lincRNA:chr6:72126162–72129969_R	lincRNA	-0,91	5,93
A_19_P00322967	lincRNA:chr20:37050934–37057222_R	processed_transcript	-0,69	5,54
A_19_P00802098	lincRNA:chr2:3579550–3585150_R	lincRNA	-0,58	5,15
A_24_P756289	SOX2OT	other	-0,86	5,08
A_33_P3613516	LOC254057	antisense	-1,10	4,94
A_19_P00318174	lincRNA:chr2:3579840–3584422_R	lincRNA	-0,73	4,52
A_33_P3287710	chr10:79,686,570–79,689,583	unassigned	-0,67	4,46
A_33_P3405043	LOC100133264	unassigned	-0,72	4,43
A_32_P88349	LOC730256	pseudogene	-0,48	4,36
A_33_P3705884	chr19:28,281,401–28,284,848	lincRNA	-0,89	4,27
A_19_P00321044	lincRNA:chr16:50682543–50683160_F	lincRNA	1,04	6,09
A_19_P00315647	lincRNA:chr11:121899032–121899389_R	other	0,65	5,63
A_32_P63013	LOC283174	unassigned	1,32	5,08
A_32_P47157	LOC92973	unassigned	0,74	4,95
A_19_P00317484	lincRNA:chr3:112315643–112316945_R	lincRNA	0,55	4,11
A_19_P00809440	lincRNA:chr11:133765815–133774297_R	other	1,23	4,06
A_33_P3789382	chr10:65,224,989–65,226,322	antisense	0,57	3,98
A_19_P00321420	lincRNA:chr11:133766329–133767054_R	unassigned	1,34	3,74
A_19_P00332120	lincRNA:chr3:156455706–156471081_R	lincRNA	0,59	3,57
A_19_P00320101	lincRNA:chr11:133767609–133771496_R	other	1,06	3,51
A_19_P00812924	lincRNA:chr11:121895965–121904065_R	other	0,53	3,46
A_19_P00326763	lincRNA:chr3:112308735–112318605_R	lincRNA	0,48	2,73
A_33_P3753757	LOC158402	other	0,53	2,46
A_33_P3393679	LOC645323	lincRNA	0,47	2,24
A_19_P00315649	lincRNA:chr11:121899032–121899389_R	other	0,55	2,23
A_19_P00809838	lincRNA:chrX:100247844–100257469_R	unassigned	0,79	2,20
A_19_P00331576	lincRNA:chr3:114043485–114052926_F	unassigned	0,36	2,00
A_33_P3259557	LOC440104	pseudogene	0,45	1,67
A_19_P00319347	lincRNA:chr2:168149680–168414843_F	lincRNA	0,77	1,67
A_19_P00320212	lincRNA:chr9:114795825–114797203_R	other	0,41	1,63
A_24_P93703	LOC440104	pseudogene	0,40	1,50
A_19_P00318878	lincRNA:chr1:247350513–247352101_R	lincRNA	0,38	1,44
A_19_P00316010	lincRNA:chr17:67547498–67549996_F	lincRNA	0,57	1,43
A_24_P349207	ENST00000380727	pseudogene	0,29	1,31
A_19_P00802064	lincRNA:chr8:2522118–2527693_R	lincRNA	0,38	1,04

The top four differentially expressed lncRNAs (Table 7) that presented chromatin marks and high abundance in brain were selected using GRCh37/hg19 assembly in UCSC Genome Browser, and their expression was validated using qRT-PCR in GSC11 cells, comparing SOX2-siRNA versus scrambled control. Our data indicated that the expression of chr19:28,281,401–28,284,848 (TCONS_00027256) was significantly down-regulated (p value = 0,018), while chr11:121899032–121899389 (TCONS_00020142) was significantly up-regulated (p value = 0.042) after SOX2 inhibition (Fig 6A). The results were consistent with the microarray data (Fig 6B).

Table 7

List of the top-four lncRNAs regulated by Sox2, organized by p-value.

lncRNA	logFC	P.Value	B
lincRNA:chr6:72126155–72129954	-0,9188744	3,40E-07	7,3327381
lincRNA:chr6:29701971–29740296	-0,78724372	3,93E-07	7,2007586
lincRNA: chr19:28,281,401–28,284,848	-0,89247792	8,26E-06	4,265306
lincRNA:chr11:121899032–121899389	0,64633267	2,07E-06	5,6334606

This selection was evaluated by the presence of histone modifications and high abundance in brain according to UCSF genome browser tool.

Fig 6

Validation of two lncRNAs regulated by SOX2 in GSC11 cells.

(A) The expression of the transcripts located in chr11:121899032–121899389 (TCONS_00020142) and chr19:28,281,401–28,284,848 (TCONS_00027256) were assessed. In both cases GSC11 cells were transfected with siRNA control or siRNA against SOX2 and three days later RNA was extracted and subject to RT-PCR. Values are normalized to GAPDH and are the mean ± SD of three replicates. (B) Comparison between microarray and qRT-PCR results. The height of each column in this graph represents the log-transformed mean fold changes in the expression of lncRNA between Scramble and siSOX2 transfected cell line.

Altogether, these results identified and confirmed the non-coding transcript profile controlled by SOX2 in GSCs. Characterizing the functional relevance of these lncRNAs will undoubtedly impact our understanding of glioblastoma biology.

Discussion

Our work provides a comprehensive view of the genome wide SOX2 regulated transcripts in GSCs, illustrating a complex scenario where SOX2 is the central player regulating different molecules and pathways in glioblastoma.

In this study, we used state-of-the art microarray technology to query the SOX2 coding and non-coding RNA transcriptome in GSCs. It is interesting to note that among the down-regulated genes following SOX2 knockdown, F11R has been shown to be overexpressed in glioblastoma cells [42,43]. F11R is necessary and sufficient for GSC maintenance and self-renewal and of clinical significance is associated with increased malignancy and poor patient prognosis [42,43]. On the other hand, we found several interesting over-expressed candidates controlled by SOX2; for example, PPP1R1B is a well-known striatal projection neuron signature marker [44]. The fact that its expression increases following SOX2 inhibition is in line with its role in neuronal differentiation.

In a previous work where the SOX2 response program in a glioblastoma cell line was analyzed [15], authors identified 489 genes whose expression were altered in response to SOX2 knockdown, using several genomic technologies. Interestingly several of these genes are also differentially expressed in our array data, such as NGFR, CEBPA, BCL2, BNIP3, EBF4, ALCAM, protocadherins and solute carrier family members. Overall these results highlight the strength of our array data and the link between SOX2 and GSCs biology. The work of Fang and colleagues exhibits some similarities with our study, such as the analysis of SOX2 regulated-coding genes. However, we focused in the molecular circuitries controlled by SOX2 in GSCs, meanwhile they performed their study in an established glioma cell line. Additionally, our study addressed the SOX2-regulated lncRNAs in GSCs. Altogether both works provide clues regarding SOX2 functions in glioblastoma.

Gene-set enrichment analysis shows SOX2 is involved in regulating “cell adhesion”, “biological adhesion”, “cell-cell signaling”, and “calcium ion binding” pathways, undercovering its key function as a driver of the glioma stem-like phenotype [45-48].

We also analyzed the canonical pathways regulated by SOX2 in GSCs. Pathways related with amino-acid metabolism were among the most deregulated, illustrating that SOX2 expression is critical for maintaining metabolic homeostasis in the GSC population, and plays important role in different tumor microenvironment conditions, such as hypoxic stress conditions [49]. Other enriched pathway altered in our analysis was the NOTCH pathway, where Hes5 and Hey1 had the most significantly down-modulated expression. Hes5 is a marker of neural multipotent progenitors with stem cell properties [50] where it sustains progenitors proliferative state inhibiting their differentiation into neurons [51]. On the other hand Hey1 has been related to a subset of molecules directly associated with hypoxia in glioblastoma tumors [52]; and might be used as a marker to distinguish glioblastoma patients with a relative good prognosis (negative Hey1 expression) [53]. Furthermore, Hey 1 is up-regulated in glioma samples correlating with tumor grade, and functionally its down-regulation results in a proliferation reduction [54], suggesting a role in the progression of glioblastoma. Taking all this into account, the canonical pathways more significantly altered after SOX2 inhibition are those related with intracellular signaling cascades and amino-acid metabolism pathways associated with tumor propagation.

Consistent with what is known about SOX2 biological function, our data-set is enriched with genes involved in morphology determination, development and cellular proliferation and migration. Interestingly, “proliferation of tumor cells” (S2 Table), is one of the most repeated subcategories for which IPA analysis assigned an activation Z score close to -2, predicting it´s inhibition, which is in line with a putative role of SOX2 in cell proliferation.

One of the most exciting aspects of our study involved expanding our knowledge of the SOX2 transcriptome into the realm of lncRNAs. In this study we showed and classified the lncRNA landscape regulated by SOX2 in GSCs. Our microarray results showed a strong correlation with published reports, demonstrating the strength of our approach and providing confidence that we can use this data-set for de novo discovery of novel SOX2 targets, including lncRNAs. One previous study determined the differentially expressed lncRNAs between glioblastoma and brain tissues, showing 654 lncRNAs upregulated and 654 down-regulated [55]. To our knowledge, this is the first study that evaluates the differentially expression of lncRNAs in GSCs controlled by SOX2.

Among the transcripts regulated by SOX2 we found that SOX2OT was down-regulated in our data set, even though we did not validate it. SOX2OT is a lncRNA which harbors SOX2 gene in its intronic region and is transcribed in the same orientation as SOX2 [56]. Several studies have demonstrated a role of SOX2OT in the regulation of SOX2 gene in human stem cells [57,58] although little is known about the exact role of this non-coding RNA. SOX2OT has been associated with carcinogenesis and, for example in breast cancer is involved in the induction and/or maintenance of SOX2 expression [59], in esophageal squamous cell carcinoma has been shown to play a role in tumor initiation and/or progression as well as in regulation of the pluripotent state of stem cells [58], and proliferation in lung cancer cells [60]. Askarian-Amiri et al demonstrated that SOX2OT has a positive effect on SOX2 expression [59]. Published data suggest the mediation of lncRNA SOX2OT in pluripotency and tumorigenesis events, probably through regulation of SOX2 expression. These data together with our own results suggest a possible role of SOX2OT in the malignant phenotype of glioblastoma, however further functional and mechanistic studies will be necessary to elucidate the precise role of SOX2OT and other lncRNA candidates in the tumorigenicity of glioblastoma.

Although advance in managing and treating glioblastoma have been made, tumor recurrence and treatment resistance remains the major cause of glioblastoma mortality. Understanding how factors such as SOX2 drive the glioblastoma tumor phenotype will aid in the development of new therapeutic approaches based on targeting GSCs. Our study integrates for the first time the coding and non-coding transcriptome controlled by SOX2 in GSCs, gaining new insights about the molecular circuitries governing glioblastoma biology.

Conclusion

In conclusion we have performed a comprehensive analysis of differential expression of coding and non-coding transcripts controlled by SOX2 in GSCs. We performed gene set enrichment analysis to find the most relevant pathways and biological functions altered in our data set. This integrated analysis allows for a better understanding of the SOX2 transcriptome in GSCs.

Differentially expressed coding transcripts from microarray data.

(XLS)

Click here for additional data file.

List of disease and biofunctions data from Ingenuity analysis.

(XLS)

Click here for additional data file.