Long Noncoding RNAs and Messenger RNAs Expression Profiles Potentially Regulated by ZBTB7A in Nasopharyngeal Carcinoma

Our previous studies showed that ZBTB7A played an important role in promoting nasopharyngeal carcinoma (NPC) progression. However, molecular mechanisms of different levels of ZBTB7A are still unclear. It is necessary to search molecular markers which are closely connected with ZBTB7A. We selected NPC sublines CNE2 with stably transfecting empty plasmid (negative control, NC) and short hair RNA (shRNA) plasmid targeting ZBTB7A as research objectives. Microarray was used to screen differentially expressed long noncoding RNAs (lncRNAs) and messenger RNAs (mRNAs) via shRNA-CNE2 versus NC-CNE2. Quantitative PCR (qPCR) was used to validate lncRNAs and mRNAs from the sublines, chronic rhinitis, and NPC tissues. Bioinformatics was used to analyze regulatory pathways which were connected with ZBTB7A. The 1501 lncRNAs (long noncoding RNAs) and 1275 differentially expressed mRNAs were upregulated or downregulated over 2-fold. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed that the upregulated or downregulated carbohydrate and lipid metabolisms probably involved in carcinogenicity of shRNA-CNE2 (P-value cut-off was 0.05). In order to find the molecular mechanisms of ZBTB7A, we validated 12 differentially expressed lncRNAs and their nearby mRNAs by qPCR. Most of the differentially expressed mRNAs are closely connected with carbohydrate and lipid metabolisms in multiply cancers. Furthermore, part of them were validated in NPC and rhinitis tissues by qPCR. As a result, NR_047538, ENST00000442852, and fatty acid synthase (FASN) were closely associated with NPC. ZBTB7A had a positive association with NR_047538 and negative associations with ENST00000442852 and FASN. The results probably provide novel candidate biomarkers for NPC progression with different levels of ZBTB7A.


Introduction
Nasopharyngeal carcinoma (NPC) is an endemic malignant head and neck tumor in southern China [1,2].Epstein-Barr virus (EBV) and tumor metastasis related factors are closely connected with NPC progression [3,4].The biomarkers of EBV include Epstein-Bar encoded small nuclear RNA (EBER) in tissues, EBV DNA, and anti-EBV antibodies in circulating plasma or serum.All of them are important predictors in early diagnosis and prognosis of NPC.Furthermore, plasma EBV DNA is a better marker in advanced NPC than others [5][6][7].However, a few patients have negative EBV DNA in plasma [8,9].The negative result indicates a limitation of EBV DNA.In order to improve detection rate of NPC, it is imperative to search novel candidate biomarkers.
Therefore, it is meaningful to discover the unknown field.We will screen differentially expressed lncRNAs and messenger RNAs (mRNAs) from NPC cell sublines with stably transfecting empty plasmid and short hair RNA (shRNA) plasmid targeting ZBTB7A by lncRNA microarray.Bioinformatics will be used to search potential connections between ZBTB7A and relevant pathways.Part of them will be selected and validated by qPCR in the sublines, chronic rhinitis, and NPC tissues.

Patient Selection.
After being approved by patients, 20 chronic rhinitis and 60 NPC tissues were obtained using biopsy (July 2016-February 2017).All of them were preserved in liquid nitrogen.The tissues were confirmed by pathology from the formalin-fixed paraffin wax-embedded samples.The pathological type of NPC was undifferentiated carcinomas of the nasopharyngeal type (UCNT).This study was approved by the ethics committee of The People's Hospital of Guangxi Zhuang Autonomous Region.

