Identification of potential long noncoding RNA associated with nasopharyngeal carcinoma using deep sequencing.

Introduction

Nasopharyngeal carcinoma (NPC) is an epithelial malignancy characterized by a unique set of geographical, aetiological and biological features distinct from those of other head and neck cancers.[1] The incidence is high in the south-eastern region of China but rare in the western world.[2] Many patients with NPC are missed by the current diagnostic approaches in the early stages, leading to treatment failure.[3] In addition, the high incidence of distant metastasis in NPC is an important risk factor of mortality.[4] To date, radiotherapy and chemotherapy are the major strategies to improve the survival rate of patients.[4] However, the molecular mechanism of NPC remains unknown.

Long noncoding RNAs (lncRNAs) are noncoding RNAs that are usually longer than 200 nucleotides and exert biological functions through their interaction with DNA and proteins.[5] LncRNAs play important roles in transcription control, post-transcriptional processing, chromatin remodelling and protein metabolism.[6] However, little is known about their associations with NPC, although they have been reported as biomarkers specifically expressed in the progression of various tumours.[6]

Deep sequencing can be used to explore a large lncRNA pool and has obvious advantages for the identification of lncRNA sequence variations and the discovery of novel lncRNAs.[7] A previous study identified 2127 differentially expressed lncRNAs (DELs) in patients with bladder urothelial cancer, four of which had potential value in the treatment and diagnosis of tumours.[8] However, knowledge of the role of lncRNAs in NPC remains limited.

The aim of the current study was to identify DELs in NPC tissues compared with chronic nasopharyngitis (CNP) tissues. The study then confirmed a subgroup of DELs using quantitative polymerase chain reaction (qPCR) in order to determine a DEL-based risk score in order to investigate the potential regulatory network of these lncRNAs.

Patients and methods

Patient samples

This prospective cohort study enrolled patients with NPC and patients with CNP in the Department of Otorhinolaryngology, Head and Neck Surgery, Shenzhen Baoan Hospital, Southern Medical University, Shenzhen, China between January 2016 and December 2017. None of the patients had received therapy before their biopsy. Fresh tissues were stored in liquid nitrogen for transportation and at –70°C afterwards. The nature of all specimens was confirmed by pathological examination. The Tumor-Node-Metastasis classification was used according to the criteria of the 7th edition of the International Union for Cancer Control and the American Joint Committee on Cancer.[9]Samples of NPC and CNP were used for deep sequencing.

The study was approved by the Ethics Committee of the Affiliated Shenzhen Baoan Hospital of Southern Medical University (no. 2016010631). All patients providing samples for the study were required to provide written informed consent.

RNA extraction and RT–qPCR

Total RNA, containing lncRNAs, was extracted from freshly frozen tissues using a QIAGEN miRNeasy Mini Kit (Qiagen, Hilden, Germany). Sample quality and integrity were assessed using an Agilent 2100 Bioanalyzer system (Agilent Technologies, Santa Clara, CA, USA). A reverse transcription (RT)–qPCR assay was carried out using a PrimeScript™ RT reagent kit (Takara Bio, Kusatsu, Japan) according to the manufacturer’s instructions. The primers were synthesized by Kangcheng (Shanghai, China) and the sequences are listed in Table 1. The cycling programme involved preliminary denaturation at 95 °C for 30 s, followed by 40 cycles of denaturation at 95 °C for 5 s and annealing at 62 °C for 34 s. Relative lncRNA levels were normalized relative to the expression of glyceraldehyde 3-phosphate dehydrogenase and calculated using the 2−ΔΔCt method. Considering the observed fold changes, P-values and the gene specificity obtained in the RNA sequencing experiments, the following six lncRNAs were selected for validation: ENSG00000250451, ENSG00000243479, ENST00000230798, ENST00000228630 (upregulated), ENSG00000230489, ENSG00000227084 (downregulated). qPCR was performed using a SYBR Green mix reagent kit (Takara Bio) with an ABI 7500 Real Time PCR system (Applied Biosystems, Foster City, CA, USA).

Table 1.

Primer sequences used in reverse transcription–quantitative polymerase chain reaction analysis of long noncoding RNAs.

