- Letter to the Editor
- Open Access
RNA-Seq profiling of circular RNAs in human laryngeal squamous cell carcinomas
Molecular Cancervolume 17, Article number: 86 (2018)
Abnormal expression of non-coding circular RNAs (circRNAs) have been reported in many types of tumors. circRNA have been suggested to be an ideal candidate biomarker for diagnostic and therapeutic implications in cancers. The aim of this study was to assess the circRNA expression profile of laryngeal squamous cell carcinomas (LSCC). The biopsy samples from patients with LSCC were obtained intra-operatively. The circRNA expression was performed using secondary sequencing. Among 10 patients with LSCC, 2 were well differentiated, 3 were moderately differentiated and 5 were adjunctive samples with normal and LSCC tissues. A total of 21,444 distinct circRNA candidates were detected. Among them, we defined the statistical criteria for selecting aberrant-expressed circRNA using a q-value of < 0.001 with a fold change of > 2.0 or < 0.5. A total of 29 circRNA were upregulated and 19 circRNA were downregulated significantly in the LSCC tissues. The intersection of these dysregulated circRNAs of normal-well differentiated set and normal-moderately differentiated set was then assessed to narrow the upregulated and downregulated circRNAs down to 18 and 5 respectively. Furthermore, an association of the circRNA-miRNA-mRNA was investigated, showing that 20 dysregulated circRNA successfully predicted an interaction with several cancer-related miRNAs. Finally, a further KEGG analysis showed that PPAR, Axon guidance, Wnt and Cell cycle signaling pathway were key putative pathways in the process of LSCC. hsa_circ:chr20:31876585–31,897,648 was found to be able to differentiate most of LSCC from the matching normal tissues. This observational study demonstrated dysregulation of circRNA in LSCC, which may have an impact on development of potential biomarkers in this disease. Validation of down-regulation of hsa_circ:chr20:31876585–31,897,648 in LSCC compared to each adjunctive tissue by Q-RT-PCR, indicating that hsa_circ:chr20:31876585–31,897,648 may be a novel promising tumor suppresser in LSCC.
Laryngeal squamous cell carcinomas (LSCC) is a malignant carcinoma derived from the epithelial tissue of laryngeal mucosa. It is an aggressive head and neck neoplasm, which makes up approximately 2.4% in all emerging malignant tumors worldwide every year . Moreover, there is an increasing incidence of LSCC in recent years. Multiply factors such as smoking and excessive alcohol consumption are contributing to the development of LSCC. Although recent research on LSCC have made some remarkable achievements, the etiology and pathogenesis underlying this disease development remains unclear. Over the last decade, despite of increasing numbers of research with a focus on identification of biomarkers to improve clinical outcomes of LSCC, LSCC is still associated with high morbidity and mortality. In particular, its prognosis remains very poor for patients with metastatic diseases, with a 5-year mortality rate of more than 50% . Therefore, better understanding of biomarkers associated with LSCC can play an important role in disease stratification and prognostication. It can also potentially inform the development of novel targeted therapies.
In recent years, identification of biomarkers in certain types of tumors have led to development of targeted therapies and improved clinical outcomes for biomarkers derived subsets of patients. To date, the dysregulated expression of non-coding RNAs including lncRNA, microRNAs have been closely associated with pathogenesis of some human cancers. Several non-coding RNA  have shown to be significantly changed in different sorts of cancers . Circular RNAs (circRNAs) are one of non-coding RNA discovered lately. They are produced by linking the 3′ and 5′ ends via the covalent bonds to form a closed loop, circRNAs. Due to this closed structure, circRNAs have been demonstrated to have a high level of stability and resistance to RNA degradation pathways . Therefore, they may be more promising and technically suitable biomarkers for cancers . Recent literatures have shown dysregulated expression of several circRNAs in human cancers . Moreover, resent studies have confirmed that the expression of circular RNAs is different between the LSCC sample and its adjacent tissues via microarray detection .
Despite several potential circRNA biomarkers have been detected in LSCC, due to the limitation of using microarray technology that could only focus on the known candidate circRNAs preset on the chip, a large number of novel promising circRNAs in LSCC beyond this microarray analysis are waiting for further screening. Moreover, to our knowledge, there are no studies to date have yet profiled the relationship of the circRNA-miRNA-mRNA expression in human LSCC in details. Therefore, the aim of this study was to profile the LSCC related dysregulated circRNA via secondary sequencing technology to identify novel circRNA biomarkers. Meanwhile, the relationship of these circRNA and their potential binding miRNA-mRNA were further analyzed to improve our understanding of the pathogenesis of this disease.
Collection and grouping of the LSCC samples based on the pathological classification
Specimens were collected from patients with laryngeal cancer who underwent laryngectomy at the Department of Otolaryngology, Head Neck Surgery, Beijing Friendship Hospital affiliated to Capital Medical University, China between August 2016 and December 2017. Secondary sequencing was used to test 10 matched samples of the LSCC tissues and corresponding adjacent non-neoplastic tissues. Ten matched cancerous and noncancerous tissues were tested by real-time qPCR. The study was approved by the Ethics Committee of Beijing friendship hospital. The patients who had received any cancer therapy before admission were excluded from this study. After immediate snap-frozen in liquid nitrogen post-resection, these specimens were stored at − 80 °C. Histopathological examination of the laryngeal tissue specimens was independently carried out by two pathologists. After histopathological vetting, a total of ten laryngeal samples–consisting of five LSCC (which included three moderately differentiated and two well differentiated specimens), and five matching contralateral normal tissue samples were finally included in the study.
RNA isolation and next generation RNA sequencing analysis
RNA isolation was performed according to standardized protocols in an RNA-dedicated work area with an access to RNase/DNase-free water and RNase-free labware. In brief, total RNA was extracted from the samples with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the kit’s instructions. The resulting RNA pellet was washed twice in 75% ethanol (1 ml), then air-dried, and re-suspended in RNase/DNase-free water (20 μl). After that, the Turbo DNase Kit (Ambion) was applied to the RNA solution to degrade any remaining DNA. The RNA was subsequently purified and concentrated with an RNeasy MinElute Cleanup kit (Qiagen, Valencia, CA, USA). The total RNA concentrations in the sample tissues and quality of each sample were then assessed with a NanoDrop ND1000 spectrophotometer (NanoDrop, Wilmington, DE, USA). Specifically, OD260/OD280 ratios between 1.8 and 2.1 deemed acceptable, while OD260/OD230 ratios of greater than 1.8 deemed acceptable. RNA integrity and DNA contamination were then assessed using electrophoresis on a denaturing agarose gel. In order to enrich pure circRNAs, linear RNAs were removed from the total RNA in each sample after treatment with RNase R (Epicentre, Inc.). The enriched circRNAs in each sample were amplified and fragmented into small pieces using divalent cations under an elevated temperature. RNA-seq library preparation was carried out using the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina. Sequencing of the libraries was conducted on an Illumina HiSeq™ 3000 system following the vendor’s recommended protocol (TruSeqR RNA Sample Preparation v2 Guide, Illumina Company Ltd).
HTSeq software was used to count the reads numbers mapped to each circRNA. RPM (Reads Per Million mapped reads, including TopHat mapping and TopHat-Fusion mapping)  was employed to calculate the expression level of individual circRNA. Differential expression between circRNAs was assessed by DEGseq algorithm. We defined the statistical criteria for selecting aberrant-expressed circRNA using a q-value < 0.05 with a fold change > 2.0 or < 0.5.
Quantitation of RNA using quantitative real-time reverse-transcription polymerase chain reaction (qRT-PCR)
Following RNA extraction from the LSCC and matching normal tissues (n = 10 pairs), cDNA was synthesized with a reverse transcriptase according to the kit’s instructions (SuperScript First-Strand Synthesis System for RT-PCR, Invitrogen, Carlsbad, CA, USA). circRNA expression was measured using qPCR (SYBR Green PCR Master Mix, Applied Biosystems, Foster City, CA, USA), conducted in a 20-μl reaction volume consisting of the following agents: 6.8 μl cDNA, 10 μl 2 x SYBR Green mix, 0.8 μl primer forward (5 μM), 0.8 μl primer reverse (5 μM), and 1.6 μl H2O. The qPCR reaction was then performed on an ABI 7300 Real-Time PCR System (Applied Biosystems): one 2-min cycle at 50 °C, one 10-min cycle at 95 °C, forty cycles of 15 s at 95 °C and 1 min at 60 °C. β-actin served as a control. PCR reactions were conducted three times. Relative circRNA expression was calculated using the 2 − ΔCt method.
Statistical analyses were performed using GraphPad Prism 5 (GraphPad Software, CA, USA), R software version 3.2.1 (http://www.r-project.org/) and Microsoft Excel (Microsoft, DC, USA). Student’s t-test was used to determine significance for differences in circRNA expression. A two-sided P value < 0.05 was considered statistically significant.
Prediction of circRNA-miRNA interactions
As circular RNAs were shown to have an impact on miRNA-mediated regulation of gene expression through miRNA sequestration , putative circRNA/miRNA interactions for the 12 differential circRNAs identified from the microarray and qRT-PCR validation experiments were predicted using the Arraystar miRNA target prediction software (Arraystar), based on TargetScan and miRanda algorithms . The Arraystar software was then used to search for MREs on the 12 differential circRNAs to construct the circRNA-miRNA network . After that, putative miRNAs were identified based on their seed-sequence complementarity using the circRNA-miRNA network. The parameters used for this Arraystar MRE analysis as previously described included : 1) seed type (seed-sequence complementarity) between nucleotides 2–7, 2) proximity of AU-richness to seed sequence, and 3) proximity to nucleotides pairing with the miRNA’s 3′ pairing sequence (nucleotides 13–16). The map of the circRNA/miRNA interaction network was then diagrammed using cytoscape 3.01.
Prediction of cancer-related circRNA-miRNA-target gene associations
Several online bioinformatics platforms were used in conjunction to predict circRNA-miRNA-target gene associations. First of all, miRNA candidates were identified in KEGG cancer-related pathways within the circRNA-miRNA network using the DIANA-miRPath v.3 platform, which has been reported in previous study , the candidate circRNA(s) were matched against the seed sequences of predicted target miRNA candidate(s) using the IntaRNA platform to predict interactions between RNA molecules. The settings that were applied to the DIANA-miRPath v.3 platform included: 0 suboptimal interactions, 5 minimum number of basepairs in seed, 0 maximum number of mismatches in seed, 37.0 temperature for energy computation, 22 folding window size target, and 22 max. Basepair distance (target). Finally, after validation of circRNA-target miRNA interaction, identification of the top ten-ranking cancer-related target genes for each of the target miRNA candidates were then performed using the DIANA-TarBase 7.0 database .
Characteristics of the human SLCC samples
Based on the histopathological vetting, a total of ten fresh tissue samples consisting of five LSCC samples (including three moderately differentiated and two well differentiated samples), and five matching adjacent normal tissue samples were finally detected in the study. The detailed information on patients’ characteristics and CONSORT flowchart describing the selection of tested tissue samples was provided in Additional file 1: Figure S1 and Additional file 1: Table S1.
Identification of differential circRNA profiles
The secondary sequencing was performed to profile circRNA expression in the laryngeal carcinoma samples. The sequencing results showed that 21,444 circRNA were detected, in which about 85% of reads were covered in the exon of genomic (Fig. 1a). After identification of differentially expressed circRNAs by using fold-change filtering (|log2(fold change)| > 1) and Student’s t-testing (Q-value < 0.05) (Fig. 1b), compared to the normal laryngeal tissue, there were 29 circRNAs that were significantly upregulated in LSCC and 19 circRNAs that were significantly downregulated in LSCC (Fig. 1d and e, Additional file 1: Table S2). Hierarchical clustering was then performed to demonstrate the circRNA expression patterns among the samples (Fig. 1c).
Cluster screening from the different group of the laryngeal carcinoma samples
Calculating the intersection of different expression circRNA from normal-well differentiated set and normal-moderately differentiated set, we were able to narrow candidate genes down to 23. Figure 1d and e was the Venn diagram of this intersection. Among them, 18 significantly downregulated circRNAs were overlapping between the foregoing comparisons, indicating the highest fold-changes (in descending order; Fig. 1d, Additional file 1: Table S3). Meanwhile, 5 significantly upregulated circRNAs that were overlapping between the foregoing comparisons were also identified (in descending order) (Fig. 1e, Additional file 1: Table S3).
Establishment of the circRNA-miRNA network
Due to interaction of circRNAs with miRNAs via miRNA response elements (MREs), putative MREs were searched for through Arraystar’s circRNA target prediction software. The predicted MREs for the 3 of 5 significantly-upregulated circRNAs were listed in Additional file 1: Table S3. As there were two unreliable predictions in upregulated circRNAs, we were unable to match them to higher frequency of targeting. The predicted MREs for the 17 significantly downregulated circRNAs were also listed in Additional file 1: Table S3. Based on these dysregulated circRNAs and their predicted MREs, a network map of circRNA-miRNA interactions with Cytoscape was established (Fig. 2a).
Prediction of the characteristics of the circRNA signature related to LSCC using gene ontology and KEGG analyses
Current literature has shown that the circRNA serves as miRNA sponges in some circumstances and regulate gene expression. Based on this, we ranked and identified miRNAs for LSCC related circRNAs. We included the top five ranking candidate miRNAs for each circRNA through specific base pairing. Moreover, we investigated the functions of their target genes using Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. It demonstrated that the target genes that were related to this circRNA signature participated in various biological processes, such as developmental process, cellular regulation and cell junction (Additional file 1: Figure S2). The altered circRNAs could also have an impact on several vital pathways, such as MAP kinase (MAPK) signaling pathway, Phosphatidylinositol 3-Kinase–AKT (PI3K-Akt) signaling pathway, and Ras signaling pathway (Fig. 2b and Additional file 1: Figure S3).
Prediction of cancer-related circRNA-miRNA-target gene associations
To better understand the underlying molecular mechanism(s) of the differential circRNAs, we first selected the miRNA from KEGG analysis that enriched in pathway in cancer (Fig. 2b, Additional file 1: Figure S3) and used the bioinformatics platform to profile the cancer-related miRNA network. Then, the predicted cancer-related circRNA-miRNA associations were showed in Fig. 3a. Additional file 1: Table S3 indicated the number of each dysregulated circRNA that were interacted with cancer-related miRNA was annotated in Fig. 3a (green nodes). Finally, we selected the top three circRNA (hsa_circ:chr17:48268229–48,277,287, hsa_circ:chr20:31876585–31,897,648, hsa_circ:chr 1:207103173–207,105,911) which could bind to 151, 112 and 208 cancer related miRNA, respectively (Additional file 1: Table S3). By comparing and overlapping the lists of these three top-ranking cancer-related circRNAs, several promising cancer-related pathways were able to be identified, which may be potential targets of these circRNA/miRNA/ mRNA axis in LSCC (Fig. 3a).
Top three candidate circRNA that were involved in the main regulation of the PPAR, AXON-guidance, Wnt and cell cycle signaling pathway in LSCC development
Three candidate circRNAs (207,105,911, 31,897,648 and 48,277,287) were selected, and their corresponding numbers of cancer related miRNAs were 151, 112 and 208, respectively. Interestingly, nine common targeted cancer related miRNAs were found in these three candidate circRNAs (Fig. 3b). Further KEGG analysis were performed using these nine miRNA-target genes (Fig. 3c), which elucidated that PPAR, Axon guidance, Wnt and Cell cycle signaling pathway were found to be a key putative signaling pathway in the process of LSCC.
At last, three candidate circRNAs were prepared for validation. As no specific primer can be designed for detecting hsa_circ:chr17:48268229–48,277,287 and hsa_circ:chr 1:207103173–207,105,911, only hsa_circ:chr20:31876585–31,897,648 was selected for qRT-PCR validation. T1 (the tumor tissues) and C1 (control, the adjunct normal tissues) were a paired sample from the first patient. All paired samples from ten patients also showed the downregulation of hsa_circ:chr20:31876585–31,897,648 in the LSCC tissues (T1-T10) compared to each of their adjunct tissues (C1-C10) (Fig. 3d). It demonstrated that hsa_circ:chr20:31876585–31,897,648 was significantly differentiated between most of the ten paired LSCC tumor tissues and the matching normal tissues. It indicated that circRNA hsa_circ:chr20:31876585–31,897,648 may have the potential to be a novel biomarker for LSCC in future research.
In this study, we have profiled the circRNA expression in LSCC in details using secondary sequencing to improve our understanding of the pathogenesis of LSCC. We have also explored a potential role of this novel circRNA biomarker used for the diagnostic and prognostic value of the disease. From sequencing and subsequent overlapping analysis, a total of 20 significantly candidate differentiated circRNAs in LSCC were identified, including 16 upregulated and four downregulated, compared to the adjacent normal tissues. Four significantly upregulated circRNAs were identified as the most potential biomarkers related in LSCC.
Furthermore, a network map of circRNA-miRNA interactions was also established for total 20 significantly differentiated circRNAs using a set of bioinformatics tools. The functions of their target genes were then assessed with the Gene Oncology (GO) analysis. It demonstrated that the target genes that were related to these 20 dysregulated circRNA participated in various biological processes, such as growth factors binding, protein phosphatase activator activity, skin development, wound healing, most of which were related to proliferation process. Moreover, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed to manifest some vital pathways, including Ras, MAPK, PI3K-AKT cancer related pathways, involved in the altered circRNAs in the LSCC samples.
Previous molecular research on LSCC has primarily focused on genomic analyses. In recent years, researchers have also assessed changes in the circRNA expression in LSCC via microarray . However, as the relationship of circRNA-miRNA-mRNA-based research is an emerging area of investigation . to our knowledge, there is no study to date has yet assessed circRNA in LSCC based on these relationships. Therefore, our study comprehensively profiled the relationship of 20 candidate circRNAs and their possible interactions with miRNAs. Based on a hypnosis that the more important circRNAs could interact with more different cancer related miRNAs and further regulated different cancer related pathways. Hence the top three dysregulated circRNAs - that have been predicted could bind nine common cancer related miRNAs in LSCC - were further validated by qRT-PCR. Then the common miRNA targeted 576 mRNA were analyzed by KEGG. The KEGG annotation indicated that these three circRNA may mainly regulate (involve in) the PPAR, Axon guidance, Wnt and Cell cycle signaling pathway. Therefore, we presumed that these pathways may be potential therapeutic target in LSCC.
Our study has several limitations. First of all, as a single centre study, it has a relatively small sample size, which has an impact on the statistical power. Further larger studies are needed to assess the role of circRNAs in LSCC. Secondly, the association of circRNAs and severity of LSCC was unable to be assessed in this study due to the limited number of LSCC samples, which also warrants further investigations.
In conclusion, the high-throughput transcriptome sequencing and bioinformatics analysis of the laryngeal squamous cell carcinoma provides useful diagnostic markers and potential therapeutic targets for the future research of this disease.
Kyoto Encyclopedia of Genes and Genomes
Laryngeal squamous cell carcinomas
Christensen A, Kristensen E, Therkildsen MH, Specht L, Reibel J, Homoe P. Ten-year retrospective study of head and neck carcinoma in situ: incidence, treatment, and clinical outcome. Oral Surg Oral Med Oral Pathol Oral Radiol. 2013;116(2):174–8.
Shah JP, Karnell LH, Hoffman HT, Ariyan S, Brown GS, Fee WE, et al. Patterns of care for cancer of the larynx in the United States. Arch Otolaryngol Head Neck Surg. 1997;123(5):475–83.
Lee KM, Choi EJ, Kim IA. microRNA-7 increases radiosensitivity of human cancer cells with activated EGFR-associated signaling. Radiother Oncol. 2011;101(1):171–6.
Zhang XO, Wang HB, Zhang Y, Lu X, Chen LL, Yang L. Complementary sequence-mediated exon circularization. Cell. 2014;159(1):134–47.
Li F, Zhang L, Li W, Deng J, Zheng J, An M, et al. Circular RNA ITCH has inhibitory effect on ESCC by suppressing the Wnt/beta-catenin pathway. Oncotarget. 2015;6(8):6001–13.
Bourguignon LY, Earle C, Wong G, Spevak CC, Krueger K. Stem cell marker (Nanog) and Stat-3 signaling promote MicroRNA-21 expression and chemoresistance in hyaluronan/CD44-activated head and neck squamous cell carcinoma cells. Oncogene. 2012;31(2):149–60.
Jeck WR, Sharpless NE. Detecting and characterizing circular RNAs. Nat Biotechnol. 2014;32(5):453–61.
Xuan L, Qu L, Zhou H, Wang P, Yu H, Wu T, et al. Circular RNA: a novel biomarker for progressive laryngeal cancer. Am J Transl Res. 2016;8(2):932–9.
Sand M, Bechara FG, Gambichler T, Sand D, Bromba M, Hahn SA, et al. Circular RNA expression in cutaneous squamous cell carcinoma. J Dermatol Sci. 2016;83(3):210–8.
Peng N, Shi L, Zhang Q, Hu Y, Wang N, Ye H. Microarray profiling of circular RNAs in human papillary thyroid carcinoma. PLoS One. 2017;12(3):e0170287.
The authors would like to thank to Prof. Weiwei Guo, the General Hospital of People’s Liberation Army (301 hospital), China for constructive criticism of the manuscript. Special thanks also to Mr. Lei, Zhang for the provision of immunohistochemistry techniques.
This project was supported by the Scientific Research Common Program of Beijing Municipal Commission of Education (KM201510025028), Research Foundation of Beijing Friendship Hospital, Capital Medical University (No. yydszx2015–02 and yyqdkt2014–23). National Nature Science Foundation of China (81470684),. Clinical Special Fund of Jiangsu Province (b12014032).
Availability of data and materials
The datasets obtained and/or analyzed during the current study were available from the corresponding authors in a reasonable request.
Ethics approval and consent to participate
The human cancer tissues used in this study were approved by the institute ethical committee of Beijing Friendship Hospital, Capital Medical University. (Approval Number: YYBB-B01–01-R02-V4.0 / 20,170,605).
Consent for publication
We have received consents from individual patients who have participated in this study. The consent forms will be provided upon request.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. Demographic characteristics of patients with LSCC involved in this study. Table S2. Individual circRNAs detected in the study (please see the attached excel spreadsheet). Table S3. LSCC specific circRNAs detected in the study (please see the attached excel spreadsheet). Figure S1. The representative images of H&E-stained normal laryngeal mucosal and LSCC specimens. Pictures included well differentiated (upper), moderately differentiated LSCC (middle), and normal tissues (lower). From left to right, image magnifications of 40×, 100×, 200×, and 400 × were displayed (scale bar = 100 μm). Figure S2. Gene Ontology annotation analysis for 20 circRNA interacted miRNA and their target gene related significant enriched biological process, cellular components and molecular function. Figure S3. KEGG analysis showing the map of pathway in LSCC using the dysregulated circRNA-miRNA-target genes. (ZIP 8463 kb)