Integrated analysis of microRNA regulatory network in nasopharyngeal carcinoma with deep sequencing

Background MicroRNAs (miRNAs) have been shown to play a critical role in the development and progression of nasopharyngeal carcinoma (NPC). Although accumulating studies have been performed on the molecular mechanisms of NPC, the miRNA regulatory networks in cancer progression remain largely unknown. Laser capture microdissection (LCM) and deep sequencing are powerful tools that can help us to detect the integrated view of miRNA-target network. Methods Illumina Hiseq2000 deep sequencing was used to screen differentially expressed miRNAs in laser-microdessected biopsies between 12 NPC and 8 chronic nasopharyngitis patients. The result was validated by real-time PCR on 201 NPC and 25 chronic nasopharyngitis patients. The potential candidate target genes of the miRNAs were predicted using published target prediction softwares (RNAhybrid, TargetScan, Miranda, PITA), and the overlay part was analyzed in Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) biological process. The miRNA regulatory network analysis was performed using the Ingenuity Pathway Analysis (IPA) software. Results Eight differentially expressed miRNAs were identified between NPC and chronic nasopharyngitis patients by deep sequencing. Further qRT-PCR assays confirmed 3 down-regulated miRNAs (miR-34c-5p, miR-375 and miR-449c-5p), 4 up-regulated miRNAs (miR-205-5p, miR-92a-3p, miR-193b-3p and miR-27a-5p). Additionally, the low level of miR-34c-5p (miR-34c) was significantly correlated with advanced TNM stage. GO and KEGG enrichment analyses showed that 914 target genes were involved in cell cycle, cytokine secretion and tumor immunology, and so on. IPA revealed that cancer was the top disease associated with those dysregulated miRNAs, and the genes regulated by miR-34c were in the center of miRNA-mRNA regulatory network, including TP53, CCND1, CDK6, MET and BCL2, and the PI3K/AKT/ mTOR signaling was regarded as a significant function pathway in this network. Conclusion Our study presents the current knowledge of miRNA regulatory network in NPC with combination of bioinformatics analysis and literature research. The hypothesis of miR-34c regulatory pathway may be beneficial in guiding further studies on the molecular mechanism of NPC tumorigenesis. Electronic supplementary material The online version of this article (doi:10.1186/s13046-016-0292-4) contains supplementary material, which is available to authorized users.


Background
Nasopharyngeal carcinoma (NPC) is a highly invasive and metastatic cancer that is widely prevalent in southern China [1,2]. During tumorigenesis and progression, multiple genetic and epigenetic abnormalities synergistically disrupt normal cell function, especially playing an important role in NPC pathogenesis [3,4]. In recent decades, microRNAs (miRNAs) have been increasingly recognized as important genetic regulators in the mammalian system [5]. Dysregulation of miRNA expression in NPC can regulate tumor cell growth, differentiation and apoptosis.
MiRNAs are endogenous, small (18-25 nt), non-proteincoding RNA molecules which modulate many physiological and pathological processes through down-regulating target genes [6]. Previous studies mostly focused on single miRNA, such as miR-138, a potential tumor suppressor by targeting CCND1 [7]. But microRNA expression profiles have been demonstrated to be unique for a better understanding of different stages of tumor progression and metastasis. Hence, it is necessary to use high-tech tools to get a comprehensive miRNA analysis of NPC etiology.
Thus far, several high-throughput techniques have been applied to identify deregulated miRNAs in NPC. Most studies to date have applied microarray technique on nasopharyngeal carcinoma and non-cancer nasopharyngitis tissues, and some of them could identify the miRNAs as a prognostic factor or recurrent marker from the distinctive miRNA expression profile [8][9][10]. But the results failed to show good interplatform concordance. Both Sengupta et al. and Luo et al. recruited microarray on laser-microdissected tissues which including NPC and normal surrounding epithelial cells, only three miRNAs (miR-34c-5p, miR-29c and miR-34b-5p) were shared in the list of abnormal expressions by comparing Sengupta's data and Luo's data [11,12]. The inconsistency may be explained by heterogeneity of samples and use of different expression measuring platforms. To get the better data reproducibility, we need some more robust tools to reduce the disturbance of heterogeneity.
The advent of the deep sequencing provides a rapid and high throughput tool that can be used to explore the large miRNA pool, and possesses obvious advantages for the identification of miRNA sequence variations and the discovery of novel miRNAs [13]. Plieskatt et al. compared miRNA expression in FFPE ( formalin fixed paraffin-embedded tissue) from NPC cases and controls using both microarray and RNA-Seq technologies, showed that RNA-Seq could additionally indicate unknown miRNAs associated with NPC [14]. Deep sequencing technology was considered as a more powerful and accurate tool in expression profile study, which would allow us to generate a comprehensive insight into the networks constructed by many more deregulated miRNAs not described in previous studies, just like the study of Chang et al., which delivered a clear picture of the global miRNA regulatory characteristics in triple-negative breast cancer by deep sequencing [15].
To better characterise the specific signature of NPC cells, we expand our observations by applying laser capture microdissection (LCM) on NPC and chronic nasopharyngitis samples to compare the miRNA signatures between them. LCM can separate the truly transformed cancerous cells from those other cell types commonly present in a tumor tissue, such as immune cells, connective tissue and new vasculature [16][17][18]. Wang et al. demonstrated the feasibility and potential power of discovering miRNA biomarkers in colorectal tumor tissue using the combination of LCM with miRNA microarrays [19], and with the combination of LCM and microarray, Sengupta et al. identified miR-29c was important in regulating genes involved in metastasis of NPC [11].
Herein, we aim to examine the miRNA expression profile in clinical tissue samples of NPC patients, and get better understanding of the miRNA regulatory network. After bioinformatics analysis and literature review, we identify miR-34c as an important partner in NPC tumorigenesis, and propose that miR-34c-target genes network (including CCND1, TP53, BCL2, CDK6, MET) should be critical in mediating NPC cellular proliferation, migration and invasion.