Cell
Culture.CNE2 and CNE3 were preserved in Research Center of Medical Sciences, The People's Hospital of Guangxi Zhuang Autonomous Region (Nanning, Guangxi, China).CNE2 was established from a patient suffering from NPC in Guangdong Province.The histological type of NPC was a poorly differentiated squamous carcinoma [38].CNE3 was established from a patient suffering from liver metastatic carcinoma of primary NPC in Guangxi Province.The histological type of liver metastatic carcinoma and CNE3 was a poorly differentiated adenosquamous carcinoma [39], but that of CNE3 turned poorly differentiated adenocarcinoma after 20 years [40].As a team member of Professor Yi Zeng, Wei Jiao has cultured and detected the cell lines in our hospital [23-25, 41, 42].5-8F and 6-10B were kindly provided by Professor Musheng Zeng (State Key Laboratory of Oncology in South China, Sun Yat-Sen University Cancer Center, Guangzhou, Guangdong, China).They were grown in RPMI Medium 1640 basic with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA).NP69 was kindly provided by Professor Sai-Wah Tsao (School of Biomedical Sciences, University of Hong Kong, Hong Kong SAR).NP69 was in Keratinocyte-SFM with 5% Bovine Pituitary Extract and Recombinant Epidermal Growth Factor (Gibco).All of them were in supplementary material.Based on the prior study [24], we used CNE2 sublines with stably transfecting empty plasmid and shRNA plasmid targeting ZBTB7A to screen differentially expressed lncRNAs and mRNAs.The sublines were named as negative control-CNE2 (NC-CNE2) and shRNA-CNE2, of which tumorigenicity was stronger than that of NC-5-8F and NC-6-10B [24].
2.3.LncRNA Microarray.shRNA-CNE2 and NC-CNE2 were used in lncRNA microarray.Both of them were acquired from three different passages.Total RNA was extracted by TRIZOL (Invitrogen, Carlsbad, CA, USA).RNA quantity and quality were measured by NanoDrop ND-1000 (Thermo Fisher Scientific, Inc.).RNA integrity was assessed by standard denaturing agarose gel electrophoresis.The Arraystar Human lncRNA Microarray V3.0 was designed for the global profiling of human lncRNAs and protein-coding mRNAs.30,586 lncRNAs and 26,109 mRNAs could be detected by third-generation lncRNA microarray.

RNA Labeling and Array
Hybridization.Sample labeling and array hybridization were performed according to the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technology, Santa Clara, CA 95051, USA) with minor modifications.Briefly, mRNA was purified from total RNA after removal of ribosome RNA (rRNA; mRNA-ONLY6 Eukaryotic mRNA Isolation Kit; Epicentre, Madison, WI 53719, USA).Each sample was amplified and transcribed into fluorescent complementary RNA (cRNA) along the entire length of the mRNAs without 3  bias utilizing a random priming method (Arraystar Flash BioMed Research International 3 RNA Labeling Kit; Arraystar, Inc., Rockville, MD 20850, USA).The labeled cRNAs were purified by RNeasy Mini Kit (Qiagen, Düsseldorf, Nordrhein-Westfalen, Germany).The concentration and specific activity of the labeled cRNAs (pmol Cy3/g cRNA) were measured by NanoDrop ND-1000. 1 g of each labeled cRNA was fragmented by adding 5 l 10 x Blocking Agent and 1 l of 25 x Fragmentation Buffer, then the mixture was heated at 60 ∘ C for 30 min, finally 25 l 2 x GE hybridization buffer was added to dilute the labeled cRNA.50 l of hybridization solution was dispensed into the gasket slide and assembled to the lncRNA expression microarray slide.The slides were incubated for 17 h at 65 ∘ C in an Agilent Hybridization Oven.The hybridized arrays were washed, fixed, and scanned with using the Agilent DNA Microarray Scanner (part number G2505C; Agilent Technology).

Data Analysis.
The Agilent Feature Extraction software (version 11.0.1.1)was used to analyze acquired array images.Quantile normalization and subsequent data processing were performed with using the GeneSpring GX v12.1 software package (Agilent Technologies).Differentially expressed lncRNAs and mRNAs between the two groups were identified using fold-change >2 as the cut-off.Hierarchical clustering and combined analysis were performed using homemade scripts.Gene ontology (GO) and pathway analysis were performed in the standard enrichment computation method.

