lncRNA TMEM51-AS1 and RUSC1-AS1 function as ceRNAs for induction of laryngeal squamous cell carcinoma and prediction of prognosis

Background Long non-coding RNAs (lncRNAs) can function as competing endogenous RNAs (ceRNAs) to interact with miRNAs to regulate target genes and promote cancer initiation and progression. The expression of lncRNAs and miRNAs can be epigenetically regulated. The goal of this study was to construct an lncRNA-miRNA-mRNA ceRNA network in laryngeal squamous cell carcinoma (LSCC) and reveal their methylation patterns, which was not investigated previously. Methods Microarray datasets available from the Gene Expression Omnibus database were used to identify differentially expressed lncRNAs (DELs), miRNAs (DEMs), and genes (DEGs) between LSCC and controls, which were then overlapped with differentially methylated regions (DMRs). The ceRNA network was established by screening the interaction relationships between miRNAs and lncRNAs/mRNAs by corresponding databases. TCGA database was used to identify prognostic biomarkers. Results Five DELs (downregulated: TMEM51-AS1, SND1-IT1; upregulated: HCP5, RUSC1-AS1, LINC00324) and no DEMs were overlapped with the DMRs, but only a negative relationship occurred in the expression and methylation level of TMEM51-AS1. Five DELs could interact with 11 DEMs to regulate 242 DEGs, which was used to construct the ceRNA network, including TMEM51-AS1-miR-106b-SNX21/ TRAPPC10, LINC00324/RUSC1-AS1-miR-16-SPRY4/MICAL2/ SLC39A14, RUSC1-AS1-miR-10-SCG5 and RUSC1-AS1-miR-7-ZFP1 ceRNAs axes. Univariate Cox regression analysis showed RUSC1-AS1 and SNX21 were associated with overall survival (OS); LINC00324, miR-7 and ZFP1 correlated with recurrence-free survival (RFS); miR-16, miR-10, SCG5, SPRY4, MICAL2 and SLC39A14 were both OS and RFS-related. Furthermore, TRAPPC10 and SLC39A14 were identified as independent OS prognostic factors by multivariate Cox regression analysis. Conclusion DNA methylation-mediated TMEM51-AS1 and non-methylation-mediated RUSC1-AS1 may function as ceRNAs for induction of LSCC. They and their ceRNA axis genes (particularly TMEM51-AS1-miR-106b-TRAPPC10; RUSC1-AS1-miR-16-SLC39A14) may be potentially important prognostic biomarkers for LSCC.


INTRODUCTION
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 (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, hypermethylation 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.
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://tcgadata.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.

Data preprocessing
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 (tau 2 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., 2001Shannon et al., -2008 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 prognosisrelated 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 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).

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).

DISCUSSION
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,

Notes.
Underlined genes were recurrence free survival related; the other genes were overall survival related. Bolded genes was both recurrence free and overall survival related.
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).
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 epithelialto-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.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This work was supported by Natural Science Foundation of Liaoning Province of PR China (No. 20170541050). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosures
The following grant information was disclosed by the authors: Natural Science Foundation of Liaoning Province of PR China: 20170541050.