Abstract
Myung-Chul Kim1,2, Na-Yon Kim1,2, Yu-Ri Seo1 and Yongbaek Kim1,3
1Laboratory of Clinical Pathology, College of Veterinary Medicine, Seoul National University, Gwanak-Gu, Seoul 151-742, The Republic of Korea
2BK21 PLUS Program for Creative Veterinary Science Research, College of Veterinary Medicine, Seoul National University, Gwanak-Gu, Seoul 151-742, The Republic of Korea
3Research Institute for Veterinary Science, College of Veterinary Medicine, Seoul National University, Gwanak-Gu, Seoul 151-742, The Republic of Korea
Correspondence to:
Yongbaek Kim, email: [email protected]
Keywords: mesothelioma, side population, microRNA, intratumoral heterogeneity, microarray
Received: December 08, 2016 Accepted: March 22, 2017 Published: April 13, 2017
ABSTRACT
This study was performed to investigate the global expression profile of microRNAs in distinct subpopulations of a human malignant mesothelioma cell line. Total RNAs were isolated from the sorted side population and non-side population of MS1. The RNAs were subjected to analysis using Affymetrix GeneChip microRNA Arrays. After data extraction and normalization, a subset of microRNAs defining cell subpopulations was identified using bioinformatics softwares. Based on the criteria of 2-fold difference and the p-value of < 0.05, a total of 95 microRNAs were differentially expressed in the side population compared to the non-side population. Functional ontology revealed that target genes of the miRNAs were categorized into various gene ontology terms, such as stem cell maintenance, cell proliferation, programmed cell death, cell migration, and cellular response to stress. The Kyoto Encyclopedia of Genes and Genomes analysis showed that ErbB-2 receptor tyrosine kinases signaling pathway was the most represented. Integrated analysis of MiRTarBase and RNA-seq identified 12 target genes of microRNAs defining side population, including DDIT4 and ROCK2. The present study indicates that a distinct set of microRNAs may be critically involved in the generation and maintenance of heterogeneous subpopulations of cancer cells. They could be a plausible target for the eradication of more aggressive cancer cell subpopulations.
BACKGROUND
Harboring heterogeneous cell populations in a tumor mass is a common feature of many cancer types [1]. The intratumoral heterogeneity refers to a hierarchical organization of phenotypically and functionally distinct cancer cells within a tumor [2]. Cancer cell subpopulations with distinct biological properties are considered to be the main obstacle against effective cancer therapy [3]. Clonal evolution has been proposed to explain the development of intratumoral heterogeneity [1]. Recently, it has been reported that epigenetic modifications contribute to the genesis of tumor cell heterogeneity [4–6].
MicroRNAs (miRNAs) are a class of non-coding small RNAs that post-transcriptionally regulate gene expression. They pair with the 3’-untranslated regions or the open reading frames of their target mRNAs. The interaction between miRNA and target gene leads to either degradation of the target genes or inhibition of the protein translation [7]. Dysregulation of miRNAs has been implicated in cancer progression by acting as either oncogenes or tumor suppressor genes [8].
Human malignant mesothelioma (HMM) is an aggressive cancer arising from the surface of body cavities [9]. Many factors including environmental contaminants and viruses have been implicated in the development of HMM [10]. Although significant progress has been made in terms of etiology and pathogenesis, the prognosis of HMM patients remains dismal with the survival time of less than 12 months from the diagnosis [9, 10]. HMMs are markedly heterogeneous in morphology as well as in biology, which makes HMM a paradigmatic model for the intratumoral heterogeneity [11]. Formation of three different histopathologic subtypes from single cell in vitro indicated the potency of HMM cells in generating tumor heterogeneity [11–13]. Using side population (SP) assay, subpopulation cells harboring stem cell properties were identified in HMM cell lines [13].
In an attempt to elucidate the involvement of epigenetics in the generation and maintenance of intratumoral heterogeneity, the miRNA expression profiles of HMM cell subpopulations were investigated using microarray analysis. A distinct subset of the miRNAs was identified from cancer cell subpopulations, and potential signaling pathways regulated by these miRNAs were determined. The present study provides background information to advance our understanding about the regulatory function of miRNAs in the generation of intratumoral heterogeneity.
RESULTS
Isolation of RNA from the sorted SP and NSP cells of MS1
SP assay composed of Hoechst 33342 dye staining and subsequent flow cytometry illustrated a distinct SP cells in the MS1 cell line as a tail in flow cytometry plot (Figure 1A). The SP fraction was significantly decreased by the treatment of 50 μM verapamil hydrochloride (Figure 1B). The isolated RNAs were determined to be of good quality with no degradation by A260/280 ratio greater than 1.8 determined using Agilent’s 2100 Bioanalyzer and the RNA Integrity Number (RIN) value higher than 8 measured using Nanodrop.
Figure 1: Identification of side population (SP) cells in the MS1 cell line. (A) Side population assay revealed that the MS1 cell line contained a distinct region of SP cells indicated by a trapezoid on each panel. (B) Treatment with verapamil hydrochloride significantly reduced the fraction of SP cells.
Identification of differentially expressed miRNAs (DEMs) by microarray analysis
Based on the criteria of over 2-fold difference and p-value less than 0.05, a total of 95 DEMs were identified to be differentially expressed in SP cells compared to NSP cells. Among the 95 DEMs, 42 DEMs were significantly up-regulated and 53 DEMs were significantly down-regulated in the SP cells compared to the NSP cells (Supplementary Table 1). Top 10 up-regulated and down-regulated miRNAs were presented, respectively (Table 1) The microarray data are available at the National Center for Biotechnology Information Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) under the accession number of GSE69910.
Table 1: Top 10 up- and down-regulated miRNAs defining SP cells identified by miRNA microarray
microRNA | Fold change (log2 ratio) | p-value |
|---|---|---|
hsa-mir-3198-1 | 4.129 | 2.19E-10 |
hsa-mir-3198-2 | 4.129 | 2.19E-10 |
hsa-mir-4497 | 2.666 | 8.72E-16 |
hsa-mir-138-1 | 2.420 | 2.90E-15 |
hsa-mir-4304 | 2.413 | 1.71E-08 |
hsa-mir-1281 | 1.882 | 3.12E-08 |
hsa-mir-489 | 1.859 | 4.83E-07 |
hsa-mir-4745 | 1.850 | 1.63E-14 |
hsa-mir-301a | 1.782 | 1.54E-09 |
hsa-mir-3935 | 1.642 | 9.57E-09 |
hsa-mir-148b | -2.231 | 1.42E-06 |
hsa-mir-484 | -2.288 | 0.000582039 |
hsa-mir-584 | -2.290 | 0.007460966 |
hsa-mir-425 | -2.290 | 0.0332853 |
hsa-mir-197 | -2.432 | 4.96E-07 |
hsa-mir-629 | -2.646 | 0.01597814 |
hsa-mir-183 | -2.835 | 0.01547796 |
hsa-mir-4485 | -2.884 | 1.43E-07 |
hsa-mir-4443 | -3.221 | 2.84E-06 |
hsa-mir-1246 | -4.673 | 0.000618761 |
A total of 20 miRNAs that were most significantly altered defining SP cells in MS1 cell line were presented.
GO and pathway analyses of the predicted target genes of DEMs
The identification of target genes for miRNAs was performed by several computational algorithms. The TargetScan database was used to predict target genes of DEMs. A total of 1,743 target genes generated by the target prediction software were subjected to the GO analysis in order to determine key-regulatory components and functional relationships of the predicted target genes [14]. The predicted target genes were categorized into biological processes, molecular functions, and cellular components. GO analysis revealed that 287 GO terms were involved in the domain of biological processes including 57 GO terms in the domain of molecular functions and 52 GO terms in the domain of cellular components. To minimize redundancy in the lists of enriched GO categories, the GO terms were merged and replaced with their representative subset of terms based on semantic similarity measures using REVIGO. The initial list of GO terms of biological processes was reduced to 123, eliminating 164 largely redundant terms. The GO terms of 57 molecular functions and 52 cellular components were reduced to 48 and 39 non-redundant terms, respectively. The non-redundant GO terms with higher than 1% frequency were visualized based on their semantic similarities in a semantic space (Figure 2). The cluster representatives of biological processes included regulation of cell proliferation, negative regulation of gene expression, regulation of cell migration, regulation of cellular response to stress, regulation of apoptotic process, and regulation of cell communication. The detailed information about the non-redundant GO terms including frequency, EASE score, and uniqueness, was presented (Supplementary Table 2).
Figure 2: Gene Ontology scatterplot constructed by REVIGO. GO terms enriched in predicted target genes of DEMs in SP cells were visualized using REVIGO, which allows to remove functionally redundant GO terms. Individual circle indicates cluster representatives. The color of circle indicates the EASE score. More functionally similar GO terms were closer in the scatterplot, but the semantic space units have no intrinsic meaning. GO terms particularly relevant to more aggressiveness of SP cells were labeled. The full list of GO terms is presented in Supplementary Table 2.
The predicted targets of DEMs were subjected to the KEGG pathway annotation to elucidate signaling networks involved in the maintenance of subpopulations of the HMM cells. The functional analysis using KEGG revealed 37 signal transduction pathways significantly involved in the SP cells compared to NSP cells (Supplementary Table 3). One of the most over-represented pathways was ErbB signaling pathway (Figure 3), and some of the other key pathways involved in the tumorigenesis of HMM included MAPK, Wnt, insulin, mTOR, and VEGF signaling pathways.
Figure 3: Target prediction of miRNAs and KEGG pathway analysis. The ErbB signaling pathway was the highest over-represented pathway among key pathways. Colors represent miRNA target genes and their expression status in KEGG pathway map. Red square: up-regulated target genes. Green square: down-regulated target genes.
Identification of differentially expressed genes (DEGs) by RNA-seq
The NGS data of the transcriptome in SP and NSP cells in the MS1 cell line was previously published [15]. By comparing the RNA-seq data of SP and NSP, as a result, the differential expression of 1,130 genes from a total of 17,122 mRNAs was identified. Among these genes, 795 DEGs were significantly up-regulated and 335 DEGs were down-regulated in the SP cells compared to the NSP cells. The NGS data of the transcriptome in the present study are available in the National Center for Biotechnology Information Sequence Read Archive (http://www.ncbi.nlm.nih.gov/Traces/sra/) under the accession number SRR2062223 and SRR2064655.
Predicted target genes of DEMs matched to mRNAs from RNA-seq analysis
To validate the microarray analysis of miRNAs, RNA-seq data were integrated to the predicted target genes of the DEMs. The RNA-seq data were used for the validation of miRNA’s target genes in SP cells in the present study. Of 1,743 predicted target genes of DEMs generated from the TargetScan database, a total of 88 gene pairs of DEM target genes and DEGs were identified. Based on the criterion of 2-fold difference, 8 up-regulated mRNAs were targeted by 3 down-regulated DEMs, and 4 down-regulated mRNAs were targeted by 3 up-regulated DEMs (Table 2). Of note, down-regulated miR-22 and up-regulated miR-138 were found to have significant association of simultaneous inverse expression in their target genes, DDIT4 (DNA-damage-inducible transcript 4) and ROCK2 (Rho-associated, coiled-coil containing protein kinase 2), respectively (p-value < 0.05). In GO and KEGG analyses, the ROCK2 gene was found to implicate cell migration through regulation of actin cytoskeleton organization and Wnt signaling pathway. In the same vein, the DDIT4 gene was shown to involve cell survival and death through the regulation of programmed cell death and mTOR signaling pathway (Table 3).
Table 2: Potential target genes of miRNAs defining SP cells identified by RNA-seq
microRNAs | Log2 FC (microarray) | miRNA-target pairs | Log2 FC (RNA-seq)a | p-valueb |
|---|---|---|---|---|
hsa-mir-7-1 | -1.831 | hsa-mir-7-1:FNDC4 | 1.693 | 0.151 |
hsa-mir-7-1:PPIF | 1.586 | 0.166 | ||
hsa-mir-7-1:PTAR1 | 1.586 | 0.166 | ||
hsa-mir-7-1:POLE4 | 1.395 | 0.166 | ||
hsa-mir-22 | -1.904 | hsa-mir-22:DDIT4 | 2.296 | 0.047 |
hsa-mir-22:FRAT2 | 1.965 | 0.116 | ||
hsa-mir-22:IRF5 | 1.639 | 0.192 | ||
hsa-mir-183 | -1.151 | hsa-mir-183:IDH2 | 1.258 | 0.269 |
hsa-mir-183* | -2.835 | hsa-mir-183*:IDH2 | 1.258 | 0.269 |
hsa-mir-138-1 | 2.420 | hsa-mir-138-1:EID1 | -1.076 | 0.328 |
hsa-mir-138-1:ROCK2 | -5.283 | 0.001 | ||
hsa-mir-663 | 1.023 | hsa-mir-663:HSPG2 | -1.728 | 0.182 |
hsa-mir-195 | 1.287 | hsa-mir-195:SLC2A3 | -1.994 | 0.146 |
Target gene prediction using RNA-seq identified a total of 12 potential target genes, showing epigenetically concomitant changes in the miRNAs and gene expression in SP cells. The expression statuses of the NSP samples were utilized as references for those in the SP samples.
a, b Fold change and p-value were determined by RNA-seq analysis
Abbreviation: FC, fold change
Table 3: Putative biological functions and pathways of ROCK2 and DDIT4 targeted by miR-122 and miR-138-1 via GO and KEGG analyses
miRNAs | Targets | GO and KEGG terms | EASE score |
|---|---|---|---|
hsa-mir-138-1 | ROCK2 | GO:0006793~phosphorus metabolic process | 0.000743227 |
GO:0006468~protein amino acid phosphorylation | 0.001408809 | ||
GO:0016310~phosphorylation | 0.003583528 | ||
GO:0051130~positive regulation of cellular component organization | 0.002402911 | ||
GO:0033043~regulation of organelle organization | 0.039881837 | ||
GO:0007010~cytoskeleton organization | 0.013056657 | ||
GO:0030036~actin cytoskeleton organization | 0.013837722 | ||
GO:0030029~actin filament-based process | 0.016937947 | ||
GO:0007242~intracellular signaling cascade | 0.004744225 | ||
KEGG:has04810~Regulation of actin cytoskeleton | 0.000105943 | ||
KEGG:has04310~Wnt signaling pathway | 0.000224177 | ||
hsa-mir-22 | DDIT4 | GO:0012501~programmed cell death | 0.009907439 |
GO:0008219~cell death | 0.003096555 | ||
GO:0006915~apoptosis | 0.007165309 | ||
GO:0009968~negative regulation of signal transduction | 2.94131E-05 | ||
GO:0010648~negative regulation of cell communication | 7.08217E-05 | ||
KEGG:hsa04150~mTOR signaling pathway | 0.017825651 |
Potential functions of miR-122 and miR-138-1 defining SP cells were presented. These miRNAs were shown to be involved in cell migration, survival, and death, via designated target gene regulation. RNA-seq identified statistical significance of ROCK2 and DDIT4 expressions in SP cells. miRTasebase revealed that the genes were experimentally verified targets of miR-122 and miR-138-1, respectively.
DISCUSSION
In the HMM, accruing evidence has emphasized the importance of aberrantly expressed miRNAs in the etiology and pathogenesis of HMM [16–18]. In the present study, microarray analysis for global expression profile of miRNAs revealed DEMs defining SP cells in a HMM cell line. A total of 95 DEMs including 42 up-regulated DEMs and 53 down-regulated DEMs were identified in SP compared to the NSP of a HMM cell line. In many human malignancies including HMM, SP is enriched for the cells with stem-like properties [14]. To the best of our knowledge, this is the first study to report the global expression profile of miRNAs in the subpopulations of HMM cells.
The notion that miRNAs are responsible for the aggressive behavior of HMM could be supported by the DEMs and their enriched GO terms found in the present study. Among the down-regulated DEMs, miR-148b is a signature of SP cells from hepatocellular carcinoma and has been implicated in enhanced metastatic and angiogenic potential of the hepatic cancer stem cells (CSCs) [19]. Similarly, miR-93 inhibits programmed cell death, facilitates cell proliferation, and promotes the colony formation of colonic CSCs [20]. Consistent with the present study, the oncogenic activation of miR-34a in CSCs with respect to cell proliferation, metastasis, drug resistance, and in vivo tumorigenicity have been reported [21–23]. Menges et al., [24] reported that aberrant expression of miR-34a was essential for the maintenance of CSC population and metastatic potential of malignant mesothelial cells. They showed that CSC populations lacking miR-34a expression were highly heterogenic and invasive in the genetically modified mice harboring NF2 and CDKN2A gene suppression [25]. One of the upregulated DEMs in our study, miR-138, has been reported to enhance cell survival, in vitro tumor sphere formation and in vivo tumorigenicity glioma [26].
The SP cells of HMM cell lines are more resistant to the chemotherapeutic drugs than the NSP cells [27]. However, the underlying mechanism of the drug resistance is not comprehensively understood. A few miRNAs identified in this study may be attributable to the enhanced drug resistance of SP cells. Up-regulation of miR-125b-2 increases the drug resistance of glioblastoma CSCs through silencing the Bcl-2 family and inhibition of mitochondria-dependent apoptosis [28]. Likewise, hypomethylation and concordant overexpression of miR-663 reduces the drug sensitivity of human breast cancer cells by repressing the expression of heparan sulfate proteoglycan 2 (HSPG2) [29]. The present study also showed that the expression of HSPG2 was suppressed by up-regulation of miR-663 in HMM cells.
Cancer cells cope with stress by stimulating the expression levels of miRNAs and their target genes [30–32]. The GO analysis in the present study showed that approximately 40% of miRNAs defining SP cells were categorized into biological processes of cellular response to stress. Among the DEMs, miR-7-1 and let-7d may be critical with regard to the stress regulation of SP cells in HMM. Both miRNAs are commonly dysregulated in HMM, promoting malignancy through the activation of EGF, PDGFA, and RAS oncogenes [17, 18, 33]. Flies with mutated miR-7 fail to develop their eyes when grown upon temperature fluctuations [34]. The let-7 family that is commonly down-regulated in many tumor types triggers stress phenotypes of tumor cell [35]. Although further studies are warranted for the detailed mechanisms, SP cells may be quite versatile in adapting to hostile microenvironments.
The ErbB signaling was the one of the most over-represented pathways identified by the KEGG pathway mapping of DEMs and their target genes in this study. The data indicate that the ErbB signaling pathway mediated via dysregulated miRNAs may be the key oncogenic activation of HMM. ErbB tyrosine kinase receptors have been associated with cell proliferation, survival, and differentiation in human solid tumors [36]. The activation of EGFR (ErbB1) and ErbB2 increases the long term survival of asbestos-exposed mesothelial cells in vivo and dissemination of HMM [37]. Among the multiple signaling pathways recognized in the present study, the ErbB signaling may be significantly involved in self-renewal and/or survival of more aggressive cancer cell subpopulation in HMM.
Integrated analysis of miRTarBase and RNA-seq revealed a number of potential target genes. Among them, ROCK2 and DDIT4 were found to be statistically significant targets of miR-138 and miR-22, respectively. The ROCK2 is involved in cell migration and phosphorus metabolism in SP cells, and reduced expression of ROCK2 inhibits migratory and invasive properties of cancer cells [38]. The DDIT4 is transcribed when exposed to cellular stress such as hypoxia and DNA damage, and it enhances cell growth and survival by inhibiting the cascades of the mammalian target of rapamycin complex 1 [39]. Although further studies for functional validation of target genes of DEMs are necessary, the present study indicates that altered expression of miR-138 and miR-22 may be associated with the maintenance of tumor heterogeneity HMM by regulating their target gene expression.
The failure of the first-line therapy and limited availability of second-line treatment options call for an urgent need of novel strategies to improve the prognosis of HMM patients [10]. Studying heterogeneous cancer cell population within HMM cell lines and tissue may elucidate the molecular mechanisms underlying the generation of tumor cells refractory to current therapies [12]. The DEMs identified in the present study could provide valuable information about the regulatory function of miRNAs in the development of intratumoral heterogeneity of HMM, opening an avenue to devise a novel strategy to cope with the malignant disease. The present study, however, is limited by the lack of a validated experiment. Thus, future research is warranted to scrutinize the mechanism by which miRNAs contribute to the tumor aggressiveness and heterogeneity in HMM.
MATERIALS AND METHODS
HMM cell line and culture
A HMM cell line, MS1, was kindly provided by Dr. Jablons (University of California San Francisco). The MS1 cells were determined to be free of mycoplasma contamination by using e-Myco Mycoplasma PCR detection kit (e-Myco, iNtRON Biotechnology, Sungnam, Korea). The cells were cultured in conventional RPMI 1640 medium (Mediatech Inc., Manassas, VA, USA) with 10% fetal bovine serum (FBS; Mediatech Inc.) and supplements at 37°C in a humidified atmosphere containing 5% CO2, as previously described [13].
SP assay and cell sorting
A detailed protocol of SP assay composed of Hoechst 33342 dye staining and subsequent flow cytometry analysis was previously described [13]. Briefly, 106 cells/mL in pre-warmed RPMI containing 2% FBS and 10mM HEPES were incubated with Hoechst 33342 dye (5μg/mL; Sigma–Aldrich, St. Louis, MO, USA) for 90 min at 37°C with intermittent mixing. After centrifugation at 480 × g for 5 minute and subsequent PBS washing containing 2 % FBS at 4°C, the incubated cells were subjected to SP analysis. The flow cytometer sorter (Becton-Dickinson FACS Aria III, Becton-Dickinson, San Jose, CA, USA) equipped with Hoechst Blue with a 450/50 broad pass (BP) filter and Hoechst Red with a 675/30 BP filter was used to detect Hoechst 33342 staining. Verapamil hydrochloride (50 μM, Sigma–Aldrich, St. Louis, MO, USA) which blocks Hoechst 33342 dye efflux was used to identify SP fractions. Three biological replicates of SP and NSP cells were pooled and subjected to further analyses.
RNA isolation
Total RNAs, including miRNAs, were isolated from the sorted SP and NSP cells using miRNeasy extraction kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s instructions. The total RNAs extracted from the sorted cells were pooled and subjected to the miRNA array analysis. RNA quality was assessed by Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) using the RNA 6000 Nano Chip, and quantity was determined using Nanodrop-1000 Spectrophotometer (Thermo scientific, Wilmington, DE, USA).
Expression profiling of human microRNAs by microarray chip assay
Total RNAs with high quality isolated from SP and NSP cells were subjected to microarray assay performed at DNA Link Inc., (Songpa-gu, Seoul, The Republic of Korea). Per RNA sample, 1.6 μg was used as an input into the Affymetrix procedure as recommended by manufacturer’s protocol (Affymetrix® FlashTag™ Biotin HSR RNA Labeling Kits, cat. no. HSR30FTA; Genisphere, LLC, Hatfield, PA, USA). Briefly, 1.6 μg of total RNA was tailed by poly (A) and end-labeled by FlashTag™ sequence that has biotin using peroxidase-antiperoxidase enzyme and T4 DNA Ligase, respectively. End-labeled miRNA was hybridized to the GeneChip® miRNA 3.0 arrays (Affymetrix Inc., Santa Clara, CA, USA) for 16 hours at 48°C and 60 rpm. After the hybridization, the chips were stained and washed in a GeneChip Fluidics Station 450 (Affymetrix Inc.) and scanned using a GeneChip Array scanner 3000 7G (Affymetrix Inc.). The expression intensity data were extracted from the scanned images using Affymetrix Command Console software version 1.1 and stored as CEL files.
Data analysis and miRNA target prediction
The intensity values of CEL files were normalized to remove bias between the arrays, using the Robust Multi-array Average (RMA) and Detected Above Background (DABG) algorithm implemented in the Affymetrix Expression Console software (version 1.3.1.) (http://www.affymetrix.com) [40]. The whole normalized data were imported into the programming environment R (version 3.0.2) and overall signal distributions of each array were compared by plotting using tools available from the Bioconductor Project (http://www.bioconductor.org) to check good normalization [41]. After confirming the normalization of the data, DEMs that showed over 2-fold difference between the average signal values of NSP cells and SP cells were selected. Additionally, the normalized data of the selected miRNAs were also imported into the programming environment R for the statistical t-test and genes with p-value less than 0.05 were extracted as candidate DEMs for further investigation [41]. For the candidate DEMs, targets that have more than 95 context score were subjected to computational prediction algorithm TargetScan 6.2 database (http://www.targetscan.org).
Gene Ontology and KEGG pathways
Predicted targets of the candidate DEMs identified by the TargetScan 6.2 database were functionally annotated and classified based on the gene functions listed in gene ontology (GO) databases including a web-based tool DAVID (the Database for Annotation, Visualization, and Integrated Discovery) (http://david.abcc.ncifcrf.gov) [42]. Fisher’s exact test was used to determine statistical significance between the DEMs and those on the GO annotation list. The GO terms scoring EASE value of < 0.05 were considered statistically significant [43]. The redundancy in the resulting set of GO terms were removed by using REVIGO web server (http://revigo.irb.hr) [44]. SimRel as semantic similarity measure was used, and the allowed similarity was medium (0.7) between the enriched GO terms [44].
Subsequently the predicted targets of the candidate DEMs were subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, a pathway mapping tool for molecular networks to identify the biological pathways. The KEGG pathways scoring EASE value of < 0.05 were considered statistically significant.
RNA-seq and identification of DEGs
The construction of RNA library and detailed analysis of the RNA-seq data was previously performed as described [15]. Briefly, total RNAs with high quality from the sorted SP and NSP cells were subjected to Next Generation Sequencing (NGS) assay performed at the DNA Link Incorporation (Songpa-gu, Seoul, The Republic of Korea). Sequencing libraries of mRNAs were prepared using an Illumina TruSeq RNA Prep kit v2 (Illumina Inc, San Diego, CA, USA.) according to the manufacturer’s instructions. The quality of the amplified libraries was verified using an Agilent Technologies 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). Cluster generation was carried out in the flow cells on the cBot automated cluster generation system (Illumina Inc.), and then the flow cells were loaded on a HiSeq 2000 sequencing system (Illumina Inc.) with 200 bps paired-end reads. The high quality of clean reads was mapped to the hg19 with TopHat (ver. 2.0.9). The Bam file was used as the output to store a list of read alignments and was added to the Cufflinks software package (ver. 2.0.2) to predict transcript structures and compare transcriptome profiles based on the RNA-Seq data [45]. To compare the expression level of a gene across samples, read counts obtained from RNA-seq were normalized as fragments per kilobase of transcript per million mapped fragments (FPKM) [46]. The FPKM was used to identify DEGs in SP and NSP subpopulations, and then the FPKM in each sample was compared and transformed to the Log2 ratio (log2(number of SP reads) – log2(number of NSP reads)). The gene expression of the NSP subpopulation was used as control data for the determination of up- or down-regulated genes in SP cells. Genes with a p-value of < 0.05 and a log2-transformed value smaller than -1 or greater than 1 were considered to be statistically significant DEGs.
Integrated analysis of DEM targets and DEGs
To validate potential targets of the candidate DEMs identified in the study, a combined prediction of miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/) and RNA-seq was carried out. The miRTarBase is a reliable online software program to search experimentally validated miRNA-target interactions (MTIs) [47]. The RNA-seq that is based on deep-sequencing technology provides a transcription profile in cells, tissues, and organisms [48]. The miRNAs-target mRNA pairs predicted from the TargetScan 6.2 database were analyzed with the miRTarBase data in an integrated way to identify the common miRNA-target pairs. Subsequently, the inverse expression of the potential MTIs was investigated using RNA-seq data. Targets showing over 2-fold difference and p-value less than 0.05 between SP and NSP cells were selected.
Author’s contribution
Contributed to conception or design: MCK and YK.
Contributed to acquisition, analysis, or interpretation of data: MCK, YK, and NK.
Performed the experiments: MCK and US.
Drafted the manuscript: MCK and YK
ACKNOWLEDGMENTS
This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (Grant #: NRF-2013R1A2A2A01068237) and the BK21 program and the Research Institute of Veterinary Sciences, College of Veterinary Medicine, Seoul National University. The authors are grateful to the DNA Link Incorporation (Songpa-gu, Seoul, The Republic of Korea) for their excellent technical assistance.
COMPETING INTERESTS
The authors have declared that no competing interests exist.
REFERENCES
1. Marusyk A, Polyak K. Tumor heterogeneity: causes and consequences. Biochimica et Biophysica Acta (BBA)-Reviews on Cancer. 2010; 1805: 105-17.
2. Pietras A. Cancer stem cells in tumor heterogeneity. Adv Cancer Res. 2011; 112: 255-81.
3. Saunders NA, Simpson F, Thompson EW, Hill MM, Endo-Munoz L, Leggatt G, Minchin RF, Guminski A. Role of intratumoural heterogeneity in cancer drug resistance: molecular and clinical perspectives. EMBO Mol Med. 2012; 4: 675-84.
4. Wilting RH, Dannenberg JH. Epigenetic mechanisms in tumorigenesis, tumor cell heterogeneity and drug resistance. Drug Resistance Updates. 2012; 15: 21-38.
5. Heng HH, Bremer SW, Stevens JB, Ye KJ, Liu G, Ye CJ. Genetic and epigenetic heterogeneity in cancer: A genome-centric perspective. J Cell Physiol. 2009; 220: 538-47.
6. Feinberg AP, Ohlsson R, Henikoff S. The epigenetic progenitor origin of human cancer. Nat Rev Genet. 2006; 7: 21-33.
7. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004; 116: 281-97.
8. Volinia S, Calin GA, Liu CG, Ambs S, Cimmino A, Petrocca F, Visone R, Iorio M, Roldo C, Ferracin M. A microRNA expression signature of human solid tumors defines cancer gene targets. Proc Natl Acad Sci U S A. 2006; 103: 2257-61.
9. Robinson BW, Lake RA. Advances in malignant mesothelioma. New Engl J Med. 2005; 353: 1591-603.
10. Cao S, Jin S, Cao J, Shen J, Hu J, Che D, Pan B, Zhang J, He X, Ding D. Advances in malignant peritoneal mesothelioma. Int J Colorectal Dis. 2015; 30: 1-10.
11. Sun X, Wei L, Lidén J, Hui G, Dahlman-Wright K, Hjerpe A, Dobra K. Molecular characterization of tumour heterogeneity and malignant mesothelioma cell differentiation by gene profiling. The Journal of pathology. 2005; 207: 91-101.
12. Carbone M, Gaudino G, Yang H. Recent Insights Emerging from Malignant Mesothelioma Genome Sequencing. J Thorac Oncol. 2015; 10: 409-11.
13. Kai K, D'Costa S, Yoon BI, Brody AR, Sills RC, Kim Y. Characterization of side population cells in human malignant mesothelioma cell lines. Lung Cancer. 2010; 70: 146-51.
14. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT. Gene Ontology: tool for the unification of biology. Nature genetics. 2000; 25: 25-9.
15. Kim MC, Kim NY, Seo YR, Kim Y. An Integrated Analysis of the Genome-Wide Profiles of DNA Methylation and mRNA Expression Defining the Side Population of a Human Malignant Mesothelioma Cell Line. Journal of Cancer. 2016; 7: 1668.
16. Quinn L, Finn SP, Cuffe S, Gray SG. Non-coding RNA repertoires in malignant pleural mesothelioma. Lung Cancer. 2015; 90: 417-26.
17. Balatti V, Maniero S, Ferracin M, Veronese A, Negrini M, Ferrocci G, Martini F, Tognon MG. MicroRNAs dysregulation in human malignant pleural mesothelioma. J Thorac Oncol. 2011; 6: 844-51.
18. Truini A, Coco S, Alama A, Genova C, Sini C, Dal Bello M, Barletta G, Rijavec E, Burrafato G, Boccardo F. Role of microRNAs in malignant mesothelioma. Cell Mol Life Sci. 2014; 71: 2865-78.
19. Liu Q, Xu Y, Wei S, Gao W, Chen L, Zhou T, Wang Z, Ying M, Zheng Q. microRNA-148b suppresses hepatic cancer stem cell by targeting neuropilin-1. Biosci Rep. 2015; 35: e00229.
20. Yu XF, Zou J, Bao ZJ, Dong J. miR-93 suppresses proliferation and colony formation of human colon cancer stem cells. World journal of gastroenterology. 2011; 17: 4711.
21. Ji Q, Hao X, Meng Y, Zhang M, DeSano J, Fan D, Xu L. Restoration of tumor suppressor miR-34 inhibits human p53-mutant gastric cancer tumorspheres. BMC Cancer. 2008; 8: 266.
22. Guessous F, Zhang Y, Kofman A, Catania A, Li Y, Schiff D, Purow B, Abounader R. microRNA-34a is tumor suppressive in brain tumors and glioma stem cells. Cell cycle. 2010; 9: 1031-6.
23. Bu P, Chen KY, Chen JH, Wang L, Walters J, Shin YJ, Goerger JP, Sun J, Witherspoon M, Rakhilin N. A microRNA miR-34a-regulated bimodal switch targets notch in colon cancer stem cells. Cell stem cell. 2013; 12: 602-15.
24. Menges CW, Kadariya Y, Altomare D, Talarchek J, Neumann-Domer E, Wu Y, Xiao GH, Shapiro IM, Kolev VN, Pachter JA. Tumor suppressor alterations cooperate to drive aggressive mesotheliomas with enriched cancer stem cells via a p53-miR-34a-c-Met axis. Cancer Res. 2014; 74: 1261-71.
25. Tanaka N, Toyooka S, Soh J, Tsukuda K, Shien K, Furukawa M, Muraoka T, Maki Y, Ueno T, Yamamoto H. Downregulation of microRNA-34 induces cell proliferation and invasion of human mesothelial cells. Oncol Rep. 2013; 29: 2169-74.
26. Chan XHD, Nama S, Gopal F, Rizk P, Ramasamy S, Sundaram G, Ow GS, Vladimirovna IA, Tanavde V, Haybaeck J. Targeting glioma stem cells by functional inhibition of a prosurvival oncomiR-138 in malignant gliomas. Cell reports. 2012; 2: 591-602.
27. Cortes-Dericks L, Carboni GL, Schmid RA, Karoubi G. Putative cancer stem cells in malignant pleural mesothelioma show resistance to cisplatin and pemetrexed. Int J Oncol. 2010; 37: 437-44.
28. Shi L, Zhang S, Feng K, Wu F, Wan Y, Wang Z, Zhang J, Wang Y, Yan W, Fu Z. MicroRNA-125b-2 confers human glioblastoma stem cells resistance to temozolomide through the mitochondrial pathway of apoptosis. Int J Oncol. 2012; 40: 119-29.
29. Hu H, Li S, Cui X, Lv X, Jiao Y, Yu F, Yao H, Song E, Chen Y, Wang M. The overexpression of hypomethylated miR-663 induces chemotherapy resistance in human breast cancer cells by targeting heparin sulfate proteoglycan 2 (HSPG2). J Biol Chem. 2013; 288: 10973-85.
30. Leung AK, Sharp PA. MicroRNA functions in stress responses. Mol Cell. 2010; 40: 205-15.
31. Marsit CJ, Eddy K, Kelsey KT. MicroRNA responses to cellular stress. Cancer Res. 2006; 66: 10843-8.
32. Mendell JT, Olson EN. MicroRNAs in stress signaling and human disease. Cell. 2012; 148: 1172-87.
33. Guled M, Lahti L, Lindholm PM, Salmenkivi K, Bagwan I, Nicholson AG, Knuutila S. CDKN2A, NF2, and JUN are dysregulated among other genes by miRNAs in malignant mesothelioma-a miRNA microarray analysis. Genes Chromosomes Cancer. 2009; 48: 615-23.
34. Li X, Cassidy JJ, Reinke CA, Fischboeck S, Carthew RW. A microRNA imparts robustness against environmental fluctuation during development. Cell. 2009; 137: 273-82.
35. Iliopoulos D, Hirsch HA, Struhl K. An epigenetic switch involving NF-κB, Lin28, Let-7 MicroRNA, and IL6 links inflammation to cell transformation. Cell. 2009; 139: 693-706.
36. Grant S, Qiao L, Dent P. Roles of ERBB family receptor tyrosine kinases, and downstream signaling pathways, in the control of cell growth and survival. Front Biosci. 2002; 7: 389.
37. Thirkettle I, Harvey P, Hasleton P, Ball R, Warn R. Immunoreactivity for cadherins, HGF/SF, met, and erbB-2 in pleural malignant mesotheliomas. Histopathology. 2000; 36: 522-8.
38. Jiang L, Liu X, Kolokythas A, Yu J, Wang A, Heidbreder CE, Shi F, Zhou X. Downregulation of the Rho GTPase signaling pathway is involved in the microRNA-138-mediated inhibition of cell migration and invasion in tongue squamous cell carcinoma. Int J Cancer. 2010; 127: 505-12.
39. Sofer A, Lei K, Johannessen CM, Ellisen LW. Regulation of mTOR and cell growth in response to energy stress by REDD1. Mol Cell Biol. 2005; 25: 5834-45.
40. Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003; 4: 249-64.
41. Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, Hornik K, Hothorn T, Huber W, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004; 5: R80. doi: 10.1186/gb-2004-5-10-r80.
42. Sherman BT, Huang DW, Tan Q, Guo Y, Bour S, Liu D, Stephens R, Baseler MW, Lane HC, Lempicki RA. DAVID Knowledgebase: a gene-centered database integrating heterogeneous gene annotation resources to facilitate high-throughput gene functional analysis. BMC Bioinformatics. 2007; 8: 426.
43. Hosack DA, Dennis G Jr, Sherman BT, Lane HC, Lempicki RA. Identifying biological themes within lists of genes with EASE. Genome Biol. 2003; 4: R70.
44. Supek F, Bošnjak M, Škunca N, Šmuc T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS One. 2011; 6: e21800.
45. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nature protocols. 2012; 7: 562-78.
46. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature methods. 2008; 5: 621-8.
47. Hsu SD, Lin FM, Wu WY, Liang C, Huang WC, Chan WL, Tsai WT, Chen GZ, Lee CJ, Chiu CM, Chien CH, Wu MC, Huang CY, et al. miRTarBase: a database curates experimentally validated microRNA-target interactions. Nucleic Acids Res. 2011; 39:D163-9. doi: 10.1093/nar/gkq1107.
48. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009; 10: 57-63.