Validation of the Differentially Expressed LncRNAs and
MRNAs by QPCR.The total RNA of NC-CNE2 and shRNA-CNE2 cells, chronic rhinitis, and NPC tissues was extracted by TRIZOL (Invitrogen) and reversely transcribed by Super-Script 6 III Reverse Transcriptase Kit (Invitrogen).12 lncR-NAs and their nearby mRNAs were selected basing on data analysis.All of them were differentially expressed.The quantitative PCR (qPCR) was executed with 2 × PCR master mix (Arraystar) by ViiA 7 qPCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.).Specific primers of the lncRNAs and nearby mRNAs were designed by Primer 5.0 (Table 1).The reaction program consisted of an initial denaturation step at 95 ∘ C for 10 minutes; denaturation at 95 ∘ C for 10 seconds, and annealing at 60 ∘ C for 60 seconds for 40 cycles; the dissociation stage at 95 ∘ C for 10 seconds, 60 ∘ C for 1 minute, 95 ∘ C for 15 seconds, and 60 ∘ C for 15 seconds.Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an internal control.The experiments were performed in triplicate.The fold-change of the lncRNAs and mRNAs were calculated by the 2 Delta-delta cycle threshold (Ct) method.

Statistical Analysis.
Statistical analysis was performed by SPSS 13.0 and sigmaPlot 12.5 (SPSS Inc., Chicago, IL, USA), GraphPad Prism 5.0 (GraphPad Software, San Diego, CA, USA).The results of the assays were presented as mean ± Standard Deviation (SD).Statistical differences between two groups were evaluated with independent samples t-test.In GO analysis, Fisher's exact test was used to search more overlap between the differentially expressed list and the GO annotation list than that of being expected by chance.In scatter plot of tissues by qPCR, Mann Whitney U test was used if the variances were significantly different through F test.P<0.05 was considered to be statistically significant.

Differentially Expressed LncRNAs and MRNAs in NPC
Cells with Low Expression of ZBTB7A.After quantile normalization of the raw data, the expression profiles of 23947 lncRNAs and 21058 mRNAs were obtained from shRNA-CNE2 and NC-CNE2 (Figure 1).Cluster analysis showed some lncRNAs and mRNAs were differentially expressed via shRNA-CNE2 versus NC-CNE2 (Figures 1(a) and 1(b)).Box plot showed the distribution of normalized intensity values of test and control groups was generally symmetrical (Figures 2(c) and 2(d)).Scatter plot figuratively showed relative folds of 1501 lncRNAs (738 upregulated and 763 downregulated) and 1276 mRNAs (679 upregulated and 597 downregulated), which were over 2-fold.The number of differentially expressed lncRNAs and mRNAs was listed.We divided them into 3 groups according to different foldchange.There were 2-5-fold-change, 5-10-fold-change, and more than 10-fold-change.77 lncRNAs and 44 mRNAs were over 10-fold-change (Table 2).

GO Analysis. GO is a system of standard classification of gene function, including biological process (BP), cellular component (CC), and molecular function (MF). The results
showed the significant GO terms of differentially expressed genes.We found that the highest enriched GO terms of upregulated mRNAs were response to stress (GO: BP), cytoplasmic part (GO: CC), and oxidoreductase activity (GO: MF) (Figure 2(a)).The highest enriched GO terms of downregulated mRNAs were response to stress (GO: BP), cytoplasm (GO: CC), and catalytic activity (GO: BF) (Figure 2(b)).The results indicated some differentially expressed lncRNAs and relevant mRNAs probably involved in the activities.

KEGG Pathway Analysis. Kyoto Encyclopedia of Genes
and Genomes (KEGG) is a database resource for knowing high-level biological functions and utilities.KEGG pathways were constructed to better know the biological pathways and search the molecular mechanisms of NPC progression.The results showed significant pathways of differential genes.KEGG pathways connected with upregulated (Figures 3(a) and 3(b), left) mRNAs and downregulated (Figures 3(a) and 3(b), right) differentially distributed mRNAs.The top upregulated pathway was steroid hormone biosynthesis.Steroid is a kind of lipid.The top downregulated pathway was protein digestion and absorption.The third downregulated pathway was peroxisome proliferators-activated receptor (PPAR) signaling pathway, which included important lipid metabolism.463 differentially expressed mRNAs were closely associated with carbohydrate and lipid metabolisms.The pathways of lipid metabolism probably were more meaningful than those of carbohydrate metabolism because of more enrichment scores in lipid pathways.In upregulated pathways, the scores  of steroid hormone biosynthesis and steroid biosynthesis were 5.04301 and 3.129842.In downregulated pathways, the score of PPAR signaling pathway was 2.035416 (all p<0.05) (Table 3).PPAR signaling pathway was markedly suppressed in the cells.Fatty acid binding protein 4 (FABP4) was markedly downregulated.The differentially expressed mRNAs of the pathways were associated with lipid metabolism, such as AKR1C3 and FABP4 (Figure 4).

