Analysis of Differential Expressions of Long Non-coding RNAs in Nasopharyngeal Carcinoma Using Next-generation Deep Sequencing.

Background: Little knowledge about long non-coding RNAs(lncRNAs) in nasopharyngeal carcinoma (NPC) has been acquired.
Methods: Next-generation sequencing was applied in 7 cases of NPC tissues and 7 cases of normal tissues in nasopharynx. PLEX, CNCI and CPAT soft-wares were used to predict novel lncRNAs. Real-time Quantitative PCR (qPCR) further validated the data in 20 cases of NPC tissues and 14 cases of normal tissues. Then the cis-regulators and trans-regulators and potential biological functions together with pathways were predicted by Bioinformatics.
Results: Totally, 4248 novel lncRNAs were found to be expressed in our samples. And 2192 lncRNAs and 23342 mRNAs were considered to be differentially expressed in NPC. Among the results, 306 lncRNAs and 4599 mRNAs were significantly up-regulated, whereas 204 lncRNAs and 2059 mRNAs were significantly down-regulated, respectively. Moreover, 62 lncRNAs trans-regulated genes were involved in Epstein-Barr virus (EBV) infection pathway in our study. Jun proto-oncogene (JUN), which was related to a cis-regulator lncRNA RP4-794H19.1, was enriched in cancers and involved in Tumor Necrosis Factor (TNF) signaling pathway, might play a key role in NPC.
Conclusion: These findings broadened the lncRNAs landscape of NPC tissues and shed light on the roles of these lncRNAs, which might be conducive to the comprehensive management of NPC.

Introduction

Globally, about 86,500 cases of nasopharyngeal carcinoma (NPC) were reported in 2012. NPC was a mortal disease prevailing in South-east Asia and North Africa 1. In China, the standardized mortality rate was 0.78 /105 issued in 20132. NPC arose as a result of its various genetic susceptibility, environmental or random factors and Epstein-Barr virus (EBV) infection3. Albeit researchers were at full tilt unraveling the tumorigenesis and the pathways, we seldom knew them yet.

Long non-coding RNAs(lncRNAs) previously were considered as “transcriptional noise”, which accounted for at least 75% of DNA transcripts4. Generally speaking, lncRNAs are longer than 200 nucleotides and short of protein-coding capacity, with lacking of open reading frame (ORF)5, 6. Encyclopedia of DNA Elements (ENCODE) Project Consortium issued at least 120,000 transcripts that served as lncRNAs in the human genome7. Now lncRNAs have been validated that it played a gorgeous role in regulating gene expression at various levels to affect cellular homeostasis, including proliferation, survival, migration and genomic stability4. Up-to-date evidences laid bare that the aberrantly expressed lncRNAs always were relative to the development and progression of cancers8-11. Given the eagerness to better understand NPC, many researchers took the cutting-edge lncRNAs to the spotlight.

The lncRNAs research in NPC usually resorted to Microarray and lncRNA-protein-based relationship technology 12-16, or at cell level17. Microarray was a high-throughput sequencing, but it was limited by countable numbers of lncRNAs, and all the probes were already known, which was the reason for refraining the technology from discovering novel lncRNAs. Despite LncRNA-protein-based relationship technology as RNA Immunoprecipitation (RIP)-chip and RIP-sequencing, its antibody specificity made the research focus on one point. The next-generation deep sequencing took its advantages in discovering massive exiting and novel lncRNAs18. Here, we applied the next-generation deep sequencing to provide a comprehensive view for insights into NPC tissues.

Materials and Methods

Samples collection

There were 20 cases of NPC tissues and 14 cases of normal tissues collected by the second Xiangya Hospital of Central South University. All sample tissues were confirmed by histopathological examination in the second Xiangya Hospital of Central South University. This study was approved by the Ethical Committee of Xiangya Hospital, Central South University. All patients had voluntarily signed a written consent before samples were collected. The acquired tissue samples were placed into liquid nitrogen immediately and then were transferred to a -80°C refrigerator.

RNA extraction

The abstraction of total RNAs derived from 20 cases of NPC tissues and 14 cases of normal tissues, using TRIzol (Life Technology, California, USA) according to the guideline protocol. The sample tissues were placed in 2ml RNA-free EP tube that grinded with four 4mm steel balls in TissueLyser LT (QIAGEN, Dusseldorf, Germany). The purity and concentration of RNAs were measured by OD260/280 readings shown on DeNovix DS-11+ (DeNovix, Delaware, USA), which ranging from 1.8-2.2 were kept. RNAs integrity was qualified by standard denaturing agarose gel electrophoresis.

cDNA synthesizing and Real-time Quantitative PCR(qPCR) analysis

