Laryngeal squamous cell carcinoma (LSCC) is one of the common malignancies of the upper respiratory tract that has been associated with a deterioration of the environment and an increase in the occupational stress. It was estimated that 13,360 new cases were diagnosed in 2017 in the United States, of which over 3,660 were fatal (Siegel, Miller & Jemal, 2017). In China, an estimated 26,400 new cases of LSCC and 14,500 cancer-related deaths also occurred in 2015 (Chen et al., 2016). Although patients with LSCC can be managed by surgical intervention, radiation therapy and chemotherapy, the overall five-year survival remains poor (approximately 60%) (Rudolph et al., 2011). Therefore, there is an urgent need to deeply understand the molecular mechanisms underlying LSCC carcinogenesis or progression in order to develop more effective therapeutic strategies.
Accumulating evidence has suggested that non-coding RNAs (ncRNAs) play crucial roles in the initiation and development of tumors. ncRNAs are loosely categorized into small ncRNAs and long non-coding RNAs (lncRNAs), both of which have regulatory functions in various biological processes. The well-documented small ncRNAs are microRNAs (miRNAs; ∼22 nucleotides long) that regulate gene expression by binding to complementary sequences in the 3′ untranslated region (UTR), leading to either inhibition of translation or degradation of the transcripts (Jean & Mihaela, 2014). Although the mechanisms remain unclear, growing evidence supports that lncRNAs could function as competing endogenous RNAs (ceRNAs) by competitively binding to miRNAs through their miRNA response elements (MRE) and subsequently regulate target RNA expression (Salmena et al., 2011). This ceRNA mechanism has generated much interest to explain tumor development and progression in many malignancies, such as gastric cancer (Song et al., 2018), thyroid carcinoma (Zhao et al., 2018) and hepatoblastoma (Liu et al., 2017a). Recent studies also have preliminarily revealed several underlying ceRNA regulatory interactions in LSCC. Luciferase reporter assay and Western blotting results suggested that AC026166.2-001 could act as a sponge of miR-24-3p and regulate the expression of p27 and cyclin D1 (Shen et al., 2018). lncRNA H19 was shown to serve as a ceRNA by sponging miR-148a-3p to upregulate the target gene DNA methyltransferase 1 (Wu et al., 2016). NEAT1 was also reported to regulate the expression of cyclin dependent kinase 6 through modulating miR-107 (Wang et al., 2016). Furthermore, a ceRNA network, including 30 genes, 21 miRNAs and 19 lncRNAs was also built based on microarray analysis of 6-paired clinical samples in LSCC (Zhang et al., 2016). However, analysis of the lncRNA-miRNA-mRNA regulatory network of LSCC with larger sample sizes and confirmation of their clinical associations are still lacking.
In addition, DNA methylation has been identified as an important mechanism to regulate gene expression in cancer cells epigenetically, which not only regulates the expression of protein-encoding genes, but also affects miRNAs and lncRNAs. For example, hyper-methylation of the promoter region was observed to lead to a loss of expression of lncRNA maternally expressed gene 3 (MEG3). Downregulated MEG3 was insufficient to sponge miR-9 and block its inhibition effects on the expressions of E-cadherin and FOXO1, consequentially resulting in poor prognosis in patients with esophageal squamous cell carcinoma (Dong et al., 2017). The study of Guo et al. (2018a) also suggested lncRNA CTC-276P9.1 was hyper-methylated in esophageal squamous cell carcinoma. Over-expression of CTC-276P9.1 inhibited cancer cell proliferation and invasion in vitro probably by regulating epithelial-mesenchymal transition. Liao et al. (2015) identified 761 lncRNA genes with DNA hyper-methylation in colorectal cancer using a free MethylCap-seq dataset. Cheung & Lee (2010) found that the loci of three miRNAs (namely miR-199a-2, miR-124a-2 and miR-184) were linked to hyper-methylated differentially methylated regions (DMRs) in human testicular cancer. However, the DNA methylation regulatory mechanisms of miRNAs and lncRNAs have rarely been reported in LSCC.
The goal of this study was to establish an lncRNA-miRNA-mRNA ceRNA network in LSCC using larger samples and to investigate their methylation patterns. Our results may provide new clues for biologists to further understand the pathogenesis of LSCC.
Material and Methods
lncRNA, miRNA, mRNA and methylation data were retrieved from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) in January 2018 according to the following inclusion criteria: (1) lncRNA, miRNA, mRNA expression or methylation profiles; (2) laryngeal tissue samples, not blood, interstitial fluid or cells; (3) inclusion of control; (4) human samples; and (5) patients with LSCC.
Two lncRNA microarray datasets were obtained under accession number GSE59652 (7 LSCC and 7 paired adjacent normal tissues) (Shen et al., 2014) and GSE84957 (9 LSCC and 9 paired adjacent non-neoplastic tissues) (Feng et al., 2016). The microarray platforms of GSE59652 and GSE84957 were Agilent-033010 (GPL13825, Arraystar Human LncRNA microarray V2.0) and Agilent-042818 (GPL17843, Agilent-042818 Human lncRNA Micorarray 8_24_v2), respectively.
Two miRNA microarray datasets were collected under accession number GSE70289 (12 LSCC tissues and 4 adjacent normal tissues) (Karatas et al., 2015) and GSE62819 (5 LSCC carcinoma and 5 paired adjacent non-neoplastic tissues). The microarray platforms of GSE70289 and GSE62819 were Agilent-031181 (GPL15018, Unrestricted_Human_miRNA_V16.0_Microarray 030840) and Affymetrix Multispecies miRNA-3 Array (GPL16384), respectively.
Four mRNA microarray datasets were available under accession number GSE51985 (10 LSCC and 10 paired adjacent normal tissues), GSE84957 (9 LSCC and 9 paired adjacent normal tissues) (Feng et al., 2016), GSE59102 (29 LSCC and 13 normal margin tissues) and GSE58911 (15 LSCC and 15 normal tissue distant to LSCC) (Sharon et al., 2015). The microarray platforms of GSE51985, GSE84957, GSE59102 and GSE58911 were Illumina HumanHT−12 V4.0 (GPL10558), Agilent-042818 (GPL17843, Human lncRNA Micorarray 8_24_v2), Agilent-014850 (GPL6480, Whole Human Genome Microarray 4x44K G4112F) and Affymetrix Human Gene 1.0 ST Array (GPL6244), respectively.
One set of DNA methylation data was acquired under accession number GSE25093 (Poage et al., 2012; Poage et al., 2011) which included 213 blood and 109 tissue samples. Among the 109 tissue samples, 56 were isolated from oral, 16 from pharyngeal, and 22 from laryngeal origin, while 15 were of unclear origin. Thus, only these 22 samples from laryngeal origin (15 LSCC tissues and 7 controls) were used in our study. The microarray platform of GSE25093 was Illumina HumanMethylation27 BeadChip (GPL8490, HumanMethylation27_270596_v.1.2).
The mRNA and miRNA Seq-data of head and neck squamous cell carcinoma (Level 3) were also downloaded from The Cancer Genome Atlas (TCGA; https://tcga-data.nci.nih.gov/). After sample barcode screening, 559 were miRNA-mRNA matched samples, of which 18 were distributed in the alveolar crest, 30 in the root of the tongue, 22 in the buccal mucosa, 67 in the mouth floor, 8 in the hard palate, 9 in the laryngeal pharynx, 138 in the larynx, 3 in the lip, 38 in the oral cavity, 156 in the tongue, 10 in the oropharynx, 45 in the tonsil and 15 from an unclear location. Only the 138 samples from the larynx were used in our study.
For the data from Affy platform, the raw data in CEL. files were preprocessed using the oligo package (version 1.41.1; Carvalho & Irizarry, 2010) in R (version 3.4.1; R Development Core Team, 2017), including data transformation, missing value imputation with median, background correction with MAS method and quantile normalization.
For the data from Agilent and Illumina platforms, the raw data in TXT. files were preprocessed using the Linear Models for Microarray Data (LIMMA) package (version 3.34.0; Ritchie et al., 2015) in R, including data log2 transformation and median normalization.
The data (FPKM, fragment per kilobase per million mapped reads) from TCGA were quantile normalized using the preprocessCore package (version 1.40.0; Bolstad, 2019) in R.
Differential expression analysis
The differentially expressed lncRNAs (DELs) and miRNAs (DEMs) between LSCC and normal controls were identified using the LIMMA method in R from their two included microarray datasets (lncRNA: GSE59652 and GSE84957; miRNA: gSE70289 and GSE62819). The p-value <0.05 and |logFC(fold change) | > 0.263 were set as the cut-off points. The overlap in the above two datasets was used for the following analysis of lncRNAs and miRNAs, respectively.
The differentially expressed genes (DEGs) between LSCC and normal controls were identified using the MetaDE.ES function in MetaDE package (version 1.0.5, https://cran.r-project.org/web/packages/MetaDE/) of R from its four included microarray datasets (GSE51985, GSE84957, GSE59102 and GSE58911). The p-value <0.05 and false discovery rate (FDR) <0.05 were set as the cut-off points. The DEGs with the same expression trend (tau2 statistic = 0, p-value of Chi-square based Q-test >0.05) in the four datasets were selected for the following analysis.
Wilcoxon signed-rank test (http://www.bioconductor.org/help/search/index.html?q=wilcox.test/) was used to screen the DMRs between LSCC and normal controls. P < 0.05 was set as the threshold value. Human annotation data were retrieved from GENCODE Release 19 (GRCh37.p13) (https://www.gencodegenes.org/human/release_19.html). The sequences of miRNAs, lncRNAs and mRNAs in the corresponding platform GPL8490 were blasted with the GRCh37.p13 to obtain the differentially methylated miRNAs, lncRNAs and mRNAs, which were then overlapped with the DELs, DEMs and DEGs to screen methylated-related DELs, DEMs and DEGs, respectively.
CeRNA regulatory network construction
Three reliable online databases, including miRcode (version 11; http://www.mircode.org/), starBase (version 2.0; http://starbase.sysu.edu.cn/index.php) (Li et al., 2014) and DIANA-LncBase (version 2.0; http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=lncbasev2/index-predicted) (Paraskevopoulou et al., 2013) were used to screen the interactions between lncRNAs and miRNAs. The union of these three datasets was used for the following analysis. The target genes of miRNAs that were linked to the lncRNAs were predicted using four frequently used algorithms, including TargetScan (version 7.2; http://www.targetscan.org/vert_71/) (Agarwal et al., 2015), miRBase (version 22; https://www.ebi.ac.uk/enright-srv/microcosm/htdocs/targets/v5/) (Griffithsjones et al., 2005), miRanda (version 1.9; http://www.microrna.org/microrna/home.do/) (John et al., 2005) and miRTarBase (version 7.0; http://mirtarbase.mbc.nctu.edu.tw/php/index.php) (Chou et al., 2017). The target genes predicted by at least two databases and a negative association with miRNAs were retained. The lncRNA-miRNA and miRNA-mRNA interactions were integrated to construct the ceRNA network, which was visualized using Cytoscape software (version 3.4; Shannon et al., 2001-2008) (Kohl, Wiese & Warscheid, 2011).
Function enrichment analysis
The Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool (version 6.8; http://david.abcc.ncifcrf.gov) (Da, Sherman & Lempicki, 2009) was used for Gene Ontology (GO) terms [including molecular function (MF), biological process (BP) and cellular component (CC) categories] and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analyses of genes in the ceRNA network. P-value <0.05 was set as the cut-off value.
Clinical associations of lncRNAs, miRNAs and mRNAs in the ceRNA network
The expression levels of lncRNAs, miRNAs and mRNAs in the ceRNA network were downloaded from the TCGA data. Univariate Cox regression analysis was performed to screen for the prognosis-related (including overall survival, OS; and recurrence-free survival, RFS) lncRNAs, miRNAs and mRNAs using the survival package (version 2.40.1; https://cran.r-project.org/package=survival), which was used to construct the prognosis-related ceRNA network. The samples were divided into two groups based on the expression of each lncRNA, miRNA and mRNA: a low expression group (<median) and a high expression (>median) group. The Kaplan–Meier method with the log-rank test was used to estimate the difference in OS and RFS between the high and low expression groups. P < 0.05 was considered statistically significant. Furthermore, multivariate Cox regression analysis was also performed using the survival package (version 2.40.1; Therneau & Lumley, 2019) to evaluate the prognostic independence of lncRNAs, miRNAs and mRNAs. The association of nodes in the prognosis-related ceRNA network with other clinical characteristics was also analyzed using the multiple linear regression model (https://stat.ethz.ch/R-manual/R-patched/library/stats/html/lm.html) in R.
Differential expression analysis
The data analysis workflow is displayed in Fig. 1. After data normalization (Supplemental Information 1–8), the DELs, DEMs and DEGs between LSCC and normal samples were screened according the stated thresholds. The results showed 306 (156 downregulated and 150 upregulated) and 396 (252 downregulated and 144 upregulated) DELs were identified in the datasets of GSE59652 (Fig. 2A) and GSE84957 (Fig. 2B) (Supplemental Information 9), respectively. After comparison, 40 DELs were found to be shared in these two datasets, including six upregulated and 20 downregulated with the consistent expression trend (Fig. 3A) (Supplemental Information 9); a total of 1,307 (765 downregulated and 542 upregulated) and 491 (126 downregulated and 365 upregulated) DEMs were identified in the datasets GSE62819 (Fig. 2C) and GSE70289 (Fig. 2D), respectively (Supplemental Information 9). After comparison, 443 DEMs were found to be common in these two datasets, among which 152 upregulated and 63 downregulated DEMs were shown to have a consistent expression trend (Fig. 3B) (Supplemental Information 9); 2,975 DEGs were found to display the similar expression trend in four mRNA expression profiles GSE51985, GSE84957, GSE59102 and GSE58911 (Fig. 2E) (Supplemental Information 9); and 4,567 DMRs were identified in the LSCC genome of GSE25093 dataset, including 1616 hypomethylated and 2,951 hypermethylated (Fig. 2F) (Supplemental Information 9). After GENCODE annotation and blast analysis, 122 lncRNAs, but no miRNAs were found to be located in DMRs. Subsequently, the lncRNAs and mRNAs in DMRs were overlapped with their expression level data above to obtain the methylation-related DELs and DEGs. Consequently, five DELs (TMEM51-AS1, HCP5, SND1-IT1, RUSC1-AS1 and LINC00324) were screened (Fig. 3C). Among these DELs, only the expression and methylation levels of lncRNA TMEM51-AS1 (Figs. 3D–3F) were opposite, indicating its expression may be regulated by methylation. These methylation-related genes were used to construct the ceRNA network.
CeRNA network construction
Twenty-four interaction pairs between five DELs and 14 DEMs were predicted using miRcode, starBase and DIANA-LncBase databases (Table 1). The expression trends of these DELs and DEMs were opposite. Subsequently, the target genes of these 14 DEMs were predicted using four algorithms, with the resultant interaction pairs of 700 in TargetScan, 486 in miRBase, 341 in miRanda and 268 in miRTarBase. A total of 404 interaction pairs were ultimately left due to prediction by at least two databases and a negative association between them. These interaction pairs between DELs and DEMs, and between DEMs and DEGs were used to construct a ceRNA network, which contained 258 nodes (five DELs, 11 DEMs and 242 DEGs) (Fig. 4). In this network, TMEM51-AS1 functioned as a ceRNA to regulate SNX21 (sorting nexin family member 21) and TRAPPC10 (trafficking protein particle complex 10) by sponging miR-106b; LINC00324 and RUSC1-AS1 acted as ceRNAs to regulate SPRY4 (sprouty RTK signaling antagonist 4), PAWR (pro-apoptotic WT1 regulator), MICAL2 (microtubule associated monooxygenase, calponin and LIM domain containing 2) and SLC39A14 (solute carrier family 39 member 14) by sponging miR-16; RUSC1-AS1 regulated SCG5 (SCG5 secretogranin V) and PRDM5 (PR/SET domain 5) by competitively binding to miR-10; RUSC1-AS1 also served as ceRNAs for ZFP1 (ZFP1 zinc finger protein) by binding to miR-7; HCP5 could interact with miR-143 to regulate RRM2 (ribonucleotide reductase regulatory subunit M2).
|HCP5||hsa-miR-10, hsa-miR-16, hsa-miR-186, hsa-miR-214, hsa-miR-7, hsa-miR-641, hsa-miR-143, hsa-miR-4770, hsa-miR-216b, hsa-miR-876|
|LINC00324||hsa-miR-143, hsa-miR-16, hsa-miR-214, hsa-miR-216b, hsa-miR-4770|
|RUSC1-AS1||hsa-miR-214, hsa-miR-10, hsa-miR-16, has-miR-216b, hsa-miR-7|
Function enrichment analysis
The DEGs in the ceRNA network was subjected to DAVID to predict their potential functions in LSCC. The results showed that 17 significant GO BP terms were enriched, including GO:0042981 ∼ regulation of apoptosis (PAWR), GO:0015031 ∼ protein transport (SNX21; SCG5), cell cycle (PRDM5) and GO:0043407 ∼ negative regulation of MAP kinase activity 4 (SPRY4). Six KEGG pathways were also enriched, including hsa05210: colorectal cancer, hsa04210:Apoptosis and hsa05205:Proteoglycans in cancer (Table 2).
|Biology process||GO:0006793 ∼ phosphorus metabolic process||0.00122||STK38, SLC20A1, ERBB3, NUAK1, MKNK1, ABI1, PIP5K1A, TRIB1, MAP3K3, SRPK2, MINPP1, ADAM10, STK24, MSH2, PRKCI, PKN2, PTPN12, GAK, MAP4K4, MTMR11, MAPK6, GSK3B, DYRK1A, PTPN1, MAP3K14, ERC1, IKBKB, DUSP7|
|GO:0006796 ∼ phosphate metabolic process||0.00122||STK38, SLC20A1, ERBB3, NUAK1, MKNK1, ABI1, PIP5K1A, TRIB1, MAP3K3, SRPK2, MINPP1, ADAM10, STK24, MSH2, PRKCI, PKN2, PTPN12, GAK, MAP4K4, MTMR11, MAPK6, GSK3B, DYRK1A, PTPN1, MAP3K14, ERC1, IKBKB, DUSP7|
|GO:0006468 ∼ protein amino acid phosphorylation||0.00498||SRPK2, ADAM10, STK38, STK24, NUAK1, ERBB3, PRKCI, PKN2, MKNK1, ABI1, TRIB1, GAK, MAP4K4, MAPK6, MAP3K3, GSK3B, DYRK1A, IKBKB, ERC1, MAP3K14|
|GO:0016310 ∼ phosphorylation||0.00801||SRPK2, ADAM10, STK38, ERBB3, MSH2, STK24, NUAK1, PRKCI, PKN2, MKNK1, ABI1, PIP5K1A, TRIB1, GAK, MAP4K4, MAPK6, MAP3K3, GSK3B, DYRK1A, IKBKB, ERC1, MAP3K14|
|GO:0008104 ∼ protein localization||0.01176||STON2, SEC23A, XPO1, AP1M1, NUP160, PRKCI, CENPF, TMSB10, TRAM2, TAP2, GSK3B, NUP210, TAP1, PIKFYVE, SNX21, RAB23, SCG5, SUPT7L, SAR1B, RAB10, ERC1, KPNA2, KPNB1|
|GO:0043407 ∼ negative regulation of MAP kinase activity||0.00991||STK38, PTPN1, SPRY4, DUSP7|
|GO:0042981 ∼ regulation of apoptosis||0.0165||DLC1, DPF2, IER3, ING3, SYVN1, ERBB3, MSH2, KLF10, PRKCI, AKAP13, CD70, PAWR, SOD2, TNFRSF10A, BAG4, TIAM1, GSK3B, GLO1, APBB2, IKBKB, MYC|
|GO:0043067 ∼ regulation of programmed cell death||0.0181||DLC1, DPF2, IER3, ING3, SYVN1, ERBB3, MSH2, KLF10, PRKCI, AKAP13, CD70, PAWR, SOD2, TNFRSF10A, BAG4, TIAM1, GSK3B, GLO1, APBB2, IKBKB, MYC|
|GO:0015031 ∼ protein transport||0.0188||STON2, SEC23A, XPO1, AP1M1, NUP160, PRKCI, CENPF, TRAM2, TAP2, GSK3B, NUP210, TAP1, SNX21, RAB23, SCG5, SAR1B, RAB10, ERC1, KPNA2, KPNB1|
|GO:0010941 ∼ regulation of cell death||0.0188||DLC1, DPF2, IER3, ING3, SYVN1, ERBB3, MSH2, KLF10, PRKCI, AKAP13, CD70, PAWR, SOD2, TNFRSF10A, BAG4, TIAM1, GSK3B, GLO1, APBB2, IKBKB, MYC|
|GO:0045184 ∼ establishment of protein localization||0.0204||STON2, SEC23A, XPO1, AP1M1, NUP160, PRKCI, CENPF, TRAM2, TAP2, GSK3B, NUP210, TAP1, SNX21, RAB23, SCG5, SAR1B, RAB10, ERC1, KPNA2, KPNB1|
|GO:0008219 ∼ cell death||0.0211||FUS, DPF2, DLC1, IER3, MICB, ERBB3, MSH2, AKAP13, RNF216, PAWR, ITPR1, SOD2, TNFRSF10A, BAG4, UNC5B, TIAM1, SIAH1, MYC, SPAST|
|GO:0010033 ∼ response to organic substance||0.0217||ADAM10, KAT2B, ERBB3, MSH2, KLF10, PRKCI, CALCOCO2, APPL1, TRIB1, B2M, HDAC4, PRKAR2A, SDC1, HDAC2, ADM, TAP2, CTSC, PTPN1, MYC|
|GO:0016265 ∼ death||0.0225||FUS, DPF2, DLC1, IER3, MICB, ERBB3, MSH2, AKAP13, RNF216, PAWR, ITPR1, SOD2, TNFRSF10A, BAG4, UNC5B, TIAM1, SIAH1, MYC, SPAST|
|GO:0044265 ∼ cellular macromolecule catabolic process||0.0228||ADAM10, SYVN1, USP1, RNH1, UBE2V2, RNF216, MYLIP, UBE2Q2, ZFP36L2, FBXW7, GMCL1, PSMD3, ZMPSTE24, SIAH1, PCYOX1, USP33, MYC, FBXO11, USP31|
|GO:0007049 ∼ cell cycle||0.0405||E2F3, KAT2B, MSH2, PAPD7, CENPF, APPL1, GAK, SASS6, MAPK6, GSK3B, PRDM5, PSMD3, ZNF318, HBP1, SIAH1, APBB2, KPNA2, MYC, SPAST|
|GO:0009057 ∼ macromolecule catabolic process||0.0427||ADAM10, SYVN1, USP1, RNH1, UBE2V2, RNF216, MYLIP, UBE2Q2, ZFP36L2, FBXW7, GMCL1, PSMD3, ZMPSTE24, SIAH1, PCYOX1, USP33, MYC, FBXO11, USP31|
|KEGG pathway||hsa05210:Colorectal cancer||0.0377||MSH2, GSK3B, APPL1, MYC, FZD7|
|hsa04210:Apoptosis||0.0421||TNFRSF10A, PRKAR2A, EXOG, IKBKB, MAP3K14|
|hsa00562:Inositol phosphate metabolism||0.0477||MINPP1, TPI1, PIKFYVE, PIP5K1A|
|hsa05169:Epstein-Barr virus infection||0.00175||POLR3F, XPO1, HDAC4, GTF2E2, HDAC2, GSK3B, PSMD3, MAP3K14, IKBKB, MYC|
|hsa05166:HTLV-I infection||0.0122||WNT5A, XPO1, E2F3, KAT2B, MAP3K3, GSK3B, MAP3K14, IKBKB, MYC, FZD7|
|hsa05205:Proteoglycans in cancer||0.0267||WNT5A, EIF4B, SDC1, TIAM1, ERBB3, MYC, FZD7, ITPR1|
Clinical associations of lncRNAs, miRNAs and mRNAs in the ceRNA network
In the 138 miRNA-mRNA matched samples of TCGA data, 114 had OS and 82 had RFS information. Univariate Cox regression analysis in these 138 samples showed that 32 RNAs were significantly associated with OS, including 1 DEL (RUSC1-AS1), 2 DEMs (hsa-miR-16 and hsa-miR-10) and 29 DEGs (i.e., PAWR, SCG5, SPRY4, MICAL2, SNX21, TRAPPC10 and SLC39A14); while 25 RNAs were associated with RFS, including one DEL (LINC00324), three DEMs (hsa-miR-16, hsa-miR-10 and hsa-miR-7) and 21 DEGs (i.e., PRDM5, SCG5, SPRY4, MICAL2 and ZFP1) (Table 3). The OS and RFS related ceRNA networks were extracted independently as shown in Figs. 5A and 5B.
Subsequently, multivariate Cox regression showed TRAPPC10 and SLC39A14 were independent factors for OS; RRM2 was an independent factor for RFS (Table 4). Although SOD2, SLC44A1 and THEM4 were also screened to be significant, their hazard ratios (HR) were not consistent with the expected according to their expression levels. Combined with the univariate results, we suggested TRAPPC10 and SLC39A14 related ceRNA axes (TMEM51-AS1-miR-106-TRAPPC10; RUSC1-AS1-miR-16-SLC39A14) may be especially important. The Kaplan–Meier curve of these lncRNAs, miRNAs and mRNAs were drawn. As expected, the low expression of miR-16 (Fig. 5D) was associated with poor prognosis and the high expression of RUSC1-AS1 (Fig. 5C), SLC39A14 (Fig. 5E) and TRAPPC10 (Fig. 6A) was associated with shorter OS.
|Overall survival||Recurrence free survival|
|ID||P-value||HR||95% CI||ID||P-value||HR||95% CI|
|Lower limit||Upper limit||Lower limit||Upper limit|
P-value < 0.05 shown in bold.
Furthermore, OS- and RFS-related DELs, DEMs and DEGs were also analyzed to investigate their associations with other clinical characteristics of LSCC to further confirm their importance. The results showed that RUSC1-AS1 was significantly associated with Pathologic N; OS- and RFS-related SPRY4 was associated with Pathologic M; OS-related MICAL2 was associated with Pathologic N and Pathologic stage; RFS-related ZFP1 and SLC39A14 were associated with Pathologic N; OS-related SNX21 and RFS-related SCG5 were associated with gender (Table 5). These findings implied SPRY4, MICAL2, ZFP1, SNX21 and SCG5 related ceRNAs (LINC00324/RUSC1-AS1-miR-16-SPRY4/MICAL2, RUSC1-AS1-miR-7-ZFP1, TMEM51-AS1-miR-106-SNX21, RUSC1-AS1-miR-10-SCG5) were also crucial for LSCC. The Kaplan–Meier curve of SNX21 is shown in Fig. 6B and the other DEGs are displayed in Figs. S1 and S2.
Although epigenetics modification has been shown to trigger silencing or overexpression of lncRNAs in cancer (Dong et al., 2017; Guo et al., 2018b; Zhou et al., 2018), the aberrant methylation-mediated expression changes of lncRNAs remain unclear in LSCC. We, for the first time, found that the downregulation of lncRNA TMEM51-AS1 may be mediated by hyper-methylation. Few studies investigated the roles of TMEM51-AS1 in cancer except one study indicated downregulated TMEM51-AS1 was significantly correlated with poor OS in chromophobe renal cell carcinoma (He et al., 2016). In the present study, we predicted that TMEM51-AS1 might function as a ceRNA to regulate SNX21 and TRAPPC10 through sponging miR-106b. Evidence demonstrated that miR-106b was up-regulated in LSCC (Lu et al., 2014; Xing et al., 2014), which was also confirmed in our microarray study. miR-106b was reported to promote the proliferation and invasion of LSCC cells by targeting RUNX3 (Ying et al., 2013), while induce cell cycle G0/G1 arrest by inhibiting tumor suppressor RB (Cai, Wang & Bao, 2011). Although no study revealed the roles of SNX21 in cancer, its family genes, such as SNX1 (Zhan et al., 2018), SNX5 (Jitsukawa et al., 2017) and SNX9 (Bendris et al., 2016) were suggested to be tumor suppressor related. Therefore, SNX21 may be theoretically downregulated in LSCC by miR-106b. Consistent with this hypothesis, our study showed that SNX21 was less expressed in LSCC tissues and patients with high expression of SNX21 had a higher OS rate. There was only one study to suggest the roles of TRAPPC10 until now and showed TRAPPC10 was an oncogenic driver to predict the poor prognosis for breast cancer patients (Pongor et al., 2015), which seemed to be contrast with our results, implying TRAPPC10 may be a new tumor suppressor gene for LSCC. The tumor inhibition effects of TRAPPC10 may be related to its potential to activate GTPase RAB11 (Milev et al., 2018) and the Rab coupling protein, the targeted deletion of which led to accelerated tumor onset (Boulay et al., 2016).
|Clinical characteristics||Significant related|
|Age( ≥60/<60 y)||–||–||ADAM10, FUS, MICAL2, LRRC8E|
|Gender(Male/Female)||–||–||DLC1, HPN, PTPN12, SNX21, ST3GAL2, LRRC8E, NXPH4, PTPN1, SCG5,|
|Alcohol use(Yes/No)||–||–||CNPY3, PYGO2, SLC39A14, TAP2|
|Pathologic_M(M0/-)||–||–||ADAM10, CNPY3, E2F3, SPRY4, LRRC8E, MAP4K4, PTBP1, RRM2|
|Pathologic_N(N0/N1/N2/N3/-)||RUSC1-AS1||–||DLC1, FUS, MICAL2, SOD2, PCYOX1, RRM2, SLC39A14, ZFP1|
|Pathologic_T(T1/T2/T3/T4/-)||–||–||LRRC40, PRPSAP2, SOD2, TSEN15, CELSR2, PTPN1|
|Pathologic_stage(I/II/III/IV/-)||–||–||MICAL2, PRPSAP2, SPRY4, CELSR2, MAP4K4, PCYOX1, PTPN1|
|Grade(G1/G2/G3/G4)||–||-||AHCYL2, DLC1 , PCYOX1, PTPN1, RRM2|
|Tobacco use(Reform/Current/Never)||–||–||HPN, RAB10, SAR1B, SOD2, TAP2|
Underlined genes were recurrence free survival related; the other genes were overall survival related. Bolded genes was both recurrence free and overall survival related.
Furthermore, we identified several other ceRNA axes, although they were not methylation-related, including LINC00324/RUSC1-AS1-miR-16-SPRY4/MICAL2/ SLC39A14, RUSC1-AS1-miR-10-SCG5 and RUSC1-AS1-miR-7-ZFP1. All these lncRNAs, miRNAs and mRNAs were significantly associated with OS and/or RFS, indicating these ceRNA axes may also be underlying therapeutic targets.
Although related report was rare, RUSC1-AS1 (Jian et al., 2015) and LINC00324 (Militello et al., 2017) had been indicated to be highly expressed in cancer cells, which were similarly confirmed in LSCC samples. Accumulating evidence also has proved the roles of miR-16, miR-7 and miR-10 in various types of cancer. miR-16 could be downregulated in tissue samples and cell lines of lung cancer (Ke et al., 2013) and osteosarcoma (Jiao, Wang & Wang, 2018). Ectopic expression of miR-16 inhibited cell proliferation, colony formation in vivo and, migration and invasion in vitro by regulating its target genes RAB23 and Smad3 (Jiao, Wang & Wang, 2018; Zhang et al., 2018). miR-10a was down-regulated in laryngeal epithelial premalignant lesions with increasing grade of dysplasia (Hu et al., 2015). Overexpression of miR-10a inhibited cell metastasis by regulating epithelial-to-mesenchymal transition (EMT) (Liu et al., 2017b). miR-7-5p was lower expressed in brain-metastatic lesions of breast cancer (Hiroshi et al., 2013) and the use of miR-7-5p mimics suppressed cell proliferation and induced apoptosis (Shi et al., 2015) via modulating the expression of Kruppel like factor 4. In agreement with these studies, we also found that these three miRNAs were less expressed (especially miR-10 and miR-7) in LSCC and negatively associated with OS and/or RFS. Although the downstream target genes of these miRNAs have been reported as above, their functions in LSCC remain poorly understood. We predicted that SPRY4/MICAL2/SLC39A14, SCG5 and ZFP1 may be the potential targets of miR-16, miR-10 and miR-7, respectively in LSCC, which had not been validated previously. Nevertheless, the studies on the molecular mechanisms of these DEGs may indirectly explain their potential interactions. The expression of SPRY4 was upregulated in testicular germ cell tumors (Tian et al., 2018). MICAL2 was a recently identified proto-oncogene, which increased cell proliferation to accelerate tumor growth, and promoted the expression of EMT-related proteins to increase cell metastasis (Mariotti et al., 2016; Wang et al., 2018). Immunohistochemical analysis showed the expression level of SLC39A14 was significantly higher in hepatocellular carcinoma tissues than that in adjacent tissues and negatively correlated with survival time (Gartmann et al., 2018). Also, the upregulation of SLC39A14 in tumor cells may be attributed to the loss of its interactive gene p53, a tumor suppressor (Zhao et al., 2017). Although there were no studies to discuss the roles of SCG5 in cancer, its family members secretogranin II and III have been seen to be overexpressed in prostate cancer (Courel et al., 2014) and small cell lung carcinoma (Togayachi et al., 2017), suggesting SCG5 may also be oncogenic for LSCC. Zinc finger proteins had also been observed to promote cell growth and metastasis in nasopharyngeal carcinoma (Li et al., 2015). In line with these findings, SPRY4, MICAL2, SLC39A14, SCG5 and ZFP1 were all upregulated in LSCC and associated with poor prognosis.
There were some limitations in this study. First, although all the known microarray or sequencing data from the public database had been included, the sample size was still not large which may influence the results. Therefore, additional clinical trials with larger samples may be essential to confirm their expression and prognosis. Second, we only preliminarily predicted that these ceRNA axes may be associated with LSCC development and prognosis. The regulatory relationships between lncRNAs and miRNAs as well as between miRNAs and mRNAs needed further experimental confirmation in vitro and in vivo (i.e., dual luciferase reporter assay or loss-of-function). Third, whether the expression of TMEM51-AS1 was regulated by methylation should be validated by using the methylation inhibitor 5-azacytidine. Fourth, although we have normalized the data from different platforms, this may still cause some underlying bias.
Our present study identifies several important mechanisms for the development and progression of LSCC: (1) methylation-mediated upregulation of lncRNA TMEM51-AS1 may function as a ceRNA for miR-106b to regulate SNX21 and TRAPPC10; (2) survival-related RUSC1-AS1/LINC00324 may function as a ceRNA to sponge miR-16, miR-10 or miR-7 and then regulate SPRY4/ MICAL2/SLC39A14, SCG5/PRDM5 and ZFP1, respectively. Altogether, these lncRNA, miRNAs or mRNAs may be potential prognostic biomarkers and therapeutic targets of LSCC.