QPCR Validation. lncRNAs mainly include antisense lncRNAs, enhancer lncRNAs, and long intergenic noncoding
RNAs (lincRNAs).They have different biological functions through regulating nearby mRNAs [43][44][45].In order to validate authenticity and reliability of lncRNA microarray, 12 differentially expressed lncRNAs and their nearby mRNAs were detected by qPCR.The differentially expressed mRNAs are closely connected with carbohydrate and/or lipid metabolisms of cancers, which means the lncRNAs probably interact with cancers via the mRNAs.The results of qPCR were consistent with those of lncRNA microarray through shRNA-CNE2 versus NC-CNE2.All of them have shown the same upregulation or downregulation trend (Figure 5).6 lncRNAs (NR 033967, ENST00000398216, uc001enh.1,TCONS 00019671, TCONS 00025256, TCONS 00029013) and 1 mRNA (integrin subunit alpha 5, ITGA5) would not be used in next assay because the results of qPCR were not considered to be statistically significant (P>0.05).

Discussion
lncRNA microarrays have been used in studies of multiple cancers [27,46,47].Differentially expressed lncRNAs and mRNAs can be found by the method, which reveals some lncRNAs are important regulatory factors in melanoma, renal tumor, hepatocellular carcinoma, and colorectal cancer [48][49][50][51].lncRNA ANRIL, also called NR 047538, promotes NPC progression via reprogramming cell glucose metabolism [31].lncRNA HOTAIR promotes tumorigenesis of NPC via upregulating FASN, which is a pivotal enzyme of lipid metabolism [52].Therefore, the technique is effectively applied in NPC study of exploring possible connections between ZBTB7A and differentially expressed lncRNAs.
We went on disclosing that carbohydrate and lipid metabolisms probably involved in NPC progression through KEGG analysis.In order to search the molecular mechanisms of ZBTB7A, we selected 12 differentially expressed lncRNAs and their nearby mRNAs to validate.12 differentially expressed mRNAs included FASN, CDKN2B, SLC2A1, ALDOC, ISG15, IER3, TXNIP, KAT5, ITGA5, DYRK1A, SIK1, and PON3.Basing on the bioinformatics analysis, the lncRNAs and mRNAs are possibly associated with cancer cells progression via carbohydrate or lipid pathway.In order to exclude the contaminated possibility by Hela cell, parts of them were validated in NPC and rhinitis tissues by qPCR.lncRNA NR 047538, lncRNA ENST00000442852, and fatty acid synthase (FASN) were closely associated with NPC.ZBTB7A had a positive association with NR 047538, while having negative associations with ENST00000442852 and FASN.The results provide novel candidate biomarkers for NPC progression with different levels of ZBTB7A.
LncRNA NR 047538 and FSAN were important markers and regulatory factor in NPC [31,74].LncRNA NR 047538 can promote NPC progression by upregulating GLUT1 expression [31].It is also associated with lipid pathways [75].Sex determining region Y-box 2 gene (SOX2) induces proliferation of NPC cells through activating lncRNA NR 047538 [76].Downregulation of lncRNA NR 047538 inhibits NPC tumorigenicity and enhances the efforts of radiotherapy and chemotherapy via regulating microRNA 125a and let-7a [77,78].FASN is an important enzyme of lipogenesis.It promotes NPC progression through remarkably providing endogenous fatty acid [52].Interestingly, we found that sterol regulatory element binding transcription factor 1 (SREBF1) also was upregulated by lncRNA microarray and qPCR (Supplementary Materials, Figure S4 and Table S1).SREBF1 is the transcript of sterol regulatory element binding protein 1 (SREBP1).Activated SREBP1 enters nucleus and transcribes the genes of lipid metabolism such as FASN [79].EBVencoded latent membrane protein 1 (LMP1) induces cell proliferation and NPC metastasis via activating SREBP1 and its downstream FASN [80].
The potential biomarker lncRNA ENST00000442852 has not been reported.We speculate it may be connected with the nearby mRNA IER3, which is also called immediate early response gene X-1 (IEX-1).It is immediately regulated by transcriptional factors, inflammatory cytokines, growth factors, and so on [81].It is connected with poor or good prognosis in different types of cancers [82,83].Therefore, the functions and characteristics of lncRNA ENST00000442852 will be explored in NPC.
Although SLC2A1, CDKN2B, and ISG15 are closely connected with NPC progression in other studies [31,84,85], our results do not show the positive results.The quantity of 20 chronic rhinitis and 60 NPC tissues is limited.The genotyping of Guangxi population probably is different from that of other region.The protein expression of the genes should be validated by Western blot or immunohistochemistry in future because of the difference between mRNA and protein expression of the same gene.