Transcript identification	Forward primer	Reverse primer
ENSG00000230798	5ʹ-GGTGGAGGAGGCGAGGATG-3ʹ	5ʹ-AGCGGACAGACAGGGATTGG-3ʹ
ENSG00000250451	5ʹ-CTGCGACACTTCCCCACC-3ʹ	5ʹ-CTACTTGCCCACGACCGA-3ʹ
ENSG00000243479	5ʹ-CGAACCTTATCTGCTATGGG-3ʹ	5ʹ-ATCTGGAGGGTAGTTCTGCTG-3ʹ
ENSG00000228630	5ʹ-CAGTGGGGAACTCTGACTCG-3ʹ	5ʹ-GTGCCTGGTGCTCTCTTACC-3ʹ
ENSG00000230489	5ʹ-TTCACCCTCGAGCTGGAATC-3ʹ	5ʹ-ATCCCATCTTGCTGTGACCC-3ʹ
ENSG00000227084	5ʹ-AGCTGGGCAAGCACCTAAAA-3ʹ	5ʹ-TGGGGTTATCTTGCGTATGCT-3ʹ
Glyceraldehyde 3-phosphate dehydrogenase	5ʹ-AGGGCTGCTTTTAACTCTGGT-3ʹ	5ʹ-CCCCACTTGATTTTGGAGGGA-3ʹ

Library preparation and lncRNA RNA sequencing

Isolated total RNA was subjected to ribosomal RNA depletion using a RiboZero Kit (Epicentre, Madison, WI, USA), cut randomly into short fragments and then reverse-transcribed into cDNA (Illumina, San Diego, CA, USA). Sequencing adapters were ligated to both ends of the cDNA fragments. After amplification of the adapter-ligated fragment by PCR, fragments with inserts size between 200 and 400 base pairs (bp) were selected and sequenced using an Illumina HiSeqTM 2500 sequencer (Illumina) according to the standard Illumina protocol. Analysis of the quality and filtering of the raw RNA-sequencing (seq) data to remove low quality reads, adaptor sequences, contaminant DNA and PCR duplicates was performed using the Trimmomatic software version 0.32 (Illumina). The trimmed RNA-seq reads were mapped to the human reference genome hg 19 using TopHat2 software version 2.0.13.[10] Cufflinks software version 1.0.3 was used to process the reads and calculate the transcription level of each gene (University of Washington, Seattle, WA, USA). To obtain the differential expression profiles, the fragment per kilobase of transcript per million mapped reads (FPKM) values were calculated for each lncRNA.

KEGG enrichment analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG) was used to understand the biological significance of the differentially expressed genes. Fisher’s exact test and χ2-test were used as statistical tests and the threshold of significance was defined based on the false discovery rate and P-value. Macrogen Korea (Seoul, Korea) performed the bioinformatic analyses.

Statistical analyses

All statistical analyses were performed using IBM SPSS Statistics for Windows, Version 22.0 (IBM Corp., Armonk, NY, USA) and GraphPad Prism (GraphPad, La Jolla, CA, USA). Heat maps were obtained using the R package.[11] Student’s t-test was used to analyse the differences in lncRNA expression levels between the NPC and CNP groups. The associations between clinicopathological parameters were determined using χ2-test and Fisher’s exact test. A P-value < 0.05 was considered statistically significant. The diagnostic value was evaluated with a receiver operating characteristic (ROC) curve analysis.

Results

A cohort of patients, 30 with NPC (7 females and 23 males; mean ± SD age 38.5 ± 14.5 years) and 27 with CNP (8 females and 19 males; mean ± SD age 41.0 ± 16.0 years), were enrolled in the study.

