Integrated analysis of lncRNA-miRNA-mRNA ceRNA network in squamous cell carcinoma of tongue

1Head and neck squamous cell carcinoma (HNSCC), which is a disease that causes serious harm to humans, is highly correlated with alcohol consumption, tobacco smoking, and betel nut chewing, and human papillomavirus infection. Squamous cell carcinoma of the tongue (SCCT) is a particular subtype and the main cause of patient mortality and morbidity from HNSCC. In general, the clinical features and treatment strategies for SCCT are similar to those of other HNSCCs, with surgical resection being the primary treatment choice. However, due to late diagnosis of locally advanced malignancies, in many cases of SCCT, surgery is either no longer an option, or should be avoided to maintain the patient’s quality of life. Despite the advances in treatment options, the prognosis of patients with advanced SCCT remains poor. In China, although pingyangmycin and/or cisplatin-based chemotherapies have shown good results, chemotherapy resistance always develops later and causes the therapy to fail. In the past three decades, the 5-year survival rate of patients with SCCT was less than 50%. Therefore, the main goal of our research has been to obtain more knowledge about SCCT cells and to identify novel therapeutic targets for treating the disease.

2Long non-coding RNAs (lncRNAs), which do not have protein-coding functions, have recently attracted increasing research attention. These RNAs play a significant role in different cellular processes, particularly in numerous kinds of tumors. For example, lncRNAs can act as biomarkers for the prognosis and diagnosis of lung adenocarcinoma. MicroRNAs (miRNAs) are small, endogenous, non-coding RNAs composed of 19–25 nucleotides. They exert the important function of regulating gene expression, and their regulatory networks are involved in many biological processes. In 2011, Salmena et al. proposed the competitive endogenous RNA (ceRNA) hypothesis, which was subsequently supported by several lines of evidence. This hypothesis describes the competitive activity of some RNAs (as ceRNAs) for common binding sites of target miRNAs, thereby altering the function of the target miRNA. The core concept is that ceRNAs interact with target miRNAs through miRNA response elements to control the transcriptome on a large scale. In the past several years, lncRNAs and SCCT were confirmed to be closely related. For instance, expression of the lncRNA SNHG6 is significantly increased in tongue cancer, and interference with SNHG6 expression can inhibit the proliferation and epithelial–mesenchymal transition (EMT) of tongue cancer cells. Zhang et al. found that the oncogenic lncRNA KCNQ1OT1 plays a vital role in SCCT growth and chemoresistance, and can be used as a new target for SCCT treatment. However, previous studies had focused on the mechanism of a single lncRNA-miRNA-mRNA axis, and there is currently no reported ceRNA network in SCCT. Consequently, it is extremely important to investigate the role of ceRNA networks in the poor prognosis of SCCT. By further learning how lncRNAs function in the pathogenesis of SCCT, we may find solutions to the most pressing challenges faced in treating this disease.

3In this study, the mRNA, miRNA, and lncRNA expression profiles of SCCT and normal tissues were downloaded from The Cancer Genome Atlas (TCGA). In addition, through comprehensive analysis, the ceRNA network for SCCT was builted, which will serve to find new targets and pathways for the development of treatments to prolong patient survival times. Finally, we conducted a prognostic analysis with several important lncRNAs and found a biomarker that could predict survival in patients with SCCT.

4The SCCT cases data of clinical and RNA expression were collected from TCGA database. The exclusion criteria were including: (i) histological diagnosis was not SCCT; (ii) no complete data (including gender, age, survival status, stage, and survival time) for analysis. 118 SCCT patients were enrolled in the study. The number of patients aged < 68 years was 49, 69 patients were ≥ 68 years old. 43 patients were female and 75 patients were male. The number of stage I, II, III, IVa and IVb patients were 13, 20, 27, 56 and 2, respectively. The number of patients, who were white, Asian, black or African American and not available, were 106, 5, 5 and 2, respectively. The number of patients, who were hispanic or latino, were 104. 9 patients were not hispanic or latino, and 5 patients were not reported. 46 patients were dead, and 72 patients were alive. SCCT characteristics and clinical data of the patients are showed in Table 1 and Additional file 3: Table S1.

