Introduction
Vulvar squamous cell carcinoma (VSCC) is a relatively rare malignant tumor in women which accounts for 80–90% of female genital tract cancers (1). Early stage VSCC is primarily treated with surgery, however, postoperative complications can occur. Several of these complications, such as lower limb lymphatic obstruction and vulvar morphologic changes, can have a deep impact on the quality of life of patients. The treatment for advanced vulvar cancer includes radiotherapy. The survival rate for patients with advanced VSCC is only ~30% due to the reduced tolerance of vulvar tissues to radiation. Thus, it is important to better understand the pathogenesis of VSCC and to identify possible effective biomarkers to aid diagnosis.
Over the past decade, research into the non-coding RNA (ncRNA) domain has attracted significant attention. Long non-coding RNAs (lncRNAs) >200 nucleotides in length do not encode any proteins (2,3) but are involved in regulating gene functions at transcriptional and post-transcriptional levels (4). These gene functions include epigenetic regulation, chromatin modification, cell cycle regulation, nuclear trafficking, transcription and splicing (5–8). The dysregulation of lncRNAs is linked to various human diseases, particularly cancer (9–11).
In the present study, an lncRNA expression profile was established from samples obtained from four paired VSCC subjects and their adjacent non-tumor (NT) tissues through a microarray platform. The aim of the present study was to identify dysregulated lncRNAs and mRNAs in VSCC patients. Six of these lncRNAs were evaluated by real-time reverse transcription-polymerase chain reaction (RT-PCR) in 35 cases with benign vulvar diseases, VSCC and their pericarcinoma tissues according to the top 30 upregulated and top 30 down-regulated lncRNAs.
Materials and methods
Patients and sample collection
The present study protocol was approved by the Institutional Review Board of Shengjing Hospital Affiliated to China Medical University. Written informed consent was obtained from all patients. A total of 16 VSCC patients, diagnosed with VSCC between January 2013 and June 2015, were enrolled in the present study. Adjacent tissues in 4 VSCC patients, 7 cases of vulvar lichen sclerosus, 6 cases of vulvar leukoplakia and 2 cases of vulvar intraepithelial neoplasia (VIN) grade I were obtained as the control group. The clinical characteristics of the patients and their tumors are summarized in Table I. None of the patients had received prior radiotherapy or other anticancer treatment. All histopathological samples were obtained at surgery and were confirmed as VSCC. All samples were stored at −80°C until use.
Total RNA extraction
Total RNA was isolated from the frozen samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), according to the manufacturer's protocols. All samples were quantified using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA. USA). The purity of total RNA was examined by its absorbance ratio at 260–280 nm. The absorbance ratios at A260/280 were between 1.8 and 2.0.
lncRNA and mRNA microarray analysis
Four paired samples were chosen for the experiments performed at CapitalBio Corporation Laboratories (CapitalBio Corp., Beijing, China). Generally, total RNA which was extracted using TRIzol reagent for microarray was reverse-transcripted to double-stranded complementary DNA (cDNA). cDNA synthesis, purification, labeling and hybridization were carried out according to the manufacturer's instructions. lncRNA expression profiling was performed using lncRNAs and mRNA 4×180K Human Gene Expression Microarray V4.0 (CapitalBio Corp.). The data extracted were analyzed by Agilent Feature Extraction version 10.7 (Agilent Technologies, Santa Clara, CA, USA), and were summarized and normalized using Agilent GeneSpring software version 11.5. Differentially expressed lncRNAs and mRNAs were identified through an absolute fold-change >2 at a p<0.05.
Functional analysis
Gene ontology (GO) and pathway analyses were performed in order to better understand the functions of differentially expressed lncRNAs and mRNAs in VSCC. GO terms consisted of biological processes, cellular components, and molecular functions that were used to annotate and classify gene function. Pathway analysis was based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database that place differentially expressed mRNAs. In brief, the Fisher's exact test and t-test were used. A p<0.05 was considered to indicate a statistically significant result.
lncRNA-mRNA co-expression network analysis
A correlation analysis was constructed between differentially expressed lncRNAs and mRNA by calculating the Pearson correlation coefficients for each dysregulated lncRNA and mRNA. A significant correlation was defined as a correlation >0.99 or ≤0.99 at a p<0.05.
Real-time PCR
Single-stranded cDNA was reverse transcribed using the PrimeScript RT reagent kit with gDNA Eraser Perfect Real-Time (Takara Bio, Kyoto, Japan), according to the manufacturer's instructions. Real-time PCR was performed using SYBR Premix Ex Taq Tli RNase H Plus (Takara Bio) and a Roche LightCycler 480 II system (Roche, Basel, Switzerland). The 20 μl real-time PCR reaction mixture contained 3 μl cDNA, 10 μl SYBR Premix Ex Taq, 4 μl RNase-Free Water, and 1.5 μl, each, of forward and reverse primers. The PCR conditions were as follows: 30 sec at 95°C, followed by 40 cycles at 95°C for 5 sec and at 60°C for 30 sec. All samples were run in triplicate for analysis. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an endogenous control. The relative abundance of lncRNA expression was calculated using the threshold cycle (Ct) method with relative quantitation (12). The primer sequences were designed in the laboratory and were synthesized by Sangon Biotech (Sangon, Shanghai, China), and are shown in Table II.
Statistical analysis
Statistical analysis was performed using Social Sciences (SPSS) 19.0 software (SPSS, Inc., Chicago, IL, USA). The independent-sample t-test was used to evaluate the expression levels of lncRNAs. Association between lncRNAs and clinicopathological factors was analyzed with the Mann-Whitney U test which was used for comparison between the two groups and the Kruskal-Wallis test which was used for comparison among the three groups. A p<0.05 (two-tailed) was regarded as statistically significant.
Results
Expression profiles of lncRNAs in VSCC tissues
From the lncRNA gene expression, dysregulated lncRNAs were identified in four paired vulvar squamous cell tumor tissues and the adjacent normal vulvar tissue samples. Among the 15,840 lncRNA transcripts analyzed, 312 lncRNAs were observed as upregulated (fold-change, >2; p<0.05) whereas 1,469 lncRNAs were detected as downregulated (fold-change >2; p<0.05) in the VSCC samples. Of these, AI769947 (fold-change, 34.2) was the most highly upregulated lncRNA in the VSCC tissues. FER1L4 (fold-change, 48.9) was the most extensively downregulated lncRNA. Hierarchical clustering analysis was performed to determine whether the expression patterns in the VSCC tissues were significantly different from these patterns in the adjacent normal tissues (Fig. 1A).
Overview of the mRNA profile in VSCC
Up to 33,045 coding transcripts were detected in the VSCC tissues and the NTs, of which 21,788 showed significantly differential expression (fold-change, >2; p<0.05). Among them, 1,521 were upregulated and 4,694 were downregulated in the VSCC tissues. Their distinct expression patterns were evaluated by hierarchical clustering analysis (Fig. 1B).
GO and pathway analysis
Go analysis (comprised of biological processes, cellular components and molecular functions) was performed in order to determine the genes and gene product enrichment. Through GO analysis, the significant differentially expressed mRNAs were principally associated with extracellular organelles (Fig. 2), but none were linked with molecular function on GO analysis. Pathway analysis indicated that 38 pathways corresponded to the most significant differential pathways.
Overview of the co-expression network
In order to verify the correlation between differentially expressed lncRNAs and mRNAs, a coding/non-coding gene co-expression network was constructed in VSCC patients compared with the control group. The co-network showed that one lncRNA correlated with multiple mRNAs and reciprocally. According to the network, we found that lncRNA, FLNB-AS1, had the most correlative mRNAs such as HOXD11 and PKM (Fig. 3). Therefore, it was proposed that the expression profile of lncRNAs and mRNAs was significantly correlated.
Real-time PCR validation of several differentially expressed lncRNAs
According to fold difference, one upregulated lncRNA (NEAT1) was initially selected and five downregulated lncRNAs (MIR31HG, MALAT1 HOTAIR, LINC00478 and MEG3) were also selected for verification of expression in two sets of tissue samples. A general consistency between the micro-array and real-time PCR showed that the six selected lncRNAs were dysregulated between VSCC tissues and the control group (Fig. 4). The levels of HOTAIR, MALAT1, MEG3 and LINC00478 (p=0.000, p=0.001, p=0.023 and p=0.009, respectively) expression were decreased and NEAT1 was upregulated as indicated by the microarray analysis (p=0.002); whereas high expression of MIR31HG in the 35 samples, was detected, in opposition to the microarray data (p=0.007).
Evaluation of MEG3 and MALAT1 as diagnostic biomarkers for VSCC
We constructed a receiver operating characteristic (ROC) curve to evaluate the diagnostic value of our findings. Both MEG3 and MALAT1 were downregulated in the VSCC samples. Among them, MEG3 was only detected in 19% (3/16) of the tumor samples, meanwhile, MALAT1 was found downregulated in 87.5% (14/16) of the VSCC samples. Furthermore, the areas under the ROC curve (AUCs) were 0.707 and 0.839, respectively (p<0.05) (Fig. 5).
Correlation between NEAT1, LINC00478 and MIR31HG expression and clinicopathological features in VSCC patients
Lymphatic metastasis is the important measure of prognosis. As shown in Table I, we found that the level of NEAT1 expression in VSCC tissues was associated with lymph node metastasis (p=0.002). We also observed that there was a statistical relationship between LINC00478 and MIR31HG expression and tumor differentiation (p=0.003, p=0.007 respectively). However, no significant correlations were found between lncRNAs levels and other clinicopathological factors studied, including the International Federation of Gynecology and Obstetrics (FIGO) stage.
Discussion
Recent research has demonstrated the key role of lncRNAs in regulating embryogenesis and gene expression, and lncRNAs have emerged as drivers of tumor suppressive and oncogenic functions in various human solid tumors (13–16). Our results revealed that aberrantly expressed lncRNAs may be a factor in VSCC pathogenesis, potentially providing new biomarkers and therapeutic targets for VSCC.
VSCC is a rare malignant tumor, however, its morbidity and mortality rates have risen in the past few decades. It is generally believed that VSCC has two etiological pathways, i.e., a high-risk human papillomavirus (HPV)-dependent route and genetic alterations such as p53 mutations and silencing of p16INK4a (17,18).
Homeobox transcript antisense RNA (HOTAIR) is a well-known lncRNA transcribed from the HOXC cluster located on chromosome 12q13.3 (19). It has been proven that HOTAIR is involved in cancer progression and prognosis and may be an important target for cancer diagnosis and therapy (20–22). According to microarray analysis, the present study revealed that the level of HOTAIR expression was decreased in 75% (3/4) of the paired samples of tumor tissues. According to Sharma et al, there is a significant linear trend towards progressive HOTAIR downregulation among HPV-negative controls, HPV-positive non-malignant tumors and cervical carcinoma samples. They speculated that HOTAIR could be a potential target of E7, in HPV16-related cervical cancers (23).
Metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) is another long lncRNA correlated with HPV. Jiang et al found that MALAT1 was expressed in HPV-positive cervical squamous cells, but not in HPV-negative normal cervical squamous cells, which suggests that HPV correlates with MALAT1 deregulation in cervical cancer (24). Various previous studies have indicated that HOTAIR and MALAT1 are overexpressed in tumor tissues compared with the normal tissues (25–28). It was thought that HOTAIR and MALAT1 contribute to carcinogenesis in most tumors whereas other studies proved the role of HOTAIR and MALAT1 as tumor suppressors. HOTAIR expression was found to be increased in breast cancer patients; nonetheless, knockdown of HOTAIR could inhibit invasion (15). Lu et al found that breast cancer patients with HOTAIR high-expression had lower risk of recurrence and poor prognosis (29). It has also been demonstrated that the level of expression of MALAT1 in gliomas was lower than that in normal tissues (30). These discrepant results suggest that lncRNAs can differentially express depending on the cancer types and that lncRNAs have active roles in tumorigenesis. Moreover, our data revealed the opposite results when compared to many other studies, namely that HOTAIR and MALAT1 were down-regulated in VSCC tissues. Firstly, we inferred that HPV infection may be involved. Both HOTAIR and MALAT1 are correlated with HPV. Unfortunately, some of the patients we evaluated refused HPV examination and we were unable to determine the number of patients originally infected with HPV. HPV infection is not the most crucial factor in the level of expression of lncRNAs, since it is irrelevant in other types of tumors. However, HPV examination is still necessary for VSCC patients. In addition, it is meaningful that HPV infection be considered in lncRNA expression analysis as it is a useful factor. We hypothesized that lncRNAs such as HOTAIR and MALAT1 are expressed differentially in HPV-positive and HPV-negative VSCC samples and that these two lncRNAs may be involved in VSCC tumorigenesis by targeting HPV. Secondly, the number of cases is one limitation in our research. Thus, larger-scale studies are needed to validate our findings for further research.
The TP53-associated route and other genetic changes, such as CDKN2A (p16), also play a key role in the development of VSCC. Genetic mutations in TP53 have been detected in HPV-negative VSCCs. They may represent early changes in HPV-independent vulvar carcinogenesis (31). Expression of p53 was significantly increased in VSCC patient samples compared with the control group (32).
Maternally expressed gene 3 (MEG3) is a lncRNA which has been recently associated with p53. Lu et al found that p53 protein levels are affected by MEG3 overexpression in vitro. MEG3 was significantly downregulated in non-small cell lung cancer (NSCLC) tissues that may be affected by DNA methylation, and it partially regulated NSCLC cell proliferation and apoptosis via the activation of p53 (33). Zhu et al found that the level of MEG3 expression was reduced in hepatoma samples compared with adjacent non-tumor samples and interacts with the p53 DNA binding domain. Moreover, various p53 target genes were found to be deregulated after overexpression of MEG3 in hepatoma cells (34). Our data showed that MEG3 was downregulated in the VSCC patient samples. We speculate that MEG3 may operate as a tumor-suppressor in VSCC, whose mechanism may involve activation of p53.
Epigenetic silencing of p16 is also an early and important event in vulvar neoplasia (19). Montes et al showed that the lncRNA, MIR31HG, is upregulated in oncogene-induced senescence (OIS) which is considered an important mechanism of tumor suppression. They also found that knockdown of MIR31HG promoted a strong p16INK4A-dependent senescence phenotype. These results suggest that MIR31HG may repress p16INK4A expression in human cancer (35). MIR31HG is underactive in most cancers, but in lung and breast cancers, it appears to be upregulated by an, as yet, unknown mechanism (36,37). According to the microarray analysis, the present study revealed that the level of MIR31HG expression was decreased in 75% (3/4) of the paired samples of tumor tissues while it was upregulated in VSCC samples by RT-PCR. These findings imply that the microarray profiling results that are based on a small sample size may not entirely be reliable and need to be further verified by RT-PCR analysis using a greater number of samples.
The lncRNA MIR31HG is also correlated with poor prognosis and may be an important target for cancer diagnosis and therapy in gastric cancer (38). According to the present study, we found that the increased level of MIR31HG expression in VSCC was significantly associated with tumor differentitation (Table I), suggesting that MIR31HG may play a critical role, not only in the tumorigenesis of VSCC, but also in its prognosis. Patients with higher expression levels of MIR31HG may have poorer prognosis. However, not all patients who were diagnosed with VSCC within the last three years at our institution underwent surgery, and only two patients died of VSCC, therefore, we were unable to obtain reliable survival analysis.
Nuclear paraspeckle assembly transcript 1 (NEAT1) is a prognostic biomarker as well as MIR31HG (39). NEAT1 was upregulated in 75% (3/4) of VSCC tissues in accordance with the results of the microarray profile (Fig. 4). According to the study of Li et al, NEAT1 may play an oncogenic role in colorectal cancer which is involved in differentiation, invasion and metastasis (40). Our results revealed a greater association between NEAT1 expression and lymph node metastasis (Table I). Lymphatic metastasis usually indicates poor prognosis, suggesting that NEAT1 may play an important role in the prognosis of VSCC.
lncRNAs are emerging as biomarkers based upon ncRNA biology. Similar to microRNAs, various lncRNAs have been found in both urine and plasma, which can be easily obtained and provide for readily-available and stable diagnostics in order to detect cancers and cancer subtypes (41,42).
lncRNAs, such as MEG3 and MALAT1, have been shown to be potential biomarkers in other types of cancers. However, there are some critical issues that still need to be resolved before they can be applied clinically. Most recent studies have been performed using only small-scale samples. Therefore, further large-scale studies are needed to validate their findings. In addition, since conventional methods such as RNA extraction take too long and the quality of lncRNAs can be influenced by multiple factors during the analytical process, more stable platforms and more rapid analytical methods are required. Finally, the expression of lncRNAs could be confounded by multiple factors in tumorigenesis procedures such as pathology. Thus, our results must be validated in larger independent cohorts to confirm that lncRNAs can be used as reliable biomarkers for cancer detection.
The present study had several limitations including its small sample size. Although our results suggested that several lncRNAs were potentially important biomarkers in VSCC, we collected only 16 samples from VSCC patients. Thus, our results require validation using larger prospective studies involving a larger number of samples. Furthermore, not all patients who were diagnosed with VSCC within the last three years underwent surgery; some were lost in the follow-up and only two patients died from VSCC. Thus, it was difficult to arrive at a reliable survival curve due to the small number of cases available. In addition, there were no stage IV patients in our cases since most patients with stage IV chose not to undergo surgery. Therefore, future experiments targeting these aspects are needed to verify our results.
In conclusion, the present study was the first to determine global aberrant lncRNA expression in VSCC tissues compared with the adjacent non-tumorous tissues using microarray analysis. We also identified a panel of dysregulated lncRNAs that may be potential biomarkers as they were also correlated with VSCC carcinogenesis. These dysregulated lncRNAs may be VSCC-specific lncRNAs and may form the basis for future diagnostic, therapeutic and functional research on VSCC.