Four patients with NPC and four with CNP were randomly selected and their RNA was extracted for deep sequencing. The sequencing data have been deposited in the Sequence Read Archive database of the National Center for Biotechnology Information, US National Library of Medicine (BioProject: PRJNA451367; SRA: SRP142570). The Cufflinks method was used to identify significantly deregulated genes in NPC (tumour, T) and CNP (inflammation, I) samples. The gene expression level was normalized as FPKM. A correlation analysis of the gene expression among samples was performed using the Pearson’s coefficient (r) of the Log2 (FPKM+1) value. The gene expression was generally highly correlated between the NPC and CNP groups (lncRNAs Pearson’s correlation = 0.9857, Figure 1a). Using the criteria of fold-change of FPKM (absolute value by Log2 ≥ 1) and P < 0.05, a total of 296 DELs were identified. Of the 296 lncRNAs, 65 were upregulated and 231 were downregulated in NPC (Figure 1a). The top 20 deregulated genes are listed in Table 2.

Figure 1.

Long noncoding RNA (lncRNA) expression profile in tissue samples. (a) LncRNA expression profiles in patients with nasopharyngeal carcinoma (NPC-T in the graph) and chronic nasopharyngitis (CNP-I in the graph) using deep sequencing technology. (b) Heat map showing distinguishable lncRNA expression profiles between patients with NPC (T1–T4) and chronic nasopharyngitis CNP (I1–I4). The colour version of this figure is available at: http://imr.sagepub.com.

Table 2.

The top 20 differentially expressed long noncoding RNAs (lncRNAs) in nasopharyngeal carcinoma compared with chronic nasopharyngitis.

Number	lncRNA	Fold change	Change
1	ENSG00000230798	12.00668	up
2	ENSG00000249069	11.90877	up
3	ENSG00000250451	11.90313	up
4	ENSG00000259508	–11.13089	down
5	ENSG00000226779	11.04815	up
6	ENSG00000243479	10.75864	up
7	ENSG00000203372	–10.67287	down
8	ENSG00000226510	10.50432	up
9	ENSG00000227084	–10.4226	down
10	ENSG00000230838	10.12606	up
11	ENSG00000236858	10.09011	up
12	ENSG00000269954	–10.0119	down
13	ENSG00000253665	9.593391	up
14	ENSG00000230489	–9.54593	down
15	ENSG00000228630	9.476746	up
16	ENSG00000251151	9.083479	up
17	ENSG00000237923	8.965784	up
18	ENSG00000225826	8.813781	up
19	ENSG00000244321	8.550747	up
20	ENSG00000258077	8.501837	up

To visualize the DELs and further explore the difference between the NPC and CNP samples, hierarchical cluster analysis was undertaken. The results of this analysis showed that lncRNAs could be clearly classified into two groups: the tumour and chronic inflammation groups (Figure 1b).

The KEGG pathway analysis showed that the DELs annotated to the Rap1 signalling pathway, natural killer cell-mediated cytotoxicity and the thyroid hormone signalling pathway were statistically significant (Table 3) (P-value < 0.05 for all comparisons). LncRNA ENSG00000227084 and LncRNA ENSG00000230489 were related to the Rap1 signalling pathway and natural killer cell-mediated cytotoxicity, respectively.

Table 3.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway of enriched differentially expressed target genes.

KEGG pathway	Number of genes	Statistical significance[a]
Rap1 signalling pathway	17	P < 0.001
Natural killer cell-mediated cytotoxicity	12	P < 0.001
Thyroid hormone signalling pathway	10	P < 0.01
Malaria	6	P < 0.01
Protein digestion and absorption	8	P < 0.01
Aldosterone-regulated sodium reabsorption	5	P < 0.01
Ras signalling pathway	14	P < 0.05
B cell receptor signalling pathway	6	P < 0.05
FoxO signalling pathway	9	P < 0.05
Pathways in cancer	19	P < 0.05

aHypergeometric distribution.

Rap1, Ras-related protein 1; FoxO, Forkhead box protein O.

The antisense prediction method is recognized as an important strategy for the prediction of the association between lncRNAs and mRNAs.[12] Based on the antisense prediction method, 18 lncRNAs were found to be complementarily paired with mRNAs (Table 4). LncRNA ENSG00000230489 was in the top 20 deregulated lncRNAs and was paired with VAV3 mRNA.

Table 4.

Antisense prediction of long noncoding RNAs (lncRNAs) and mRNA associations.