The first strand of cDNA was synthesized by Promega GoScript (Promega, Madison, USA) with 2 μg total RNA. The qPCR setting 40 cycles of amplification was performed on GoTaq qPCR Master Mix (Promega, Madison, USA) on LightCycler 480 Ⅱ QPCR System (Roch, Basel, Switzerland). β-actin served as a reference, and relative-expression levels of lncRNAs were quantified through 2(-∆∆Ct) methods.

Library preparation for lncRNA sequencing

There were 7 cases of NPC tissues and 7 cases of normal tissues for RNA sequencing. Total RNAs were abstracted as previously, strand-specific RNA-seq Libraries were generated with TruSeq® RNA LT Sample Prep Kit v2 (Illumina, CA, USA) following the handbook, subsequently the Ribo-Zero Gold kit (Illumina, CA, USA) was used for the removing of rRNA, and then the fragmentations of RNAs were harvested. RNA fragmentations were used to synthesize double-stranded cDNA with substituting dUTPs for dTTPs in the second strand cDNA synthesis. Next, double-stranded cDNA were added to A-tailed and ligated to adaptors. Finally, PCR amplification was performed to enrich cDNA libraries. PCR products were qualified by Quant-iT™ PicoGreen® dsDNA Assay Kit (Life Technology, MA, USA) for RNA sequencing.

RNA sequencing

Sequencing was performed on an Illumina Hiseq 2500 instrument (Illumina, CA, USA) with the TruSeq PE Cluster Kit v3-cBot-HS (Illumina, CA, USA).

Results

The prediction of novel lncRNAs

Transcripts were aligned with databases, such as the Kyoto Encyclopedia of Genes and Genomes (KEGG), Clusters of Orthologous Groups (COG) and Swissprot Database, then the mapped transcripts (mRNA) were dropped. Subsequently, transcripts annotated in the lncRNA library (Ensemble) were also cast away. The remaining transcripts were inputted into three Coding Potential Calculator soft-wares (PLEX, CNCI and CPAT) to predict coding and non-coding sequences. We kept the 4248 transcripts which overlapped in three calculation methods, serving as novel lncRNAs. CNCI-CPAT, CPAT-PLEK and PLEK-CNCI overlapped numbers were 237, 26 and 401, respectively. Venn Diagram of predicted novel lncRNAs was shown in Figure 1.

Figure 1

Venn Diagram of predicted novel lncRNAs with three the Coding Potential Calculator soft-wares.

The landscape of differential expressions of lncRNAs and mRNAs sequencing in NPC and normal tissues

In total, 2192 lncRNAs and 23342 mRNAs were considered to be differentially expressed in NPC. Among the results, 204 lncRNAs and 2059 mRNAs were significantly down-regulated, and 306 lncRNAs and 4599 mRNAs were significantly up-regulated (≥2 fold change, FDR<0.05). Compared with normal tissues, OSER1-AS1 (fold change: -4.54, FDR<0.05) and MUC4 (fold change: -8.00, FDR<0.05) were the most down-regulated lncRNA and mRNA, as opposed to the most up-regulated lncRNA and mRNA, which were RP11-624L4.1 (fold change: 4.35, FDR<0.05) and XLOC.42541.1 (fold change:13.20, FDR<0.05), respectively. The volcano plot of differential expressions of lncRNAs and mRNAs were displayed in Figure 2. Heat maps of differential expressions and hierarchical clustering of lncRNAs and mRNAs in NPC and normal tissues were demonstrated in Figure 3.

Figure 2

(A) The volcano plot of differential expressions of lncRNAs in NPC and normal tissues. Dark spot: log2(fold change)<1, P value>0.05; Red spot: log2(fold change)≥1, P value≤0.05; Green spot: log2(fold change)≤-1, P value≤0.05. (B) The volcano plot of differential expressions of mRNAs in NPC and normal tissues. Dark spot: log2(fold change)<1, P value>0.05; Red spot: log2(fold change)≥1, P value≤0.05; Green spot: log2(fold change)≤-1, P value≤0.05.

Figure 3

(A) Heat maps of differential expressions of and hierarchical clustering of lncRNAsin NPC and normal tissues. (B) Heat maps of differential expressions of and hierarchical clustering of mRNAs in NPC and normal tissues.

RNA sequencing data validation by qPCR