Patient samples and laser capture microdissection
A total of 213 NPC patients and 33 chronic nasopharyngitis patients were recruited prospectively from Nanfang Hospital (Southern Medical University, Guangzhou, China). None of the patients had received anticancer treatments before undergoing biopsy. Fresh tissues were snap frozen in liquid nitrogen and stored at −160°C. All samples were pathologically confirmed. The patients were informed about the sample collection and had signed informed consent forms. TNM classification was according to the definitions of the seventh edition of the UICC-American Joint Committee on Cancer staging criteria. 12 NPC biopsy specimens and 8 non-cancer nasopharyngitis biopsy samples were chose for laser capture microdissection (LCM). The controls were carefully selected to match the gender and age distribution of NPC patients. Because the sample of stage I failed to pass the quality control, we chose 2 cases of stage II, and 5 cases of each stage III IV (Table 1). Specimens were first frozen-sectioned by using a LEICA CM 1900 cryomicrotome. After hematoxylin and eosin (H&E) staining, LCM was performed on a MMI Cellcut Microdissection Instrument (Molecular Machines & Industries, Swiss). Phase contrast images were acquired using OlympusIX71 microscope. The research protocols were approved by the Ethics Committee of Nanfang Hospital and registered in Clinical.trials.gov (ID: NCT01171235).

Small RNA deep sequencing
Total RNA, containing miRNA, was extracted using Trizol Reagent (Invitrogen, CA) from laser-microdessected samples, and passed the RNA quality control for sequencing. The quality and integrity of total RNA was assessed with Agilent 2100 Bioanalyzer (Agilent Technologies, USA). To achieve optimal tissue miRNA profiles, we carried out high-throughput next-generation sequencing (Illumina, BGI, Shenzhen) of 12 NPC samples and 8 chronic nasopharyngitis samples by following the manufacturer's recommended protocols. We screened the high quality clean read sequences by the alignment to NCBI GenBank data and miRBase 21.0 for the further analysis. The raw data have been submitted to NCBI under BioProject accession No.PRJNA 289899.

Differential miRNA expression analysis
To identify miRNAs differentially expressed between NPC patients and controls, we applied transcripts per million (TPM) to normalize the expression of miRNA in two groups (NPC and controls). Then we calculated fold change (FC) and P-value via T-test, and corrected P-value into false discovery rate (FDR) using the Benjamin and Hochberg method [20]. FDR ≤ 0.05 and | log2FC | ≥ 2 were set as the cut-offs to screen out differentially expressed miRNAs.