Number	lncRNA	mRNA
1	ENSG00000179428	IL6
2	ENSG00000224875	EML4
3	ENSG00000238133	MAP3K20
4	ENSG00000230530	LIMD1
5	ENSG00000227117	MTMR3
6	ENSG00000232415	ELN
7	ENSG00000246263	UBR5
8	ENSG00000245105	A2M
9	ENSG00000245648	KLRK1
10	ENSG00000250696	KLRK1
11	ENSG00000250850	PAX5
12	ENSG00000248309	MEF2C
13	ENSG00000261448	XYLT1
14	ENSG00000180139	ACTA2
15	ENSG00000272182	PDHB
16	ENSG00000258168	PTPRB
17	ENSG00000245552	SESN3
18	ENSG00000230489	VAV3

Reverse transcription–qPCR assays were used to confirm the deep sequencing results for the following six lncRNAs: ENSG00000250451, ENSG00000243479, ENST00000230798, ENST00000228630 (upregulated), ENSG00000230489 and ENSG00000227084 (downregulated). Their expression was consistent with the sequencing results (Figure 2a). LncRNAs ENSG00000230489 and ENSG00000227084 were decreased in NPC tissues compared with their levels in CNP tissues (Figures 2b and 2c).

Figure 2.

Validation of the RNA sequencing results. (a) The levels of six long noncoding RNAs (lncRNAs) were evaluated by quantitative polymerase chain reaction in samples from patients with nasopharyngeal carcinoma (NPC; n = 30) and chronic nasopharyngitis (CNP; n = 27) and compared with the RNA sequencing data. (b, c) The levels of lncRNA ENSG00000230489 and ENSG00000227084 in the NPC and CNP groups. All data are represented as the mean ± SD; *P < 0.05 for NPC versus CNP samples; Student’s t-test.

The lncRNAs ENSG00000227084 and ENSG00000230489 were in the top 20 downregulated lncRNAs in NPC (Table 2). Moreover, lncRNAs ENSG00000227084 and ENSG00000230489 might play essential roles in the development of NPC by interfering with the Rap1 signalling pathway and natural killer cell-mediated cytotoxicity, respectively. The potential diagnostic value of lncRNAs ENSG00000230489 and ENSG00000227084 for NPC was evaluated with CNP patients as the control subjects. The areas under the ROC curves for lncRNAs ENSG00000230489 and ENSG00000227084 were 0.9138 and 0.8037, respectively. These results suggested that they have potential diagnostic value for NPC (Figure 3).

Figure 3.

Receiver operating characteristic (ROC) curves for the nasopharyngeal carcinoma samples. (a) ROC curve of ENSG00000230489 (area under the ROC curve = 0.9138); (b) ROC curve of ENSG00000227084 (area under the ROC curve = 0.8037).

The potential diagnostic value of the two lncRNAs was evaluated by a correlation analysis between clinical characteristics and the expression of lncRNAs ENSG00000230489 and ENSG00000227084. This analysis found that a low level of lncRNA ENSG00000230489 but not lncRNA ENSG00000227084 was positively associated with distant metastasis (P = 0.014) (Table 5).

Table 5.

Relationship between the levels of two long noncoding RNAs (lncRNAs) and the clinicopathological characteristics of patients with nasopharyngeal carcinoma (NPC; n = 30).

Characteristics	Patients with NPC n = 30	Levels of ENSG00000230489		Statistical significance[a]	Levels of ENSG00000227084		Statistical significance[a]
		Low n = 15	High n = 15		Low n = 19	High n = 11
Age, years	38.5 ± 14.5
Interquartile range	24–52
≥38.5 years	18	10	8	NS	12	6	NS
<38.5 years	12	5	7		7	5
Sex
Female	7	4	3	NS	5	2	NS
Male	23	11	12		14	9
T classification
T1+T2	16	6	10	NS	10	6	NS
T3+T4	14	9	5		9	5
N classification
N0+N1	15	7	8	NS	9	6	NS
N2+N3	15	8	7		10	5
Distant metastasis
Yes	9	8	1	P = 0.014	7	2	NS
No	21	7	14		12	9

Data presented as mean ± SD, interquartile range and n of patients.

aχ2-test and Fisher’s exact test; NS, no significant association (P ≥0.05).