5RNA expression data of SCCT patients were available from TCGA database. The raw reads of lncRNA and mRNA data were post-treated and normalized in R software (Additional file 1: Figure S1). The miRNA expression data from TCGA database were normalized in R software (Additional file 2: Figure S2). The tumor tissue and adjacent non-tumor tissue of SCCT patients were facilitated differential expressions of mRNA, lncRNA, and miRNA. Furthermore, intersection of lncRNA, miRNA and mRNA was selected.

6Compared to the normal group with SCCT, “limma” package in R software was used to identify the differentially expressed mRNAs (DEmRNAs) with thresholds of |fold Change (FC)| > 2.0 and P value < 0.01 and differentially expressed miRNAs with |FC| > 2.5 and P value < 0.01.

7“ClusterProfiler” package in R software was used for functional enrichment analysis, and GO biological processes and KEGG pathways at the significant level (q-value < 0.01) were employed.

8To determine the prognostic characteristics of DERNAs, combining the clinical data the survival curves of these samples with differentially expressed mRNA, lncRNA and miRNA were plotted by using the “survival” package in R based on Kaplan-Meier curve analysis. P values < 0.05 were regarded as significant.

9The lncRNA-miRNA-mRNA ceRNA network was based on the theory that lncRNAs can directly interact by invoking miRNA sponges to regulate mRNA activity. “GDCRNATools” (http://bioconductor.org/packages/devel/bioc/html/GDCRNATools.html) package in R software were used to establish ceRNA network. The ceRNA network was plotted with Cytoscape v3.6.0. The plug-in BinGO of Cytoscape is an APP for BF network of the hub genes.

10We explored 1943 SCCT-specific mRNAs (1007 downregulated and 936 upregulated; Table 2 and Fig. 1) and 107 lncRNAs (34 downregulated and 73 upregulated; Fig. 1, Table 2, and Table 3). The differentially expressed genes (DEGs) are shown in Fig. 2a. Additionally, 100 miRNAs (44 upregulated and 56 downregulated; Fig. 2b, and Table 4) were found.

11GO analysis results showed that changes in biological processes (BP) of DEGs were significantly enriched in extracellular structure organization, extracellular matrix organization, urogenital system development, muscle contraction, collagen metabolic process, mitotic nuclear division, renal system development, collagen catabolic process, sister chromatid segregation, and collagen metabolic process (Fig. 3a). Changes in cell component (CC) of DEGs were mainly enriched in proteinaceous extracellular matrix, endoplasmic reticulum lumen, apical part of cell, contractile fiber, myofibril, contractile fiber part, sarcomere, extracellular matrix component, basement membrane, basal lamina (Fig. 3b). Changes in molecular function (MF) were mainly enriched in actin binding, growth factor binding, coenzyme binding, microtubule binding, iron ion binding, glycosaminoglycan binding, collagen binding, structural constituent of muscle, extracellular matrix structural constituent, platelet−derived growth factor binding (Fig. 3c). KEGG pathway analysis revealed that the DEGs were mainly enriched in focal adhesion, human papillomavirus infection, ECM − receptor interaction, protein digestion and absorption, small cell lung cancer, arginine and proline metabolism, PI3K-Akt signaling pathway, dilated cardiomyopathy (DCM), valine, leucine and isoleucine degradation, cell cycle (Fig. 4).

12We studied the association of the DEGs and DEmRNAs with patient’ survival to identify the key genes and mRNAs that were related to the prognosis of patients with SCCT. We identified 10 mRNAs (CSRP2, CKS2, ADGRG6, MB21D1, GMNN, RIPOR3, RAD51, PCLAF, ORC1, NAGS), 9 lncRNAs (LINC02560, HOXC13 − AS, FOXD2 − AS1, AC105277.1, AC099850.3, STARD4 − AS1, SLC16A1 − AS1, MIR503HG, MIR100HG) and 8 miRNAs (miR − 654, miR − 503, miR − 450a, miR − 379, miR − 369, miR − 190a, miR − 101, let−7c) that were significantly differentially expressed in the survival analyses (Fig. 5a-c).

13We built the ceRNA network on the basis of the miRNA, lncRNA, and mRNA the expression profiles in patients with SCCT. In total, 27 miRNA nodes, 53 mRNA nodes, 6 lncRNA nodes, and 152 edges were identified as differentially expressed profiles. The network is showed in Fig. 6. It is well known that lncRNAs and mRNAs have co-expression patterns in ceRNA networks. Thus, we chose a hub lncRNA (degree> 5, Additional file 4: Table S2) and its linked mRNAs and miRNAs in the triple global network and then reconstructed the sub-network. As shown in Fig. 7, the lncRNA KCNQ1OT1-miRNA-mRNA sub-network was composed of 1 lncRNA node, 7 miRNA nodes, 11 mRNA nodes, and 41 edges.

14SCCT, a major type of HNSCC, is a refractory cancer under current therapeutics. Studies have demonstrated that lncRNAs regulate gene expression through a variety of pathways, contributing to tumorigenesis and tumor metastasis. The ceRNA hypothesis proposes a new regulatory mechanism mediated by lncRNAs that are used as endogenous miRNA sponges. In this study, we found the genes and mRNAs that were differentially expressed between normal and tumor tissue. Through GO and KEGG analyses, we further analyzed the pathways and functions in which the DEGs are involved. The GO biological processes results suggested that specific genes may be concentrated in several process areas, such as extracellular structures, muscle contraction, and mitotic nuclear division. Some of the annotated pathways have been shown to be associated with cancer in previous studies. PI3K-Akt signaling is involved in cell proliferation and growth as well as down-regulating cell apoptosis. Recent preclinical and clinical studies of highly selective agents that target various regulators of the mammalian cell cycle have demonstrated cell-cycle arrest in models of human cancer. Through survival analysis, we identified 10 mRNAs (CSRP2, CKS2, ADGRG6, MB21D1, GMNN, RIPOR3, RAD51, PCLAF, ORC1, NAGS), 9 lncRNAs (LINC02560, HOXC13 − AS, FOXD2 − AS1, AC105277.1, AC099850.3, STARD4 − AS1, SLC16A1 − AS1, MIR503HG, MIR100HG) and 8 miRNAs (miR − 654, miR − 503, miR − 450a, miR − 379, miR − 369, miR − 190a, miR − 101, let−7c) that were significantly related to the overall survival of patients with SCCT. Next, by using bioinformatics tools, we builted a ceRNA network with SCCT-specific miRNA and lncRNA expression and selected the hub lncRNA KCNQ1OT1 to construct a sub-network.

15KCNQ1OT1, also known as KCNQ1 overlapping transcript 1, is an imprinted antisense lncRNA in the KCNQ1 locus. Early studies have shown that KCNQ1OT1 is up-regulated and involved in the tumorigenesis of breast cancer and hepatocellular carcinoma. Zhang et al. found that KCNQ1OT1 could induce SCCT cell growth and inhibit the sensitivity of the tumor to cisplatin. Previous studies have shown that KCNQ1OT1 acts as an oncogene and plays a key role in promoting SCCT cell growth and chemotherapy resistance.

16We constructed a SCCT-specific ceRNA network and chose a hub lncRNA for SCCT by bioinformatics analysis. To the best of our knowledge, only a limited number of studies have analyzed lncRNA obtained from large-scale samples. We provide a method for identifying potential lncRNA biomarkers. Furthermore, we found the ceRNA network in SCCT, which should help further our understanding of the mechanism underlying the pathogenesis of this disease.