MOLONC664_proof ■ 30 March 2015 ■ 1/18 M O L E C U L A R O N C O L O G Y X X X ( 2 0 1 5 ) 1 e1 8
available at www.sciencedirect.com 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
ScienceDirect www.elsevier.com/locate/molonc
Q5
Q4
Comprehensive analysis of microRNA expression profile in malignant glioma tissues Monika Piweckaa,1, Katarzyna Rollea,1, Agnieszka Beltera, c _ Anna Maria Barciszewskab, Marek Zywicki , Marcin Michalakd, Stanis1aw Nowakb, Miros1awa Z. Naskre˛t-Barciszewskaa, Jan Barciszewskia,* a
Institute of Bioorganic Chemistry, Polish Academy of Sciences, Poznan, Poland Department of Neurosurgery and Neurotraumatology, Poznan University of Medical Sciences, Poznan, Poland c Bioinformatics Laboratory, Institute of Molecular Biology and Biotechnology, Faculty of Biology, A. Mickiewicz University, Poznan, Poland d Institute of Dendrology, Polish Academy of Sciences, Kornik, Poland b
A R T I C L E
I N F O
A B S T R A C T
Article history:
Malignant gliomas represent the most devastating group of brain tumors in adults, among
Received 16 October 2014
which glioblastoma multiforme (GBM) exhibits the highest malignancy rate. Despite com-
Received in revised form
bined modality treatment, GBM recurs and is invariably fatal. A further insight into the mo-
26 February 2015
lecular background of gliomagenesis is required to improve patient outcomes. The primary
Accepted 16 March 2015
aim of this study was to gain broad information on the miRNA expression pattern in ma-
Available online -
lignant gliomas, mainly GBM. We investigated the global miRNA profile of malignant glioma tissues with miRNA microarrays, deep sequencing and meta-analysis. We selected
Keywords:
miRNAs that were most frequently deregulated in glioblastoma tissues, as well as in peri-
Brain tumors
tumoral areas, in comparison with normal human brain. We identified candidate miRNAs
Glioblastoma
associated with the progression from glioma grade III to glioma grade IV. The meta-
Gliomas
analysis of miRNA profiling studies in GBM tissues summarizes the past and recent ad-
Meta-Analysis
vances in the investigation of the miRNA signature in GBM versus noncancerous human
microRNA
brain and provides a comprehensive overview. We propose a list of 35 miRNAs whose
miRNA profiling
expression is most frequently deregulated in GBM patients and of 30 miRNA candidates recognized as novel GBM biomarkers. ª 2015 Published by Elsevier B.V. on behalf of Federation of European Biochemical Societies.
1.
Introduction
Malignant gliomas are the most common and aggressive primary brain tumors in adults. WHO grade IV glioblastoma multiforme (GBM) is the most frequent type and is characterized
by strong vascular proliferation, invasiveness, diminished apoptosis, radio- and chemoresistance. Despite treatment, the median survival of GBM patients is approximately 12 months after diagnosis (Cheng et al., 2010; Stupp et al., 2005).
* Corresponding author. Tel.: þ48 61 8528503x132; fax: þ48 61 8528919. E-mail address: [email protected] (J. Barciszewski). 1 These authors contributed equally to this work. http://dx.doi.org/10.1016/j.molonc.2015.03.007 1574-7891/ª 2015 Published by Elsevier B.V. on behalf of Federation of European Biochemical Societies.
Please cite this article in press as: Piwecka, M., et al., Comprehensive analysis of microRNA expression profile in malignant glioma tissues, Molecular Oncology (2015), http://dx.doi.org/10.1016/j.molonc.2015.03.007
Q1
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
MOLONC664_proof ■ 30 March 2015 ■ 2/18
2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
M O L E C U L A R O N C O L O G Y X X X ( 2 0 1 5 ) 1 e1 8
Treatment of GBM remains one of the most challenging tasks in clinical oncology. It is the consequence of its highly invasive nature and of its extensive heterogeneity at the cellular and molecular levels (Bonavia et al., 2011). Recently, large-scale analysis of GBM tissues carried out within The Cancer Genome Atlas Research project confirmed an enormous complexity of epigenetic and genetic alterations in GBM. Drivers of glioblastoma initiation, progression and therapeutic resistance have been attributed to changes in the tumor metabolism and microenvironment, activation of stem pathways, epi- and genetic alterations including somatic mutations, copy number variations, genomic gains and losses and transcriptional rearrangements (Brennan et al., 2013; Schonberg et al., 2013). Currently, it seems to be very unlikely to distinguish a single genetic or cellular event that can effectively be targeted for all patients. Future therapies may require some individualization according to each patient’s tumor genotype profile or glioblastoma molecular subtype (Purow and Schiff, 2009; Verhaak et al., 2010). Recently, it has become evident that microRNAs (miRNAs) play a crucial role in the regulation of gene expression both in normal and cancer cells. miRNAs are a class of naturally occurring, small (20e23 nt) non-coding RNAs that target mRNAs at the posttranscriptional level. At sites with extensive pairing complementarity, metazoan miRNAs can direct Argonaute-catalyzed mRNA cleavage or, more frequently, direct translational repression, mRNA destabilization or a combination of the two. Animal miRNAs bind their target sites mostly in the 30 untranslated region (30 -UTR) of mRNA transcripts, but an increasing amount of evidence has confirmed many functional miRNA target sites located in protein coding regions of mammalian transcripts as well (Fang and Rajewsky, 2011; Gu et al., 2013). It is assumed that miRNAs are involved in the regulation of approximately 60% of all protein-coding genes in mammals (Friedman et al., 2009). A single miRNA can be involved in the regulation of many target RNAs and thus modulate multiple pathways. This class of non-coding RNAs is perceived as a group of crucial regulators of virtually all cellular processes, including proliferation, differentiation, apoptosis and growth. miRNAs can function either as tumor promoters (oncogenic miRNAs or ‘oncomirs’) or as tumor suppressors, and contribute to carcinogenesis. They are considered to be potentially useful targets or tools for therapeutic applications, as well as promising biomarkers of human cancers (Di Leva et al., 2014; Iorio and Croce, 2012). Alterations in the specific miRNA expression profile have been identified in a number of cancers, including gliomas. To date, several miRNA profiling studies of GBM tissues based on microarrays, PCR arrays and deep sequencing have been published (Ciafre et al., 2005; Dong et al., 2010; Hua et al., 2012; Lang et al., 2012; Rao et al., 2010; Skalsky and Cullen, 2011). A combined analysis of the mRNA and miRNA expression profiling signature has been used to identify five GBM subclasses with concordant miRNA and mRNA expression patterns corresponding to each major stage of neural stem cell differentiation (Kim et al., 2011). Integration of miRNA and mRNA signature data is an emerging potent tool for revealing subtype-specific regulators of differential expression in GBM (Setty et al., 2012). We aimed here to gain broad information on the miRNA expression pattern in malignant gliomas based on microarray
and deep sequencing approaches. We compared global miRNA expression in adult malignant gliomas (mainly glioblastoma), tumor margin and non-tumor brain tissues. Furthermore, we sought to identify miRNAs that can be associated with glioma progression by investigating the miRNA expression patterns in glioma specimens of WHO grade III and IV. The other objective was to perform a meta-analysis of miRNA profiling studies in GBM tissues. We found out that over 290 miRNAs have been reported to be deregulated in GBM tissues in comparison to normal human brain, out of which we selected a set of 35 miRNAs that are most frequently deregulated in glioma patients worldwide. Finally, we propose 30 potentially novel miRNA biomarkers of glioblastoma. To our knowledge, the meta-analysis of miRNA expression in glioblastoma has not been reported previously. Such an integrative approach is crucial for the identification of gliomaspecific miRNAs that demonstrate a potential to develop miRNA-based therapies or new diagnostic applications.
2.
Material and methods
2.1.
Tissue collection
Malignant glioma tissues and adjacent peritumoral brain tissues (glioma borders) were obtained at the time of surgery from 20 patients operated in the Department of Neurosurgery and Neurotraumatology of the Poznan University of Medical Sciences, Poland between 2010 and 2011. Tissues were flash frozen after surgery. Prior to the procedure, donors’ consent and approval from the Bioethics Council of the Poznan University of Medical Sciences had been obtained. Total RNA from normal human brain was obtained commercially (FirstChoice Human Brain Reference RNA, Ambion). According to the manufacturer’s information, this RNA sample was pooled from multiple, healthy donors and several brain regions. Additional array data on miRNA expression in three normal, postmortem brain tissues had been published previously (Hu et al., 2011) and deposited in the Gene Expression Omnibus (GEO) database. The miRNA expression array data was downloaded from GEO datasets GSM652746, GSM652747, GSM652748. It was obtained from total RNA isolated from postmortem superior frontal gyrus and analyzed using Agilent Human miRNA Microarrays (G4471A, Agilent Technologies).
2.2.
RNA isolation and quantification
Total RNA from frozen glioma and peritumoral brain tissues was extracted using a standard Trizol (Invitrogen) method, according to the manufacturer’s protocol. RNA samples were treated with DNase I using DNA-free DNase Treatment and Removal Reagent (Ambion). RNA was quantified and assessed using Agilent 2100 Bioanalyzer and RNA 6000 Nano Kit (Agilent Technologies). Only samples with RNA integrity over 8 were used for further studies.
2.3.
MicroRNA microarrays
100 ng of each RNA sample were hybridized to Agilent Human microRNA Microarrays 14.0, 8 15 K (G4471A, Agilent
Please cite this article in press as: Piwecka, M., et al., Comprehensive analysis of microRNA expression profile in malignant glioma tissues, Molecular Oncology (2015), http://dx.doi.org/10.1016/j.molonc.2015.03.007
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
MOLONC664_proof ■ 30 March 2015 ■ 3/18 M O L E C U L A R O N C O L O G Y X X X ( 2 0 1 5 ) 1 e1 8
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Technologies). MicroRNA labeling, microarray hybridization and washing, as well as data generation were performed at Warsaw University of Life Sciences, following standard Agilent’s instructions. Agilent miRNA assays used in this study contained eight 15 K microarrays on a single glass slide, each including probes for 866 human and 89 human viral miRNAs from the Sanger miRbase v14.0. The GeneSpring GX 9 software (Agilent Technologies) was used for value extraction. Prior to analysis, non-specific filtering was done to remove probes and samples with an excessive number of missing values. Probes with missing signal in at least 60% of the samples were removed, then samples in which the signal was missing in at least half of the probes were also removed. This reduced the number of probes from initial 871 to 303 and the number of samples to 20. Then the data was median-centered. The raw probe data from human gliomas (20 samples) was simultaneously analyzed with 3 datasets on miRNA expression in normal human brain (GEO GSM652746, GSM652747, GSM652748). The analysis of differential expression was performed with linear methods for microarrays using R/Bioconductor limma package. We used BenjaminieHochberg (FDR) multiple testing adjustment. Probes declared as differentially expressed realized adjusted p-values below 0.05. All original microarray data is deposited in the NCBI GEO database [GSE61710].
3
(“mirna” mode) was used. The score threshold was set to 80% of the best possible match for a given read and only the best strata of hits were reported. The resulting SAM files were processed using the in-house perl scripts in order to extract the read regions aligning to the human genome, excluding the non-aligning ends. Then, sequences were collapsed to unique representatives with the occurrence count. The process was repeated for all libraries separately.
2.4.3.
miRNA prediction and expression profiling
Sequencing has been carried out on SOLiD V4 instrument by Sequomics Biotechnology Ltd., Szeged, Hungary.
For prediction of novel miRNA species, we employed the miRDeep2 software, version 2.0.0.5 (Friedlander et al., 2012). In order to increase the predictive power of the method, we used the collapsed reads from all libraries as an input for a single analysis. The minimum read stack was set to 20. Next, the sequences of mature predicted putative miRNAs and their precursors were extracted and complemented by known miRNAs from miRBase (rel. 19). Such a complete set was used as a reference for expression profiling using the miRDeep2 quantifier.pl script. The candidates with miRDeep2 score higher than 4 or with seed region identical to known mouse miRNAs, together with a read count from all libraries higher than 100, were selected for further analyses. The miRDeep2 score cut-off was chosen based on the distribution of the calculated false positive predictions. Statistica Software, 1998 Edition, was used for the statistical analysis of miRNA differential expression. Analysis of variance (ANOVA) was used to determine the significance of the difference between means and the Tukey’s test was used for pairwise comparisons, in which p 0.05 was considered as statistically significant.
2.4.1.
2.5.
2.4.
Deep sequencing analysis
Library preparation for small RNA sequencing
Library preparation for small RNA sequencing was performed using Applied Biosystems Incorporated’s (ABI) small RNA sequencing protocol and the SOLiD V4 System. Total RNA was subjected to miRNA enrichment using PureLink miRNA Isolation Kit (Invitrogen). The samples containing small RNAs were converted into cDNA libraries using SOLiD Small RNA Expression Kit (ABI). Briefly, after hybridization and ligation to the adaptor mix, the samples were reverse transcribed to generate cDNAs. The cDNA libraries were amplified through 18 cycles of PCR using one of the supplied primer sets containing 6-nt-long sequence-specific barcodes. The individual libraries of PCR products, containing cDNAs with barcodes, were purified and size-selected in the range of 60e80 nt by electrophoresis on 10% (w/v) TBE-Urea polyacrylamide gels. The individual cDNA libraries were quantitated and titrated for emulsion PCR, which was carried out according to SOLiD 4 System Templated Bead Preparation Protocol. From the small RNA sequencing, 36-nt-long sequencing reads were obtained.
2.4.2.
Pre-processing and mapping of the sequencing reads
From the raw sequencing reads, 30 adapter sequences were removed using the cutadapt software [http://journal.embnet.org/index.php/embnetjournal/article/view/200]. Next, cleaned reads were aligned to the human genome (hg19 assembly) in color space using the SHRIMP2 software, version 2.2.3 (David et al., 2011). The alignment mode was set to local and the default set of options optimized for miRNA mapping
Functional analysis with IPA
The selected miRNAs were analyzed to identify the networks and pathways. For this purpose, we used Ingenuity Pathway Analysis software (IPA, Ingenuity Systems; http://www.ingenuity.com). A core analysis was employed to identify the most relevant miRNA targets, canonical pathways, biological functions and physiological processes from the interactions reported in the IPA database.
2.6. Meta-analysis: identification and selection of relevant studies PubMed was used to search for glioblastoma/gliomas miRNA expression profiling studies published up to January 2014. The database was searched by means of the MeSH terms: ‘glioblastoma’, ‘GBM’, ‘gliomas’ and ‘malignant gliomas’ in combination with ‘miRNA’ and ‘microRNA’. Potentially associated publications were assessed by looking through their abstracts and the most relevant publications were subjected to closer examination. Eligible studies had to meet the following criteria: (i) they were miRNA expression profiling studies in malignant glioma/glioblastoma patients, (ii) they used tissue samples obtained from surgically resected glial tumors and corresponding noncancerous or normal tissues for comparison and differential expression analysis, and (iii) the validation method and validation sample set were reported. Accordingly, studies were excluded based on the following
Please cite this article in press as: Piwecka, M., et al., Comprehensive analysis of microRNA expression profile in malignant glioma tissues, Molecular Oncology (2015), http://dx.doi.org/10.1016/j.molonc.2015.03.007
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
MOLONC664_proof ■ 30 March 2015 ■ 4/18
4 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
M O L E C U L A R O N C O L O G Y X X X ( 2 0 1 5 ) 1 e1 8
criteria: (i) they were not conducted on human samples, (ii) they lacked normal or noncancerous tissue as a reference in comparison analysis, and (iii) they were review articles or profiling studies in the serum, cerebrospinal fluid or glioma cell lines.
2.6.1.
Data extraction
From the full text and corresponding supplement information, the following items were collected for each study: author, journal and year of publication, location of the study, selection and characteristics of tumor and non-tumor samples used in the study, platform of miRNA profiling, cut-off criteria for differentially expressed (DE) miRNAs and for the list of upand downregulated miRNA features and their corresponding fold change, if available. The extraction was performed by two reviewers independently. The differentially expressed miRNAs were ranked according to the number of the studies that consistently reported the miRNA as DE and with a consistent direction of change.
borders of tumors in comparison to normal brain, but did not turn out to be significantly deregulated in tumors. To test whether the discovered glioblastoma-specific upand down-regulated miRNAs were known to be associated with other human pathologies, we used the miR2Disease (Jiang et al., 2009) and miRBase databases (Supporting Information Tables S2 and S3). An inspection of the current version of miRBase showed that some entries conflict with other ncRNAs, in particular snoRNAs and tRNA. When it came to our analysis, miR-720, miR-1274a, miR-1308, miR1826 and miR-1280 were recognized as misannotated. Mir2Disease inspection revealed that miR-155, miR-21, miR-210 and miR-93 were the most widely identified as up-regulated in human cancers. It is interesting to note that from the group of miRNAs down-regulated in GBM, miR-381, miR-379 and miR299-5p, appeared to be consistently up-regulated in a spectrum of primary muscular disorders, including Duchenne muscular dystrophy, myopathies, nemaline myopathy and dermatomyositis.
3.2.
3.
Results
3.1. MicroRNA expression pattern in malignant glioma tissues Using miRNA microarrays, we investigated the miRNA expression pattern in malignant glioma tissues, tumor margins (adjacent peritumoral brain tissues) and normal human brain. For the latter, we used commercially available brain total RNA from multiple, healthy donors and postmortem nonneoplastic brains. A brief characterization of tumor samples is presented in Supporting Information Table S1. By differential expression (DE) analysis with normal human brain, we sought to determine the miRNA signature both in glioblastoma tissues and in peritumoral tissues. Specifically, normal-versus-peritumoral comparison was conducted to identify those miRNAs that could be associated with glioblastoma tumor cell infiltration. DE analysis revealed distinct expression patterns of miRNAs in tumors and normal brain samples, whereas borders of tumors did not constitute a separate cluster (Figure 1A and B). We observed a significant change in the expression of 97 miRNAs derived from tumor tissues in comparison to normal human brain, including 56 up-regulated and 41 downregulated miRNAs, respectively (Supporting Information Tables S2 and S3 and Figure 1A). The comparison of adjacent peritumoral tissues with normal brain samples revealed 6 miRNAs with elevated expression and 16 miRNAs downregulated within the borders of tumors (Supporting Information Tables S4 and S5 and Figure 1B). Twenty one out of 22 miRNAs identified in borders of gliomas were also found in GBMs (Figure 1C). The group of miRNAs that were similarly deregulated both in tumor and peritumoral tissues included: up-regulated miR-21, miR-630, miR-155, miR-1260, miR-5425p, miR-142-5p and down-regulated miR-137, miR-124, miR129-3p, miR-769-5p, miR-132, miR-128, miR-7, miR-410, miR136, miR-153, miR-323-3p, miR-330-3p, miR-598, among others. Only miR-625 was exclusively identified within the
Small RNA deep sequencing
Since miRNA pattern may change during malignant glioma progression, we sought to test the miRNA profiles at different stages of glioma malignancy. Therefore, we used SOLiD sequencingbased miRNA expression profiling to identify the most abundant miRNAs in malignant glioma of WHO grade III and IV (glioblastoma multiforme) and we performed differential expression analysis. Totally, over 19,257 million effective reads were obtained from 4 libraries analyzed (Supporting Information Table S6). After filtration to eliminate rRNA, tRNA, snRNA and snoRNA sequences, the remaining effective reads were mapped to miRBase release 19. The stringency of miRNA detection was lowered to cut-off of 100 counts across the four libraries. The matches corresponded to 489 of the known miRNAs deposited in miRBase release 19 (Supporting Information Table S7). The 20 most abundant miRNAs that were identified in our study are shown in Table 1. They represent as much as 61.7% of all sequence reads and over 76% of sequence reads mapped to known, previously described miRNAs. After identifying the known miRNAs, we found out that over 25,000 reads (0.13%) corresponded to 25 potentially novel miRNAs (Supporting Information Table S8). Eighteen top scoring miRNAs possessed miRDeep2 score over 5 and the estimated probability that the miRNA candidate was a true positive accounted for 64 14%. Another 7 miRNA candidates with miRDeep2 score below 5 were predicted as potentially novel human miRNAs based on their high sequence similarity to known miRNAs from mice. Differentially expressed miRNAs were determined by fold change (FC) and log2 fold change in normalized read counts between WHO grade III and grade IV gliomas. DE analysis revealed 25 miRNAs that were differentially expressed (p-value 1.5 or FC < 0.5
29 10
29 5
ND 5
FC
71
43
28
FC
14
12
2
ND
17
9
8
p-value
97
30
67
FC
94
55
29
FC
8
2
6
p-value 0.01
10
7
3
p-value 5
56
18
38
Cancer/Normal ratio Self-made (arrays) AND Taqman MicroRNA Assay (Applied Biosystems) Agilent Technologies ND
26
15
11
22
15
miRNA microarray Agilent Technologies (TCGA data) miRNA microarray Mercury LNA array version 10.0 (Exiqon) Real-time PCR TaqMan MicroRNA Assay (Applied Biosystems) miRNA microarray Agilent Technologies (TCGA data) Deep sequencing Illumina Membrane- array hybridization AND real-time PCR miRNA microarray
Total Up-regulated Down-regulated miRNAs miRNAs
7 (continued on next page)
MOLONC664_proof ■ 30 March 2015 ■ 9/18
Ciafre
Method
M O L E C U L A R O N C O L O G Y X X X ( 2 0 1 5 ) 1 e1 8
9
Please cite this article in press as: Piwecka, M., et al., Comprehensive analysis of microRNA expression profile in malignant glioma tissues, Molecular Oncology (2015), http://dx.doi.org/10.1016/j.molonc.2015.03.007
a Table 3 e MicroRNA expression profiling studies of human malignant glioma tissues included in the meta-analysis.
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
MOLONC664_proof ■ 30 March 2015 ■ 10/18
10
40 73 Deep sequencing
Illumina
p-value
available at www.sciencedirect.com 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
ScienceDirect www.elsevier.com/locate/molonc
Q5
Q4
Comprehensive analysis of microRNA expression profile in malignant glioma tissues Monika Piweckaa,1, Katarzyna Rollea,1, Agnieszka Beltera, c _ Anna Maria Barciszewskab, Marek Zywicki , Marcin Michalakd, Stanis1aw Nowakb, Miros1awa Z. Naskre˛t-Barciszewskaa, Jan Barciszewskia,* a
Institute of Bioorganic Chemistry, Polish Academy of Sciences, Poznan, Poland Department of Neurosurgery and Neurotraumatology, Poznan University of Medical Sciences, Poznan, Poland c Bioinformatics Laboratory, Institute of Molecular Biology and Biotechnology, Faculty of Biology, A. Mickiewicz University, Poznan, Poland d Institute of Dendrology, Polish Academy of Sciences, Kornik, Poland b
A R T I C L E
I N F O
A B S T R A C T
Article history:
Malignant gliomas represent the most devastating group of brain tumors in adults, among
Received 16 October 2014
which glioblastoma multiforme (GBM) exhibits the highest malignancy rate. Despite com-
Received in revised form
bined modality treatment, GBM recurs and is invariably fatal. A further insight into the mo-
26 February 2015
lecular background of gliomagenesis is required to improve patient outcomes. The primary
Accepted 16 March 2015
aim of this study was to gain broad information on the miRNA expression pattern in ma-
Available online -
lignant gliomas, mainly GBM. We investigated the global miRNA profile of malignant glioma tissues with miRNA microarrays, deep sequencing and meta-analysis. We selected
Keywords:
miRNAs that were most frequently deregulated in glioblastoma tissues, as well as in peri-
Brain tumors
tumoral areas, in comparison with normal human brain. We identified candidate miRNAs
Glioblastoma
associated with the progression from glioma grade III to glioma grade IV. The meta-
Gliomas
analysis of miRNA profiling studies in GBM tissues summarizes the past and recent ad-
Meta-Analysis
vances in the investigation of the miRNA signature in GBM versus noncancerous human
microRNA
brain and provides a comprehensive overview. We propose a list of 35 miRNAs whose
miRNA profiling
expression is most frequently deregulated in GBM patients and of 30 miRNA candidates recognized as novel GBM biomarkers. ª 2015 Published by Elsevier B.V. on behalf of Federation of European Biochemical Societies.
1.
Introduction
Malignant gliomas are the most common and aggressive primary brain tumors in adults. WHO grade IV glioblastoma multiforme (GBM) is the most frequent type and is characterized
by strong vascular proliferation, invasiveness, diminished apoptosis, radio- and chemoresistance. Despite treatment, the median survival of GBM patients is approximately 12 months after diagnosis (Cheng et al., 2010; Stupp et al., 2005).
* Corresponding author. Tel.: þ48 61 8528503x132; fax: þ48 61 8528919. E-mail address: [email protected] (J. Barciszewski). 1 These authors contributed equally to this work. http://dx.doi.org/10.1016/j.molonc.2015.03.007 1574-7891/ª 2015 Published by Elsevier B.V. on behalf of Federation of European Biochemical Societies.
Please cite this article in press as: Piwecka, M., et al., Comprehensive analysis of microRNA expression profile in malignant glioma tissues, Molecular Oncology (2015), http://dx.doi.org/10.1016/j.molonc.2015.03.007
Q1
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
MOLONC664_proof ■ 30 March 2015 ■ 2/18
2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
M O L E C U L A R O N C O L O G Y X X X ( 2 0 1 5 ) 1 e1 8
Treatment of GBM remains one of the most challenging tasks in clinical oncology. It is the consequence of its highly invasive nature and of its extensive heterogeneity at the cellular and molecular levels (Bonavia et al., 2011). Recently, large-scale analysis of GBM tissues carried out within The Cancer Genome Atlas Research project confirmed an enormous complexity of epigenetic and genetic alterations in GBM. Drivers of glioblastoma initiation, progression and therapeutic resistance have been attributed to changes in the tumor metabolism and microenvironment, activation of stem pathways, epi- and genetic alterations including somatic mutations, copy number variations, genomic gains and losses and transcriptional rearrangements (Brennan et al., 2013; Schonberg et al., 2013). Currently, it seems to be very unlikely to distinguish a single genetic or cellular event that can effectively be targeted for all patients. Future therapies may require some individualization according to each patient’s tumor genotype profile or glioblastoma molecular subtype (Purow and Schiff, 2009; Verhaak et al., 2010). Recently, it has become evident that microRNAs (miRNAs) play a crucial role in the regulation of gene expression both in normal and cancer cells. miRNAs are a class of naturally occurring, small (20e23 nt) non-coding RNAs that target mRNAs at the posttranscriptional level. At sites with extensive pairing complementarity, metazoan miRNAs can direct Argonaute-catalyzed mRNA cleavage or, more frequently, direct translational repression, mRNA destabilization or a combination of the two. Animal miRNAs bind their target sites mostly in the 30 untranslated region (30 -UTR) of mRNA transcripts, but an increasing amount of evidence has confirmed many functional miRNA target sites located in protein coding regions of mammalian transcripts as well (Fang and Rajewsky, 2011; Gu et al., 2013). It is assumed that miRNAs are involved in the regulation of approximately 60% of all protein-coding genes in mammals (Friedman et al., 2009). A single miRNA can be involved in the regulation of many target RNAs and thus modulate multiple pathways. This class of non-coding RNAs is perceived as a group of crucial regulators of virtually all cellular processes, including proliferation, differentiation, apoptosis and growth. miRNAs can function either as tumor promoters (oncogenic miRNAs or ‘oncomirs’) or as tumor suppressors, and contribute to carcinogenesis. They are considered to be potentially useful targets or tools for therapeutic applications, as well as promising biomarkers of human cancers (Di Leva et al., 2014; Iorio and Croce, 2012). Alterations in the specific miRNA expression profile have been identified in a number of cancers, including gliomas. To date, several miRNA profiling studies of GBM tissues based on microarrays, PCR arrays and deep sequencing have been published (Ciafre et al., 2005; Dong et al., 2010; Hua et al., 2012; Lang et al., 2012; Rao et al., 2010; Skalsky and Cullen, 2011). A combined analysis of the mRNA and miRNA expression profiling signature has been used to identify five GBM subclasses with concordant miRNA and mRNA expression patterns corresponding to each major stage of neural stem cell differentiation (Kim et al., 2011). Integration of miRNA and mRNA signature data is an emerging potent tool for revealing subtype-specific regulators of differential expression in GBM (Setty et al., 2012). We aimed here to gain broad information on the miRNA expression pattern in malignant gliomas based on microarray
and deep sequencing approaches. We compared global miRNA expression in adult malignant gliomas (mainly glioblastoma), tumor margin and non-tumor brain tissues. Furthermore, we sought to identify miRNAs that can be associated with glioma progression by investigating the miRNA expression patterns in glioma specimens of WHO grade III and IV. The other objective was to perform a meta-analysis of miRNA profiling studies in GBM tissues. We found out that over 290 miRNAs have been reported to be deregulated in GBM tissues in comparison to normal human brain, out of which we selected a set of 35 miRNAs that are most frequently deregulated in glioma patients worldwide. Finally, we propose 30 potentially novel miRNA biomarkers of glioblastoma. To our knowledge, the meta-analysis of miRNA expression in glioblastoma has not been reported previously. Such an integrative approach is crucial for the identification of gliomaspecific miRNAs that demonstrate a potential to develop miRNA-based therapies or new diagnostic applications.
2.
Material and methods
2.1.
Tissue collection
Malignant glioma tissues and adjacent peritumoral brain tissues (glioma borders) were obtained at the time of surgery from 20 patients operated in the Department of Neurosurgery and Neurotraumatology of the Poznan University of Medical Sciences, Poland between 2010 and 2011. Tissues were flash frozen after surgery. Prior to the procedure, donors’ consent and approval from the Bioethics Council of the Poznan University of Medical Sciences had been obtained. Total RNA from normal human brain was obtained commercially (FirstChoice Human Brain Reference RNA, Ambion). According to the manufacturer’s information, this RNA sample was pooled from multiple, healthy donors and several brain regions. Additional array data on miRNA expression in three normal, postmortem brain tissues had been published previously (Hu et al., 2011) and deposited in the Gene Expression Omnibus (GEO) database. The miRNA expression array data was downloaded from GEO datasets GSM652746, GSM652747, GSM652748. It was obtained from total RNA isolated from postmortem superior frontal gyrus and analyzed using Agilent Human miRNA Microarrays (G4471A, Agilent Technologies).
2.2.
RNA isolation and quantification
Total RNA from frozen glioma and peritumoral brain tissues was extracted using a standard Trizol (Invitrogen) method, according to the manufacturer’s protocol. RNA samples were treated with DNase I using DNA-free DNase Treatment and Removal Reagent (Ambion). RNA was quantified and assessed using Agilent 2100 Bioanalyzer and RNA 6000 Nano Kit (Agilent Technologies). Only samples with RNA integrity over 8 were used for further studies.
2.3.
MicroRNA microarrays
100 ng of each RNA sample were hybridized to Agilent Human microRNA Microarrays 14.0, 8 15 K (G4471A, Agilent
Please cite this article in press as: Piwecka, M., et al., Comprehensive analysis of microRNA expression profile in malignant glioma tissues, Molecular Oncology (2015), http://dx.doi.org/10.1016/j.molonc.2015.03.007
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
MOLONC664_proof ■ 30 March 2015 ■ 3/18 M O L E C U L A R O N C O L O G Y X X X ( 2 0 1 5 ) 1 e1 8
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Technologies). MicroRNA labeling, microarray hybridization and washing, as well as data generation were performed at Warsaw University of Life Sciences, following standard Agilent’s instructions. Agilent miRNA assays used in this study contained eight 15 K microarrays on a single glass slide, each including probes for 866 human and 89 human viral miRNAs from the Sanger miRbase v14.0. The GeneSpring GX 9 software (Agilent Technologies) was used for value extraction. Prior to analysis, non-specific filtering was done to remove probes and samples with an excessive number of missing values. Probes with missing signal in at least 60% of the samples were removed, then samples in which the signal was missing in at least half of the probes were also removed. This reduced the number of probes from initial 871 to 303 and the number of samples to 20. Then the data was median-centered. The raw probe data from human gliomas (20 samples) was simultaneously analyzed with 3 datasets on miRNA expression in normal human brain (GEO GSM652746, GSM652747, GSM652748). The analysis of differential expression was performed with linear methods for microarrays using R/Bioconductor limma package. We used BenjaminieHochberg (FDR) multiple testing adjustment. Probes declared as differentially expressed realized adjusted p-values below 0.05. All original microarray data is deposited in the NCBI GEO database [GSE61710].
3
(“mirna” mode) was used. The score threshold was set to 80% of the best possible match for a given read and only the best strata of hits were reported. The resulting SAM files were processed using the in-house perl scripts in order to extract the read regions aligning to the human genome, excluding the non-aligning ends. Then, sequences were collapsed to unique representatives with the occurrence count. The process was repeated for all libraries separately.
2.4.3.
miRNA prediction and expression profiling
Sequencing has been carried out on SOLiD V4 instrument by Sequomics Biotechnology Ltd., Szeged, Hungary.
For prediction of novel miRNA species, we employed the miRDeep2 software, version 2.0.0.5 (Friedlander et al., 2012). In order to increase the predictive power of the method, we used the collapsed reads from all libraries as an input for a single analysis. The minimum read stack was set to 20. Next, the sequences of mature predicted putative miRNAs and their precursors were extracted and complemented by known miRNAs from miRBase (rel. 19). Such a complete set was used as a reference for expression profiling using the miRDeep2 quantifier.pl script. The candidates with miRDeep2 score higher than 4 or with seed region identical to known mouse miRNAs, together with a read count from all libraries higher than 100, were selected for further analyses. The miRDeep2 score cut-off was chosen based on the distribution of the calculated false positive predictions. Statistica Software, 1998 Edition, was used for the statistical analysis of miRNA differential expression. Analysis of variance (ANOVA) was used to determine the significance of the difference between means and the Tukey’s test was used for pairwise comparisons, in which p 0.05 was considered as statistically significant.
2.4.1.
2.5.
2.4.
Deep sequencing analysis
Library preparation for small RNA sequencing
Library preparation for small RNA sequencing was performed using Applied Biosystems Incorporated’s (ABI) small RNA sequencing protocol and the SOLiD V4 System. Total RNA was subjected to miRNA enrichment using PureLink miRNA Isolation Kit (Invitrogen). The samples containing small RNAs were converted into cDNA libraries using SOLiD Small RNA Expression Kit (ABI). Briefly, after hybridization and ligation to the adaptor mix, the samples were reverse transcribed to generate cDNAs. The cDNA libraries were amplified through 18 cycles of PCR using one of the supplied primer sets containing 6-nt-long sequence-specific barcodes. The individual libraries of PCR products, containing cDNAs with barcodes, were purified and size-selected in the range of 60e80 nt by electrophoresis on 10% (w/v) TBE-Urea polyacrylamide gels. The individual cDNA libraries were quantitated and titrated for emulsion PCR, which was carried out according to SOLiD 4 System Templated Bead Preparation Protocol. From the small RNA sequencing, 36-nt-long sequencing reads were obtained.
2.4.2.
Pre-processing and mapping of the sequencing reads
From the raw sequencing reads, 30 adapter sequences were removed using the cutadapt software [http://journal.embnet.org/index.php/embnetjournal/article/view/200]. Next, cleaned reads were aligned to the human genome (hg19 assembly) in color space using the SHRIMP2 software, version 2.2.3 (David et al., 2011). The alignment mode was set to local and the default set of options optimized for miRNA mapping
Functional analysis with IPA
The selected miRNAs were analyzed to identify the networks and pathways. For this purpose, we used Ingenuity Pathway Analysis software (IPA, Ingenuity Systems; http://www.ingenuity.com). A core analysis was employed to identify the most relevant miRNA targets, canonical pathways, biological functions and physiological processes from the interactions reported in the IPA database.
2.6. Meta-analysis: identification and selection of relevant studies PubMed was used to search for glioblastoma/gliomas miRNA expression profiling studies published up to January 2014. The database was searched by means of the MeSH terms: ‘glioblastoma’, ‘GBM’, ‘gliomas’ and ‘malignant gliomas’ in combination with ‘miRNA’ and ‘microRNA’. Potentially associated publications were assessed by looking through their abstracts and the most relevant publications were subjected to closer examination. Eligible studies had to meet the following criteria: (i) they were miRNA expression profiling studies in malignant glioma/glioblastoma patients, (ii) they used tissue samples obtained from surgically resected glial tumors and corresponding noncancerous or normal tissues for comparison and differential expression analysis, and (iii) the validation method and validation sample set were reported. Accordingly, studies were excluded based on the following
Please cite this article in press as: Piwecka, M., et al., Comprehensive analysis of microRNA expression profile in malignant glioma tissues, Molecular Oncology (2015), http://dx.doi.org/10.1016/j.molonc.2015.03.007
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
MOLONC664_proof ■ 30 March 2015 ■ 4/18
4 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
M O L E C U L A R O N C O L O G Y X X X ( 2 0 1 5 ) 1 e1 8
criteria: (i) they were not conducted on human samples, (ii) they lacked normal or noncancerous tissue as a reference in comparison analysis, and (iii) they were review articles or profiling studies in the serum, cerebrospinal fluid or glioma cell lines.
2.6.1.
Data extraction
From the full text and corresponding supplement information, the following items were collected for each study: author, journal and year of publication, location of the study, selection and characteristics of tumor and non-tumor samples used in the study, platform of miRNA profiling, cut-off criteria for differentially expressed (DE) miRNAs and for the list of upand downregulated miRNA features and their corresponding fold change, if available. The extraction was performed by two reviewers independently. The differentially expressed miRNAs were ranked according to the number of the studies that consistently reported the miRNA as DE and with a consistent direction of change.
borders of tumors in comparison to normal brain, but did not turn out to be significantly deregulated in tumors. To test whether the discovered glioblastoma-specific upand down-regulated miRNAs were known to be associated with other human pathologies, we used the miR2Disease (Jiang et al., 2009) and miRBase databases (Supporting Information Tables S2 and S3). An inspection of the current version of miRBase showed that some entries conflict with other ncRNAs, in particular snoRNAs and tRNA. When it came to our analysis, miR-720, miR-1274a, miR-1308, miR1826 and miR-1280 were recognized as misannotated. Mir2Disease inspection revealed that miR-155, miR-21, miR-210 and miR-93 were the most widely identified as up-regulated in human cancers. It is interesting to note that from the group of miRNAs down-regulated in GBM, miR-381, miR-379 and miR299-5p, appeared to be consistently up-regulated in a spectrum of primary muscular disorders, including Duchenne muscular dystrophy, myopathies, nemaline myopathy and dermatomyositis.
3.2.
3.
Results
3.1. MicroRNA expression pattern in malignant glioma tissues Using miRNA microarrays, we investigated the miRNA expression pattern in malignant glioma tissues, tumor margins (adjacent peritumoral brain tissues) and normal human brain. For the latter, we used commercially available brain total RNA from multiple, healthy donors and postmortem nonneoplastic brains. A brief characterization of tumor samples is presented in Supporting Information Table S1. By differential expression (DE) analysis with normal human brain, we sought to determine the miRNA signature both in glioblastoma tissues and in peritumoral tissues. Specifically, normal-versus-peritumoral comparison was conducted to identify those miRNAs that could be associated with glioblastoma tumor cell infiltration. DE analysis revealed distinct expression patterns of miRNAs in tumors and normal brain samples, whereas borders of tumors did not constitute a separate cluster (Figure 1A and B). We observed a significant change in the expression of 97 miRNAs derived from tumor tissues in comparison to normal human brain, including 56 up-regulated and 41 downregulated miRNAs, respectively (Supporting Information Tables S2 and S3 and Figure 1A). The comparison of adjacent peritumoral tissues with normal brain samples revealed 6 miRNAs with elevated expression and 16 miRNAs downregulated within the borders of tumors (Supporting Information Tables S4 and S5 and Figure 1B). Twenty one out of 22 miRNAs identified in borders of gliomas were also found in GBMs (Figure 1C). The group of miRNAs that were similarly deregulated both in tumor and peritumoral tissues included: up-regulated miR-21, miR-630, miR-155, miR-1260, miR-5425p, miR-142-5p and down-regulated miR-137, miR-124, miR129-3p, miR-769-5p, miR-132, miR-128, miR-7, miR-410, miR136, miR-153, miR-323-3p, miR-330-3p, miR-598, among others. Only miR-625 was exclusively identified within the
Small RNA deep sequencing
Since miRNA pattern may change during malignant glioma progression, we sought to test the miRNA profiles at different stages of glioma malignancy. Therefore, we used SOLiD sequencingbased miRNA expression profiling to identify the most abundant miRNAs in malignant glioma of WHO grade III and IV (glioblastoma multiforme) and we performed differential expression analysis. Totally, over 19,257 million effective reads were obtained from 4 libraries analyzed (Supporting Information Table S6). After filtration to eliminate rRNA, tRNA, snRNA and snoRNA sequences, the remaining effective reads were mapped to miRBase release 19. The stringency of miRNA detection was lowered to cut-off of 100 counts across the four libraries. The matches corresponded to 489 of the known miRNAs deposited in miRBase release 19 (Supporting Information Table S7). The 20 most abundant miRNAs that were identified in our study are shown in Table 1. They represent as much as 61.7% of all sequence reads and over 76% of sequence reads mapped to known, previously described miRNAs. After identifying the known miRNAs, we found out that over 25,000 reads (0.13%) corresponded to 25 potentially novel miRNAs (Supporting Information Table S8). Eighteen top scoring miRNAs possessed miRDeep2 score over 5 and the estimated probability that the miRNA candidate was a true positive accounted for 64 14%. Another 7 miRNA candidates with miRDeep2 score below 5 were predicted as potentially novel human miRNAs based on their high sequence similarity to known miRNAs from mice. Differentially expressed miRNAs were determined by fold change (FC) and log2 fold change in normalized read counts between WHO grade III and grade IV gliomas. DE analysis revealed 25 miRNAs that were differentially expressed (p-value 1.5 or FC < 0.5
29 10
29 5
ND 5
FC
71
43
28
FC
14
12
2
ND
17
9
8
p-value
97
30
67
FC
94
55
29
FC
8
2
6
p-value 0.01
10
7
3
p-value 5
56
18
38
Cancer/Normal ratio Self-made (arrays) AND Taqman MicroRNA Assay (Applied Biosystems) Agilent Technologies ND
26
15
11
22
15
miRNA microarray Agilent Technologies (TCGA data) miRNA microarray Mercury LNA array version 10.0 (Exiqon) Real-time PCR TaqMan MicroRNA Assay (Applied Biosystems) miRNA microarray Agilent Technologies (TCGA data) Deep sequencing Illumina Membrane- array hybridization AND real-time PCR miRNA microarray
Total Up-regulated Down-regulated miRNAs miRNAs
7 (continued on next page)
MOLONC664_proof ■ 30 March 2015 ■ 9/18
Ciafre
Method
M O L E C U L A R O N C O L O G Y X X X ( 2 0 1 5 ) 1 e1 8
9
Please cite this article in press as: Piwecka, M., et al., Comprehensive analysis of microRNA expression profile in malignant glioma tissues, Molecular Oncology (2015), http://dx.doi.org/10.1016/j.molonc.2015.03.007
a Table 3 e MicroRNA expression profiling studies of human malignant glioma tissues included in the meta-analysis.
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
MOLONC664_proof ■ 30 March 2015 ■ 10/18
10
40 73 Deep sequencing
Illumina
p-value