Since the biological function and clinical value of these two lncRNAs has never been studied, further research is required to elucidate the role of these two lncRNAs in NPC, for example, by determining an lncRNA-mRNA co-expression network (Figure 4).

Figure 4.

Gene network of the correlative genes of long noncoding RNAs (lncRNAs) ENSG00000230489 and ENSG00000227084. The LncRNAs are shown as red nodes while mRNAs are shown as green nodes.

Discussion

Nasopharyngeal carcinoma is one of the most common carcinomas of the head and neck in China.[2] The molecular features underlying the pathogenesis of NPC remain largely unknown. In addition, early detection of NPC is not easy due the lack of obvious clinical features in the early stages.[13] Increasing evidence suggests that lncRNAs play an important role in the progression of cancer by affecting chromatin remodelling and gene expression.[14] At the same time, lncRNAs have been reported to be biomarkers for predicting tumorigenesis and development in the diagnosis of multiple diseases.[15] Previous studies have demonstrated that a few dysregulated lncRNAs, such as AFAP1-AS, ANRIL and LET, are linked to NPC.[16-18] However, the expression signatures and the diagnostic value of lncRNAs in NPC development remain largely unknown.

Deep sequencing has been widely used for clinical detection and has greatly accelerated the discovery and characterization of novel ncRNAs. For example, a previous study demonstrated the expression profiles of microRNAs between NPC and CNP patients by deep sequencing and identified miR-34c as an important gene in NPC tumorigenesis.[19] However, to the best of our knowledge, this current study is the first to use deep sequencing technology to evaluate the genome-wide expression of lncRNAs in tumours and inflammatory tissues derived from the nasopharynx.

The current deep sequencing results were consistent with those of previous studies on the known NPC-related lncRNA ENSG00000228630 (HOTAIR),[20] which is aberrantly upregulated in various types of cancer, including breast, hepatocellular and bladder cancer.[21-23] Moreover, lncRNAs ENSG00000250451 (HOXC-AS1), ENSG00000230798 (FOXD3-AS1) and ENSG00000243479 (MNX1-AS1) have been reported to contribute to tumorigenesis in several tumors;[24-26] and this current study is the first to confirm their expression in NPC, further supporting their potential diagnostic value in tumorigenesis.

The lncRNAs ENSG00000227084 and ENSG00000230489 (VAV3-AS1) are related to the Rap1 signalling pathway and natural killer cell-mediated cytotoxicity, respectively. A previous study showed that increased and aberrant activation of Rap1 signalling can lead to tumour formation and progression in human breast epithelial cells.[27] Since Rap1 signalling has been increasingly recognized as an important pathway linked to cancer biogenesis and metastasis,[28] it is possible to speculate that lncRNA ENSG00000227084 may be a key factor promoting the development of NPC.

Natural killer cell-mediated cytotoxicity plays an important role in the natural immune defence against tumours.[29] According to KEGG pathway analysis, lncRNA ENSG00000230489 might play important roles in the development of NPC by interfering with natural killer cell-mediated cytotoxicity. In addition, antisense prediction showed that lncRNA ENSG00000230489 was associated with VAV3. A previous study reported that some specific noncoding RNAs can regulate the proliferation and metastasis of tumour cells by targeting VAV3.[30]Furthermore, this current study found that NPC patients with distant metastases had lower ENSG00000230489 levels than those without metastases, indicating that lncRNA ENSG00000230489 may function as a potential upstream regulator of VAV3 as a tumour suppressor.

Using ROC curve analysis to predict the potential diagnostic value of the lncRNAs ENSG00000227084 and ENSG00000230489, the current study found that both of these have diagnostic value for NPC. Moreover, the diagnostic value of ENSG00000230489 was higher than that of ENSG00000227084. Whether lncRNAs ENSG00000230489 and ENSG00000227084 could be potential targets for the treatment of NPC requires further investigation.

In conclusion, this current study provides, for the first time, a global lncRNA expression profile in NPC by deep sequencing. Furthermore, this study suggests that lncRNAs ENSG00000230489 and ENSG00000227084 may be potential diagnostic and treatment biomarkers of NPC.