Seven lncRNAs were sifted through on the basis of differentially expressed lncRNAs prediction for validating RNA sequencing data, which were lncRNA-NR2F2 antisense RNA1(NR2F2-AS1) (fold change:-3.2), lncRNA-family with sequence similarity 95 member C(FAM95C) (fold change:-3.3), lncRNA- long intergenic non-protein coding RNA 1106 (LINC01106) (fold change:-3.7), lncRNA-CH507- 513H4.6 (fold change:3.8), lncRNA-THAP9 antisense RNA 1(THAP9-AS1) (fold change:3.7), lncRNA- CH507-513H4.3 (fold change:3.8) and lncRNA-RP4- 794H19.1 (fold change:1.8). QPCR results were consistent with RNA sequencing data(Table 1). The lncRNAs in NPC reported in former research also credited our RNA-seq data(Table 2). QPCR results were shown in Figure 4.

Table 1

Primers used for qPCR analysis of lncRNA levels.

LncRNA	Primer	Expression	Product Length (bp)
NR2F2-AS1	F: GGCTTGGTTGGACCGTACTCR:ACTGCCTGCCACTCTTGTGA	Down	121
FAM95C	F: CCTCATCCAGAGCAGGACAAR: GCCTGCATTAGCATCTGTGC	Down	128
LINC01106	F:GAGAGAGGAATGGCTGCATGTR:TCAGGAGGTTGCAGGTTGAC	Down	121
CH507-513H4.6	F: CATGCCAGAGTCTCGTTCGTTR: GGACACGGACAGGATTGACA	Up	125
THAP9-AS1	F: GAAGGTCCACAGCCAGTTCCR: CGCATTGCTGTAGGCATAGG	Up	121
CH507-513H4.3	F:TTGCTCAATCTCGGGTGGCTR:TGGTGGAGCGATTTGTCTGG	Up	138
RP4-794H19.1	F:CCTTCCGCCTTGGCAAAGATR:TGTGGGGTGAAGGTCTGGTT	UP	120

Table 2

LncRNAs reported in NPC in former research.

LncRNA	Location	Expression	Fold change
xist21	Chromosome X73820651-73852753	UP	1.5
Hotair20	Chromosome 1253962308-53974956	UP	1.2
CASC923	Chromosome 875223117-75278889	UP	1.6
MEG322	Chromosome 14100779410-100861031	DOWN	2.6
EWSAT124	Chromosome 1569072926-69095820	UP	1.7

Figure 4

The relative-expression levels of seven lncRNAs (A, NR2F2-AS1; B, FAM95C; C, LINC01106; D, CH507-513H4.6; E, THAP9-AS1, F, CH507-513H4.3; G, RP4-794H19.1) in NPC and normal tissues. Kruskal-Wallis H test was applied, *p<0.05.

Differential expressions of lncRNAs and Cis-regulation and trans-regulation genes co-expression networks

We defined the chromosomal genes (co-expression mRNA mapped), which were 100 kb upstream or downstream of differentially expressed lncRNAs, as potential lncRNA cis-regulatory genes. Totally, there were 1895 differentially expressed lncRNAs with 4848 cis-regulatory genes. Among them, LA16c-313D11.9 had 231 cis-regulatory genes(the most). One differentially expressed lncRNA for one cis-regulatory gene number was 59. One cis-regulatory gene for one differentially expressed lncRNA was 566.

Pearson correlation coefficient(R2≥0.95) was applied to predict trans-regulatory genes. In total, there were 1284 differentially expressed lncRNAs with 6437 trans-regulatory genes. Among them, AC000123.4 had 577 trans-regulatory genes(the most). One differentially expressed lncRNA for one trans- regulatory gene number was 2263. One trans- regulatory genes for one differentially expressed lncRNA was 566. The networks of cis-regulation and trans-regulation were illuminated in Figure 5.

Figure 5

(A) Cis-regulation of the networks. Red spot represented differential expression of lncRNA. Blue spot represented cis-regulation gene. (B) Trans-regulation of the networks. Red spot represented differential expression of lncRNA. Blue spot represented trans-regulation gene.

Gene ontology (GO) enrichment and pathway analysis

GO enrichment analysis was conducted by TopGO software, and associated with Fisher Exact Test. The results suggested that the significant functions mainly focused on embryonic skeletal system development, embryonic skeletal system morphogenesis, ATP-dependent microtubule motor activity and regulation of gene expression in lncRNAs cis-regulation genes between NPC and normal tissues. As for lncRNAs trans-regulation genes, they laid in intracellular part, intracellular and nuclear part. The detailed results were presented in Figure 6. A and Figure 7. A.

Figure 6

(A) GO enrichmentanalysis of differential expressions of lncRNAs cis-regulation genes. (B) KEGG pathway analysis of differential expressions of lncRNAs cis-regulation genes.

Figure 7

(A) GO enrichment of differential expressions of lncRNAs trans-regulation genes (B) KEGG pathway analysis of differential expressions of lncRNAs trans-regulation genes.