Quantitative reverse transcription-PCR analysis of miRNA expression levels
Quantitative reverse transcription-PCR was used to validate the sequencing results. Total RNA from 201 NPC and 25 chronic nasopharyngitis biopsy specimens was extracted using Trizol Reagent (Invitrogen, CA). Reverse transcription of the total RNA was performed using Allin-One First-Strand cDNA Synthesis kit (GeneCopoeia Inc., USA) according to the manufacturer's protocol. Real-time PCR was performed by using All-in-OneTM qPCR Mix (Applied GeneCopoeia Inc., USA) on Roche Lightcycler 480 System. U6 snRNA was used as miRNA endogenous control respectively. All samples were normalized to internal control and fold changes were calculated through relative quantification [21].

Bioinformatics analysis
Four softwares (RNAhybrid, TargetScan, Miranda, PITA) were used for target gene prediction and only the genes identified by all four approaches were selected out [22][23][24][25]. And we chose the overlapped genes targeted by oncogenic miRNAs for further study, as well as the tumor suppressor miRNAs. To understand the functions of predicted target genes, Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis were performed with package GO stats (http://www.geneontology.org/) of P value <0.05 was set as the cut-off to select out significantly enriched terms [26]. Then the miRNA regulatory network analysis was performed using the IPA software (http://www.ingenuity.com). The created genetic networks describe functional relationships among miR-NAs and genes based on known associations in the databases [27]. Networks related were ranked according to their biological relevance to the gene list provided.

Statistical analysis
Statistical analyses were conducted using spss19.0 software. The data are shown as the mean ± SEM unless otherwise noted. Two-tailed Student's t test was used for identify differentially expressed miRNAs between NPC and controls. One-way ANOVA was used to show significant associations between the miR-34c-5p level and clinicopathological parameters. P values of <0.05 were considered statistically significant.

Validation of expression of miRNAs
To confirm the deep sequencing results, we used qRT-PCR to assess the expression of 8 miRNAs with an independent cohort (including 65 NPC patients and 20 chronic nasopharyngitis patients). The relative level of each miRNA was shown in Fig. 2a. Expression of 7 miRNAs measured by qRT-PCR was dramatically different between NPC and chronic nasopharyngitis tissues and significantly correlated with their sequencing data. Because miR-92b-3p was not reliably measured by qRT-PCR in the tissue specimens, it was excluded from further analysis. The dynamic expression levels of miRNAs were revealed that the patterns were classified into 2 groups, and the varying tendencies between control and NPC were consistent with the sequencing result (Fig. 2b, c). Moreover, the expression level of miR-34c decreased with the ascent of clinical stage, which implicated that miR-34c might be more important in NPC progression.
Prediction and function analysis of target genes of miRNAs Functional analysis was conducted on genes predicted as targets of 7 differentially expressed miRNAs by taking the overlay part of 4 software results (RNAhybrid, TargetScan, Miranda, PITA). 666 genes targeted by downregulated miR-NAs and 248 genes targeted by upregulated miRNAs were acquired after screening (Additional file 1: Table S1, S2).
The selected genes were analyzed by the Gene ontology (GO) enrichment and Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis. The enrichment analysis of GO categories included biological process (BP), molecular function (MF), and cellular component (CC) three parts. The highest enrichment were the terms like "Interleukin-1 secretion", "Mast cell activation involved in immune response", "Fibroblast growth factor binding", "Lipopolysaccharide binding" and "Cortical actin cytoskeleton" (Table 3). These results showed the function of target genes mostly focused on cell growth, secretion, migration and immune system process.
The result of KEGG pathway analysis was enumerated in Table 4. Chemokine signaling pathway, Cytokinecytokine receptor interaction and Cell cycle, Natural killer cell mediated cytotoxicity, etc. were remarkable pathways which could further confirm the target genes' function in tumor cell proliferation, secretion and tumor immunology. So the differentially expressed miRNAs might be involved in the development of NPC by targeting these genes.

Regulatory networks for differentially expressed miRNAs
The biological function analysis of target genes from Ingenuity pathway analysis (IPA) database showed that cancer was the top functional category which was significantly related to the ectopic miRNAs, along with embryonic and organismal development, which also could be associated with cancer ( Table 5). The filtering used in IPA allowed us to connect 6 miRNAs with target genes. In the regulation network of miRNA, miR-34 family  (including miR-34c-5p, miR-449c-5p) targeted the most genes, and 5 directions of miR-34c come from experimental observations were highlighted (Fig. 3a). That indicated the relationships between miR-34c and the 5 genes (TP53, CCND1, BCL2, CDK6, MET) in cancer have already been confirmed. The binding sites of miR-34c on these 5 genes could be found through TargetScan software (Additional file 2: Figure S1). CCND1, TP53 and BCL2 were clearly in the center of network, so their related networks were picked out for pathway construction. TP53 has been largely confirmed as a tumor suppressor gene that mediated cell cycle regulation with its downstream factors and could be treated as diagnostic marker of NPC (Fig. 3b). And CCND1 and TP53 were both involved in PI3K/AKT/mTOR signaling that modulated tumor cells growth and proliferation (Fig. 3b, c). Sometimes they even induced apoptosis or autophagy. Base on the IPA data, BCL2 could also be predicted to be the diagnostic marker of NPC as TP53, while it acted significantly in molecular mechanisms of cancer and functioned through apoptosis signaling (Fig. 3d).

Comparison of miR-34c levels during NPC progression
As the previous results showed miR-34c might be more valuable than others, miR-34c levels in more tissue samples using qRT-PCR were further validated. A total of 201 NPC tissue samples were used to validate miR-34c via qRT-PCR, it was significantly reduced in NPC patients as compared with controls (P < 0.01, Fig. 4a), and showed the significant correlation with clinical stages respectively (P < 0.05, Fig. 4b). All these data suggested that low concentration of miR-34c was associated with higher tumour stage.

MiR-34c inhibited MET, CCND1, CDK6, BCL2 expression in NPC cells
To further characterize whether these target genes (MET, CCND1, CDK6, BCL2) respond to miR-34c in NPC cells, CNE-2 cells were transfected with miR-34c mimic or miR-  Ctrl. Western blot was performed to detect the expression level of MET, CCND1, CDK6 and BCL2. As shown in Fig. 4c, the overexpression of miR-34c with miR-34c mimics led to a notable decrease in these four genes protein levels compared with the negative control in CNE-2 cells. These results indicated that MET, CCND1, CDK6 and BCL2 were actually targeted by miR-34c in NPC cells.

Discussion
We analyzed the miRNA expression profile from 12 NPC and 8 chronic nasopharyngitis tissues after LCM with the RNA-seq technology in this study. Deregulated miRNAs were screened from the statistical analyses and a panel of 7 significantly differentially expressed miRNAs (miR-34c-5p, miR-375, miR-449c-5p, miR-205-5p, miR-92a-3p, miR-193b-3p and miR-27a-5p) was found to be an effective regulatory factor between NPC and controls. And our result was also consistent with many studies of other tumors, including miR-34c and miR-375, which were both downregulated in colorectal carcinoma [28,29], and miR-449c was showed to inhibit gastric carcinoma growth [30]. On the contrary, miR-205, miR-92a, miR-193b and miR-27a-5p were reported to be upregulated in other kinds of cancers, including lung cancer, cervical cancer, glioma and renal cell carcinoma [31][32][33][34]. According to our expression profile of miRNA, there is no study about the deregulated expression of miR-27a-5p, miR-193-3p, miR-449c-5p or miR-92a-3p in NPC yet. Mir-375 expression was significantly reduced in NPC and inhibited tumor growth by targeting MTDH [35]. The results of Tang et al. implicated miR-205-5p may be a novel NPC candidate biomarker, which was showed to determine the radioresistance of NPC through the miR-205-PTEN-Akt pathway [36]. The reduced expression of miR-34c frequently occurred in many studies of NPC, which was believed to play important role in NPC. The recent study of Li et al. proved that miR-34c suppressed NPC cell viability, colony formation by targeting MET [37]. Our study also showed the miR-34c level in NPC tissues got a decreasing trend during NPC progression via qRT-PCR.
Nevertheless, it has been known that miRNAs always functioned not only through one single gene, but a huge network which covered multiple genes. For a deeper understanding of the functional genomic alterations induced by deregulated miRNAs, target genes were predicted by multiple softwares included RNAhybrid, TargetScan, Miranda and PITA, and the bioinformatics analysis of the functional category and integrated network were conducted, including GO, KEGG enrichment analysis and IPA. The GO enrichment analysis found that the target genes of these dysregulated miRNAs were associated with cell growth, protein binding and secretion, which were all involved in tumorigenesis, KEGG analysis also obtained similar results (Table 3, Table 4). Interestingly, the "Natural killer cell mediated cytotoxicity" was the top term in KEGG analysis and the genes involved in this process were all targeted by miR-34c, which meant tumor immunology should not be ignored in NPC pathogenesis. Recent Fig. 4 MiR-34c levels during NPC progression and target genes validation. a The expression level of miR-34c-5p in human NPC specimens compared with control biopsy samples. b miR-34c-5p expression was higher in stage I, whereas stages II-IV had lower levels. c The protein expression levels of MET, CCND1, CDK6 and BCL2 in miR-34c mimic transfected CNE-2 cells were lower than the controls. Western blot was independently repeated at least three times. MicroRNA abundance was normalised to U6 RNA. Statistical analysis was performed using the t-tests (a, c) and the one-way ANOVA (b) studies showed miR-34c could affect NK cell's activation in melanoma [38], indicating an exciting future in the field of NK cell-based cancer immunotherapy [39].
MiRNAs regulatory network from IPA data showed the relationships between 5 target genes (TP53, CCND1, BCL2, CDK6, MET) and miR-34c in various cancers were experimental confirmed and located in the center of network (Fig. 3a). TP53, CCND1 and BCL2 were demonstrated to modulate tumor cells growth and proliferation via PI3K/AKT/mTOR signaling. These results indicated the downregulation of miR-34c might play an important role in NPC.
TP53 is predicted to be a diagnostic marker of NPC in IPA, which functions as a tumor suppressor protein that can induce cell cycle arrest and apoptosis. It was deemed to positively play a synergetic role in suppressing NPC tumorigenesis with miR-34 family not like miR-125a/b or miR-BHRF-1 which directly downregulated TP53 [40][41][42]. This further demonstrates that in the pathogenesis of NPC, even a single gene can be regulated by many factors, the exploration of miRNA regulatory network is especially vital in miRNA integrated analysis of NPC.
CCND1, CDK6 and BCL2 are all cell cycle regulators. CCND1 usually functions as a regulatory subunit of CDK6, whose activity is required for cell cycle G1/S transition and thought to be an oncogene in various tumors including NPC. The expression level of CCND1 was related to radiosensitivity in NPC treatment [43]. Previous studies showed EBV could encoded LMP1 to enhance the promoter activity of CCND1 in NPC, and miR-144 suppressed the expression of PTEN to increase the expression of CCND1 to promote tumor migration and invasion [44,45]. Our results indicated that CCND1 might be a target gene of miR-34c in NPC, which should be further confirmed in future. As one target of CCND1, CDK6 is important for cell cycle G1 phase progression and G1/S transition, while miR-26a can suppress CDK6 expression through EZH2 to inhibit tumorigenesis of NPC [46]. BCL2 acted as the anti-apoptotic protein by increasing the rate of cell division and arresting cells in the G0/G1 checkpoint of the cell cycle, and it could also be treated as a diagnosis marker as TP53 [47]. Mir-21 could target BCL2 directly to suppress NPC cells proliferation and migration, while PDCD4 regulated its expression level via modulating miR-184 [48,49].
The mechanism of MET in many tumors has been well explored in various research [50,51]. In NPC, miR-34c had been clearly proved to target MET [37], and it was activated by hepatocyte growth factor (HGF), then modulated cellular proliferation, migration, invasion via inhibition of EGR-1 expression [52,53]. Furthermore, PHA-665752, MET inhibitor, could furthermore effectively affect the radiosensitivity of NPC and suppress NPC cells proliferation [54]. Hence, Li et al. showed MET protein overexpression and gene amplification were independent prognostic factors and might be treated as therapeutic biomarkers in NPC [55].
Above all mentioned, the PI3K/AKT/mTOR signaling was considered as the function pathway of deregulated miRNAs to mediate cellular growth and proliferation, apoptosis and autophagy. CCND1, CKD6, BCL2 and TP53 worked together on cell circle regulation, while MET was also significant in NPC cellular proliferation, migration and invasion, and might work through other pathways. MiR-34c was thought to be the key factor in this network. And our qRT-PCR results further confirmed the expression level of miR-34c decreased with the ascent of tumor stage (Fig. 4b). So we hypothesized miR-34c regulated NPC development through the network constructed by these factors and inferred their mechanism in NPC respectively in this model (Fig. 5).

Conclusions
In conclusion, our study disclosed a global profiling of differentially expressed miRNAs that was closely related with NPC carcinogenesis and progression. For a better understanding of the miRNA regulation on genome-wide repression and activation, we use integrated approaches including GO, KEGG and IPA. In the miRNA regulatory network, miR-34c and its target genes including TP53, CCND1, BCL2, CDK6 and MET are shown to play