In this study, we collected three pairs of invasive cervical cancer tissues (S102_TN_5_NS, S369_TN_3_NS, S372_TN_2_NS) and paracancerous tissues (Ctl_1_NS, Ctl_2_NS, and Ctl_3_NS) from the Obstetrics and Gynecology Hospital of Fudan University (Shanghai, China). These patients were histologically diagnosed with squamous cervical carcinoma (SCC) or adenocarcinoma and classified according to the National Comprehensive Cancer Network (NCCN) criteria [21]. Immunohistochemical results showed that these tumor tissues were CK7(+), P16(+), and Ki-67(+). The clinical characteristics of these cases are shown in Table 1. Tissue specimens were collected from the surgery room and snap-frozen for RNA isolation (see next subsection). This study was approved by the Ethics Committee of the Obstetrics and Gynecology Hospital of Fudan University (2018-11).

Total RNA was isolated using a RNAiso Plus kit (cat. #9109; Takara Bio, Kusatsu, Shiga Prefecture, Japan) according to the manufacturer’s instructions, quantified using a spectrophotometer, and then analyzed with an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA) for RNA integrity. The resulting RNA samples were further purified with the NucleoSpin RNA Clean-up XS kit (Macherey Nagel, Strabe, Germany) and RNase-free DNase set (Qiagen, GmBH, Germany) according to the manufacturers’ protocols.

Six human ceRNA arrays (4 $\times$ 180 K) were purchased from Shanghai Biotechnology Corporation (SBC; China) and used for the microarray analysis of these tissue samples. In brief, the quantified and quality-controlled RNA samples (see above) were labeled using the Low Input Quick Amp Labeling Kit, One-Color (Agilent Technologies) as probes following the kit’s instructions. The labeled cRNA samples were purified using a Qiagen RNease Mini kit and then hybridized to the ceRNA microarray slides in a hybridization oven (Agilent Technologies) overnight (1.65 $\mu$g Cy3-labeled cRNA for each reaction) using the Gene Expression Hybridization Kit (Agilent Technologies), according to the manufacturer’s instructions. After that, the microarrays were intensively washed with the washing buffer from the Wash Buffer Kit (Agilent Technologies) following the kit’s protocol. In the end, the microarrays were scanned using the Agilent Microarray Scanner to collect data. For data analysis, Agilent Technologies (Shanghai Biotechnology Corporation) helped us to extract and analyze the scanned microarray data using the Feature Extraction software 12.0 (Agilent Technologies, https://www.agilent.com.cn/zh-cn/product/mirna-microarray-platform/mirna-microarray-software/feature-extraction-software-228496) and limma software packages in the R program. The raw data were normalized using the Quantile algorithm with the limma packages. Changes in gene expression level with a fold-change of $\ge$2 or $\le$0.5 were used as a cut-off value for the differentially expressed lncRNAs or mRNAs (upregulated vs. downregulated) in cervical cancer vs. normal tissues.

The microarray data of differentially expressed lncRNAs and mRNAs in cervical cancer tissues vs. normal ones were first analyzed using the GO database (http://www.geneontology.org/) for their functional categories as the specific biological GO terms [22]. Then, we performed the pathway analyses of the differentially expressed lncRNAs and mRNAs using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, assisted by the Shanghai Biotechnology Corporation, to examine their potentially related gene pathways in cervical cancer tissues vs. normal ones.

We further predicted and evaluated the targeting mRNAs of the differentially expressed lncRNAs for their cis/trans-regulatory effects on mRNAs. Specifically, cis-targeted mRNAs were considered as targets of lncRNAs when the mRNA genes were localized within 10 kb of the paired lncRNAs, while trans-targeting occurred when lncRNAs contained the specific gene sequences. After selecting the complementary or similar gene sequences paired with lncRNAs using the Basic Local Alignment Search Tool (BLAST) online software, we calculated the RNA duplex energy level between the two sequences using the RNA duplex energy parameter of $\le$–30.

To assess the relationship between lncRNAs and mRNAs, we conducted the lncRNA-mRNA co-expression network analysis of differently expressed lncRNAs (total of 1446, including 242 upregulated and 1204 downregulated) and mRNAs (total of 1300, including 169 upregulated and 1131 downregulated) using Pearson’s correlation coefficients (PCCs). The lncRNA-mRNA co-expression network was constructed by Cytoscape software (Version 3.7.1; The Cytoscape Consortium, San Diego, CA, USA) and analyzed statistically. Correlation with a value of –0.99 $\ge$ x $\ge$ +0.99 between lncRNAs and mRNAs expression levels was considered meaningful. Overall, we obtained a total of 653 nodes and 400 edges to be included in the network analysis.

To verify the microarray data, we selected 15 mRNAs and 10 lncRNAs randomly and assessed their expression levels in three cases of cervical cancer and normal tissue specimens using RT-qPCR. Total RNA was isolated from tissue samples and reversely transcribed into cDNA using an RNAiso Plus kit (Takara) and HiScript II Q RT SuperMix (Vazyme Biotech, Nanjing, China), respectively, according to the kits’ instructions. These cDNA samples were subjected to qPCR amplification using the SYBR Green RT RCP master mix (Genetech, Shanghai, China) with specific primers in a 7900 HT Sequence Detection System (Applied Biosystems, Foster City, CA, USA). The primers of mRNAs and lncRNAs were designed and synthesized by Sinotech Genomics (Shanghai, China) and are listed in Table 1. The relative fold-change of each mRNA and lncRNA was calculated using the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }Ct}$ method after normalized to the level of glyceraldehyde 3-phosphate dehydrogenase (GAPDH) mRNA. The experiment was performed in triplicate and repeated at least once.

The data were statistically analyzed using SPSS (Version 17.0; SPSS Inc., Chicago, IL, USA) and are presented as the mean $\pm$ standard deviation (SD) or fold-change. The statistical analysis included Student’s t-test, where a p 0.05 is considered statistically significant.

To assess the RNA expression profiles in cervical cancer, we first conducted a whole transcriptome microarray analysis in three pairs of cervical tumor and normal tissues using the ceRNA chips (Shanghai Biotechnology, China), which contain 88,371 circRNA probes, 68,423 lncRNA probes, and 18,853 mRNA probes. A total of 1300 DE-mRNA and 1446 DE-lncRNA transcripts were identified in cervical cancer development with a threshold of fold-change $\ge$2 or $\le$–2; p $\le$ 0.05. Particularly, we identified 242 upregulated and 1204 downregulated lncRNAs and 169 upregulated and 1131 downregulated mRNAs. Among these genes, NR_001298 was the most upregulated and lnc-DDT-2:2 the most downregulated lncRNAs, while NM_021983 (HLA-DRB4) was the most upregulated and NM_052892 (PKD1L2) the most downregulated mRNAs. The heat map data and volcano plot data are shown in Fig. 1A–D, respectively.

Fig. 1.

Identification of differentially expressed mRNAs and lncRNA in cervical cancer. (A) Heat maps of differentially expressed mRNAs. (B) Heat maps of differentially expressed lncRNAs, which are based on the expression values of all lncRNAs and mRNAs detected through the microarray analysis. (C) Volcano plots of differentially expressed mRNAs. (D) Volcano plots of differentially expressed lncRNAs, and cutoff criterion was fold-change of $\ge$2 or $\le$–2; p $\le$ 0.05.

To explore the possible functions of lncRNA, we predicted its potential targets in cis-regulatory relationships by searching for protein-coding genes 10-kb upstream and downstream of all the identified lncRNAs. To gain further insight into the potential biological and pathophysiological processes that may be moderated by these DE-RNAs in cervical cancer, we performed the GO term and KEGG pathway analyses. Concerning the biological process, we found that DE-lncRNA was mostly enriched in extrinsic apoptotic-signaling pathway via death domain receptors and was involved in nuclear-transcribed mRNA poly(A) tail shortening and metalloaminopeptidase activity. The top 30 GO terms are shown in Fig. 2A. These results suggest that one of the primary roles of lncRNA may be post-transcriptional regulation of gene expression. Undistinguishably, the GO terms of DE-mRNA were mainly enriched in protein K11-linked ubiquitination, mitochondrial fusion, and the anaphase-promoting complex. The KEGG pathway analysis indicates that these up- and downregulated mRNAs function in ubiquitin-mediated proteolysis and molecular pathways involved in cancer, respectively (Fig. 2A). In particular, the lncRNAs and cis-targeted mRNAs showed to play in a role in the TGF-$\beta$ signaling, cell adhesion molecules (CAMs), and cancer pathways (Fig. 2B).

Fig. 2.

Data obtained from the Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of differentially expressed mRNAs and lncRNAs in cervical cancer. (A) GO term annotations and KEGG pathway enrichment analysis of mRNAs with top 30 enrichment scores. (B) KEGG pathway enrichment analysis of cis-targeted lncRNAs and trans-targeted lncRNAs with the top 30 enrichment scores.

The functions of lncRNA primarily depend on interaction with mRNAs, which may be an important target of the former. It has been reported that the functions of most lncRNAs include the targeting of transcription, translation, and post-translation of their targeting mRNAs. To gain a better understanding of these DE-RNAs in cervical cancer, we determined the potential targets of lncRNAs in cis-regulatory relationships and performed a co-expression analysis using Pearson’s correlation analysis with Psych packages (Corr. Test method, parameter CI = F, adjust = BH) in the R program. The screening correlation coefficient was determined to be r $\ge$ 0.7 with p 0.05, then Cytoscape (https://cytoscape.org/, version 3.7.0) was used to construct the co-expression network. In total, we were able to identify 653 nodes and 400 edges for the network (Fig. 3). It is suggested that an mRNA may be associated with one to six lncRNAs, while one lncRNA is correlated to one to three mRNA. Significantly, six mRNAs were shown to most connect with this network, namely NM_033318 (SMDT1), NM_021090 (MTMR3), NM_021964 (ZNF148), NM_003026 (SH3GL2), NM_018127 (ELAC2), and NM_017934 (PHIP) (Table 2). A functional analysis was conducted to uncover the related biological functions of these mRNAs, revealing that the co-expressed genes were enriched in 51 GO terms (31 under the biological process, 11 GO under the cellular component, and 9 under the molecular function). Particularly, some of the terms were tumor-associated terms, including autophagy (GO:0006914), regulation of vascular genesis (GO:2001214), transcriptional repressor complex (GO:0017053), and regulation of gene expression (GO:0010629) (Fig. 4A). Furthermore, the co-expressed genes were found to be enriched in Herpes simplex virus 1 infection (Fig. 4B). These findings indicate that lncRNAs also regulate trans-target genes.

Fig. 3.

Data on the differentially expressed lncRNA-mRNA co-expression network. The lncRNA-mRNA co-expression network was constructed using the Cytoscape software (https://cytoscape.org/, version 3.7.0) and analyzed statistically. Correlation values of –0.7 $\ge$ x $\ge$ +0.7.

Fig. 4.

Data on GO terms and KEGG analysis of genes in the differentially expressed lncRNA-mRNA co-expression network. (A) GO analysis genes in co-expression network. (B) KEGG analysis of genes in co-expression network.

To assess the prognostic value of these above hub genes, we searched and downloaded data from the online database GEPIA2 (http://gepia2.cancer-pku.cn/#index), which contains gene expression and overall survival data. Importantly, GEPIA2 is an enhanced web server for large-scale gene expression profiling and interactive analysis that includes 198,619 isoforms and 84 cancer subtypes from TCGA and GTEx samples [23]. Our results demonstrate that MTMR3 has a lower expression in cervical squamous cell carcinoma and endocervical adenocarcinoma. Based on the Kaplan-Meier data analysis data, patients with higher MTMR3 expressed tumor had poor overall survival (OS), although the p-value was not statistically significant (HR = 0.94, p-value = 0.78; Fig. 5A). Similarly, the low expression of SMDT1 was also correlated with poor overall survival (OS) (HR = 0.65, p-value = 0.077; Fig. 5B). As for the other hub genes, association of their expression with survival of patients was not statistically significant (Supplementary Fig. 1).

Fig. 5.

Different expressions of the hub genes in cervical cancer vs. normal tissues and association with overall survival. (A and C) Expression levels of MTMR3 or SMDT1. (B and D) Overall survival analysis using Kaplan-Meier curves and the log-rank test.

To verify the reliability of our microarray analysis results, we performed RT-qPCR validation in 15 randomly selected mRNAs and 10 randomly selected lncRNAs from cervical cancer and normal tissue samples (Primers used were shown in Table 3). We found that two lncRNAs (lnc-CLRN3-3:1 and lnc-TMEM238-2:1) were upregulated in cervical cancer tissues, whereas ENST00000369391, ENST00000547898, ENST00000561507, lnc-SEMA3A-2:1, lnc-STOM-1:1, NR_110418, and NR_026867 were verified to be downregulated in cervical cancer tissues vs. those of the control tissues (Fig. 6A). Moreover, according to the validation of differentially expressed mRNAs, we discovered that CCNE2, FN1, INHBA, and UBE2S were upregulated, while APPL1, GSN, KAT2B, MITF, PIK3R1, PRICKLE1, SAV1, SMAD4, TCF7L2, TGFBR2, and UBE2D3 were downregulated in cervical cancer tissues vs. control tissues (Fig. 6B). The qPCR results confirm that the expression patterns of these RNAs are consistent with their expression levels calculated from the RNA-seq data (The comparative results were shown in Supplementary Fig. 2), which indicates the reliability of high-throughput sequencing.

Fig. 6.

RT-qPCR validation of microarray data in cervical cancer vs. normal control tissue specimens. (A) Relative expression levels of lncRNAs in cervical cancer tissues. (B) Relative expression levels of mRNAs in cervical cancer tissues.