Conclusions
Above all, the results showed the changes of the differentially expressed lncRNAs and mRNAs with stable loss of ZBTB7A expression in shRNA-CNE2 cells.The results may provide potential biomarkers for NPC progression with different levels of ZBTB7A, such as lncRNA NR 047538, lncRNA ENST00000442852, and FASN.In the future, we will deeply dig the connections between the lncRNAs/mRNAs and ZBTB7A in NPC.

Figure 1 :
Figure 1: RNA expression profiles in shRNA-CNE2 and NC-CNE2 cells.Hierarchical clustering indicated (a) lncRNA and (b) mRNA profiles.Box plot of (c) lncRNA or (d) mRNA profiles was a traditional method for visualizing separately the distribution of dataset.Scatter plot of (e) lncRNA or (f) mRNA profiles was a convenient visualization for assessing the variation of chips.T1-T3 meant test groups which were stably transfected with shRNA plasmid targeting ZBTB7A.C1-C3 meant negative control groups which were stably transfected with blank plasmid.'Red' showed relatively high expression.'Green' showed relatively low expression.The values of the x-and y-axes were the average normalized signal values of NC-CNE2 and shRNA-CNE2 in the scatter plot (log 2 scaled).

3. 4 .Figure 2 :
Figure 2: GO (gene ontology) analysis.(a) GO upregulated.(b) GO downregulated.The P value denoted the significance of GO terms enrichment in the DE mRNAs.The lower the p-value, the more significant the GO Term (P value ≤0.05 was recommended).DE, differentially expressed; BP, biological process; CC, cellular component.MF, molecular function.* It meant oxidoreductase activity, acting on the CH-CH group of donors, NAD or NADP as acceptor.* * It meant oxidoreductase activity, acting on the CH-CH group of donors.* * * It meant oxidoreductase activity, acting on the aldehyde or oxo group of donors.

Figure 3 :
Figure 3: KEGG pathways in shRNA-CNE2 compared to NC-CNE2.(a) The most significantly enriched KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways were illustrated by bubbles.The y-axis indicated the enrichment factors, which refers to the upregulated (left) and downregulated (right) differentially expressed mRNA number to the total mRNAs in a certain pathway.The size of bubble indicated mean number of mRNA enriched in a given pathway.The colour of bubble indicated P value (adjusted p value).(b) Bar plots represented enriched KEGG pathways associated with upregulated (left) and downregulated (right) DE mRNAs.The y-axis indicated average change fold of the implicated mRNAs.Sig, significant; DE, differentially expressed.

Figure 7 :
Figure 7: ZBTB7A were separately connected with NR 047538, FASN, and ENST00000442852 by scatter plot.ZBTB7A had a positive association with (a) FASN and negative associations with (b) NR 047538 and (c) ENST00000442852.p<0.05 was considered to be statistically significant.

Table 1 :
The list of specific primers of 12 differentially expressed lncRNAs and nearby mRNA designed utilizing primer 5.0.
Seqname: the sequence identifier; GeneSymbol: the symbol of the lncRNA or mRNA.

Table 4 :
The p values of the differentially expressed lncRNAs and nearby mRNAs by microarray and qPCR.