KEGG pathway analysis concluded that the most significant pathways of lncRNAs cis-regulation genes were Proteoglycans in cancer, Vascular smooth muscle contraction and Glutamatergic synapse, while Proteasome, Spliceosome and Huntington's disease were in lncRNAs trans-regulation genes. The rest were demonstrated in Figure 6. B and Figure 7. B.

Discussion

In 2015, there were sixty thousand new NPC cases reported in China19. For NPC patients usually were diagnosed at late stage, and it's a challenge to determine the disease at an early stage. Here, we resorted to the next-generation sequencing technology in 7 cases of NPC tissues and 7 cases of nasopharyngeal normal tissues. 2192 differentially expressed lncRNAs were found, which had the making of potential biomarkers for early-stage diagnosis, progression and prognosis20. Moreover, 4248 novel lncRNAs were obtained, and might tremendously enrich the human lncRNA pool.

There was an interesting finding when we mapped differentially expressed lncRNAs to human chromosomes. Y- chromosome had no differentially expressed lncRNAs in NPC and normal tissues(Figure 8), meanwhile, no such studies have been issued. According to demographic trends, men were two to three fold more likely to acquire NPC than women were1. Maybe Y- chromosome was short of some suppressor lncRNAs inflicting a higher morbidity compared with women in NPC.

Figure 8

Distribution of differential expressions of lncRNAs and mRNAs along each chromosome. The chromosome constitution (outer circle, depicted in blue), distribution of differential expressions of mRNAs (middle circle, depicted in red: up-regulation, green: down-regulation) and distribution of differential expressionsof lncRNAs (inner circle, depicted in red: up-regulation, green: down-regulation).

Seven qPCR results were credible based on our work, and several reported studies also supported it. Increased X inactive specific transcript (XIST) expression was associated with poor prognosis and stunted cell growth21. HOX transcript antisense RNA(Hotair) might play an oncogenic role in NPC progression, and its knockdown attenuated NPC cell growth and angiogenesis20. Both DNA copy number loss and aberrant promoter methylation led to maternally expressed 3(MEG3) inactivation. As a tumor suppressor, which could be significantly inhibited in cell proliferation, colony formation22. Overexpressed cancer susceptibility 9 (CASC9) promoted the glycolysis and tumorigenesis of NPC cells via activating hypoxia inducible factor 1 alpha subunit (HIF1alpha)23. Ewing sarcoma associated transcript 1(EWSAT1) could up-regulate the expressions of miR-326/330-5p clusters targeted gene cyclin D1 through a competitive 'sponge' of miR-326/330-5p clusters, which might promote the development of NPC24.

However, only few known lncRNAs had functional annotations. The methods for lncRNAs function prediction always investigated their correlated protein coding genes and related biological pathways25. It also had been proved that lncRNAs function through cis-regulation of nearby protein coding genes was a common mechanism26. Only a few lncRNAs were refined as cis-regulators, while the majority of functionally annotated ones were trans-regulators27. Therefore, we divided the differentially expressed lncRNAs cis-regulation or trans-regulation genes into two categories to explore the potential target genes, GO enrichment and KEGG pathways analysis. 62 lncRNAs trans-regulated genes were involved in the EBV infection pathway in our study (Figure 7B), as CD4428 and interleukin 1 receptor associated kinase 1(IRAK1)29 (Table S1). In our work, we defined 100kb distance as a criterion between lncRNAs and cis-regulated genes. Jun proto-oncogene(JUN), which was located 98542b downstream to LncRNA RP4-794H19.1, was enriched in cancers and involve in Tumor Necrosis Factor (TNF) signaling pathway. For example, the expression levels of JUN family members were positively correlated to the cell proliferation in diffuse large B-cell lymphomas30. Furthermore, NAG7 through the activity of JUN pathways could promote the invasion in NPC31. For good measure, the EBV-encoded EBNA1 regulated the JUN transcription factor pathway in NPC and increased angiogenesis in vitro32. Our results underpinned that lncRNA RP4-794H19.1 in relation to JUN might play a key role in NPC. However our work strongly depended on the theoretically computerized calculation, lacking further numerous practice validation. These findings broadened the lncRNAs landscape of NPC and shed light on the roles of these lncRNAs, which might be conducive to the comprehensive management of NPC.

Supplementary table S1.

Click here for additional data file.

Supplementary table S2.

Click here for additional data file.

Supplementary table S3.

Click here for additional data file.

Supplementary table S4.

Click here for additional data file.

Supplementary table S5.

Click here for additional data file.

Supplementary table S6.

Click here for additional data file.

Supplementary table S7.

Click here for additional data file.

Supplementary table S8.

Click here for additional data file.

Supplementary table S9.

Click here for additional data file.