Uncovering the ceRNA network and DNA methylation associated with gene expression in nasopharyngeal carcinoma

Objective This study aimed to uncover abnormally expressed genes regulated by competitive endogenous RNA (ceRNA) and DNA methylation nasopharyngeal carcinoma and to validate the role of lncRNAs in the ceRNA network on nasopharyngeal carcinoma progression. Methods Based on the GSE64634 (mRNA), GSE32960 (miRNA), GSE95166 (lncRNA), and GSE126683 (lncRNA) datasets, we screened differentially expressed mRNAs, miRNAs and lncRNAs in nasopharyngeal carcinoma. A ceRNA network was subsequently constructed. Differentially methylated genes were screened using the GSE62336 dataset. The abnormally expressed genes regulated by both the ceRNA network and DNA methylation were identified. In the ceRNA network, the expression of RP11-545G3.1 lncRNA was validated in nasopharyngeal carcinoma tissues and cells by RT-qPCR. After a knockdown of RP11-545G3.1, the viability, migration, and invasion of CNE-2 and NP69 cells was assessed by CCK-8, wound healing and Transwell assays. Results This study identified abnormally expressed mRNAs, miRNAs and lncRNAs in nasopharyngeal carcinoma tissues. A ceRNA network was constructed, which contained three lncRNAs, 15 miRNAs and 129 mRNAs. Among the nodes in the PPI network based on the mRNAs in the ceRNA network, HMGA1 was assessed in relation to the overall and disease-free survival of nasopharyngeal carcinoma. We screened two up-regulated genes regulated by the ceRNA network and hypomethylation and 26 down-regulated genes regulated by the ceRNA network and hypermethylation. RP11-545G3.1 was highly expressed in the nasopharyngeal carcinoma tissues and cells. Moreover, the knockdown of RP11-545G3.1 reduced the viability, migration, and invasion of CNE-2 and NP69 cells. Conclusion Our findings uncovered the epigenetic regulation in nasopharyngeal carcinoma and identified the implications of RP11-545G3.1 on the progression of nasopharyngeal carcinoma. Supplementary Information The online version contains supplementary material available at 10.1186/s12920-023-01653-1.


Introduction
Nasopharyngeal carcinoma is a malignant tumor that occurs in the nasopharyngeal mucosa epithelium and small salivary gland, which is primarily a poorly differentiated squamous cell carcinoma [1].Nasopharyngeal carcinoma is highly malignant and prone to metastasis, as well as recurrence.This malignancy is especially prevalent in the regions of east and southeast Asia.Unfortunately, approximately 20% of patients with nasopharyngeal carcinoma continue to have local recurrence and distant metastasis after treatment, which is mainly due to the radioresistance of nasopharyngeal carcinoma [2] [3].Because nasopharyngeal carcinoma is highly malignant and progresses rapidly, most patients are diagnosed as the middle and advanced stages.The 5-year survival rate of patients with advanced nasopharyngeal carcinoma is less than 60% [3].Therefore, it is significant to investigate further treatment strategies to provide novel options for treating nasopharyngeal carcinoma.
The analysis of molecular characteristics using highthroughput sequencing has greatly improved our understanding of the pathogenesis of nasopharyngeal carcinoma.Non-coding RNAs (e.g., long-chain noncoding RNAs [lncRNAs] and microRNAs [miRNAs]) are involved in regulating physiological processes and pathogenesis of various diseases [4].Dysregulated competing endogenous RNA (ceRNA) is related to the occurrence and progression of cancers, including nasopharyngeal carcinoma [5].CeRNA acts as a sponge of miRNA via an miRNA response element, thereby forming an intricate ceRNA network [6].Nevertheless, the molecular mechanisms of ceRNA have not been fully elucidated in nasopharyngeal carcinoma.The pathogenesis of nasopharyngeal carcinoma is a biological process which includes both genetic and environmental factors.Moreover, the accumulation of genomic and epigenetic changes contributes to the occurrence of nasopharyngeal carcinoma [7].DNA methylation is the most common modification in epigenetics, and its abnormal changes represent the main factors that lead to the development of various tumors, including nasopharyngeal carcinoma [8].Abnormal methylation can affect the function of key genes by changing their expression and thus affecting various processes during tumor development [9].Genespecific changes in methylation that occur during the early stages of cancer can assist with early cancer detection and prevention [10].Thus, it is important to comprehensively uncover abnormally expressed genes that are regulated by DNA methylation.
This study aimed to uncover the ceRNA network and DNA methylation associated with gene expression in nasopharyngeal carcinoma and validate the function of RP11-545G3.1 lncRNA on the progression of nasopharyngeal carcinoma.

Nasopharyngeal carcinoma and preprocessing datasets
The mRNA, miRNA and lncRNA expression profiles and DNA methylation profiles of nasopharyngeal carcinoma samples were retrieved from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/gds/)database.The mRNA expression profiles of 12 nasopharyngeal carcinoma and four normal nasopharyngeal samples on the GPL570 platform were extracted from the GSE64634 dataset [11].The miRNA expression profiles of 312 paraffin-embedded nasopharyngeal carcinoma and 18 normal nasopharyngeal tissues were obtained from the GSE32960 dataset on the GPL14722 platform [12].Furthermore, the lncRNA expression profiles of four nasopharyngeal carcinoma and para-carcinoma tissues were retrieved from the GSE95166 dataset using the GPL15314 platform [13].In addition, this study downloaded the lncRNA expression profiles of three paired nasopharyngeal carcinoma and normal nasopharynx specimens on the GPL16956 platform [5].The expression matrix data were normalized by quartile utilizing the normalizeBetweenArrays option in the limma package [14].Moreover, the DNA methylation data were normalized using the wateRmelon package [15].Outliers were detected and removed by employing a principal component analysis (PCA) using the psych package.The correlation between samples at the gene expression level was analyzed and visualized as heat maps.R v3.6.0 were applied in this study with R packages including limma v3.1.6,Cytoscape v3.7, STRING v11, survival v2.43 were applied in this study [14,16,17].

Differential analysis
Differentially expressed mRNAs, miRNAs and lncRNAs, as well as differentially methylated genes (DMGs) were screened between nasopharyngeal carcinoma and normal nasopharyngeal tissue samples by applying the limma package (version 3.5.1)[18](S Fig. 1).The screening threshold was set as |fold change (FC)| ≥ 1.5 and adjusted P < 0.05.The results were visualized into a volcano map and clustering heatmap using the heatmap.2package [19].Quality control (QC) of raw sequences ensures the availability of RNA-seq data and reduces biases introduced by sequencing depth, read distribution, and coverage uniformity, etc. Supplementary plots on the possible batch effects and quality control of each of the analyzed datasets should be included.(S Figs. 5 and 6).

Function enrichment analysis
Based on the Gene Ontology (GO) database [20] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway database [21], the functional enrichment of differentially expressed genes (DEGs) or DMGs was analyzed.A Fisher's exact test was used to identify which group of genes exhibited the highest association with specific functional items.The P-value and false discovery rate (FDR) corresponding to each item was calculated [22].Items with FDR < 0.05 were considered to be significantly enriched.

Patients and specimens
This study included 21 patients with nasopharyngeal carcinoma admitted to the Affiliated Hospital of Youjiang Medical College for Nationalities from January 2018 to January 2019.Samples of nasopharyngeal carcinoma and matched normal tissues were collected from each patient.All patients underwent a pathological biopsy, and did not receive radiotherapy or chemotherapy prior to the biopsy.The biopsy combined with clinical manifestations and CT examination of the nasopharynx were used for a diagnosis of nasopharyngeal squamous cell carcinoma.The following exclusion criteria were used: (1) age < 18 or > 75 years; (2) combined with other primary tumors or metastatic cancer; (3) severe infectious disease and immune dysfunction; (4) pregnant and lactating women; and (5) mental disorders.All specimens were stored in liquid nitrogen.This study was approved by the Ethics Committee of Affiliated Hospital of Youjiang Medical College for Nationalities, Guangxi province, China, and all patients provided informed consent.

Transfection
CNE-2 and C666-1 cells were seeded into a 6-well plate (2 × 10 5 cells/well).After culturing the cells in antibioticfree cell culture medium for 24 h, the cell fusion rate reached 80%~85%.A Lipofectamine™2000 transfection kit (Invitrogen; USA) was performed in accordance with the manufacturer's instructions, and shRNAs against RP11-545G3.1 (GenePharma, Shanghai, China) and shnegative control were transfected into the cells.Fresh medium was applied 24 h after transfection and the cells were cultured for an additional 48 h.RT-qPCR was employed to detect the silencing effect of RP11-545G3.1.

Cell counting kit-8 (CCK-8)
The cells were seeded into a 96-well culture plate, and the cell density was adjusted to 5 × 10 3 cells/well.Then, 10 µL CCK-8 solution (Beyotime, Shanghai, China) was added to each well and the cells were incubated for 4 h at 37 °C before removing the supernatant.The samples were then incubated in 150 µL of a dimethyl sulfoxide solution for 10 min at 37 °C.The absorbance value (OD value) at a wavelength of 450 nm of the cells in each well was determined using a microplate reader (GENESYS 10 S UV-Vis, Thermo Fisher Scientific, UK).

Western blot
The total protein from the tumor samples or cells was extracted using a cell protein extraction kit (Thermo, USA).The bicinchoninic acid (Zhongshan Jinqiao Biotechnology Co., Ltd., Beijing, China) method was used to determine the total protein concentration and to perform a quantitative detection of the samples.A sodium lauryl sulfate-polyacrylamide gel was prepared, and 50 µg of protein sample was added to each well.After subjecting the samples to 75 V constant voltage electrophoresis for 30 min, 100 V constant voltage electrophoresis was performed for 2 h.The wet transfer method was applied to electrically transfer the protein to a polyvinylidene fluoride membrane.After blocking the membrane with 5% bovine serum albumin for 60 min, the corresponding primary antibody was added and incubated overnight at 4 °C.Next, a secondary antibody (goat anti-rabbit IgG H&L (HRP) preabsorbed, 1/1000; ab7090; Abcam, Cambridge, MA, USA) was incubated with the membrane for 1 h at room temperature.Using an Odyssey two-color infrared laser imaging system, the membrane was scanned.β-actin was used as an internal control to analyze the level of target protein expression.The primary antibodies included CARM1 (1/5,000; ab245466; Abcam), CCND2 (1/500; ab230883; Abcam), P21 (1/500; ab227443; Abcam) and β-actin (1/5,000; ab179467; Abcam).

Wound healing assay
The cells were seeded in a 12-well plate (4 × 10 4 /well).After 48 h, a 10 µL pipette tip was used to create a mark on the horizontal line perpendicular to the back of the culture plate.Serum-free culture medium was added and the cells were continued culturing for 24 h.At 0 h and 24 h, three fields of view were randomly selected for each well under a fluorescence microscope.

Transwell assay
Matrigel (BD Biosciences; San Jose, CA, USA) was spread evenly over the bottom of the upper Transwell chamber (Corning; NY, USA).After the Matrigel was completely solidified, the remaining liquid in the upper chamber was removed, and 50 µL of serum-free DMEM medium was added.The shRNA-treated cells were resuspended into 200 µL serum-free medium and inoculated into the upper layer of the chamber (4 × 10 4 cells/well).At the same time, 500 µL cell culture medium containing 20% FBS was added to the lower Transwell chamber.After the Transwell chamber was incubated in a cell incubator for 24 h, the medium was removed and the cells in the upper chamber were removed with a cotton swab and then fixed in 10% methanol solution for 20 min.After rinsing the chambers three times with PBS, the cells were stained with crystal violet solution for 15 min.Three fields were randomly selected for high-power imaging.

Statistical analysis
Statistical analysis was performed using R language and GraphPad Prism software (version 8.0, San Diego, CA, USA, www.graphpad.com).The measurement data were displayed as the mean ± standard deviation.Kaplan-Meier curves for the overall survival and disease-free survival of the genes with degree ≥ 2 in the PPI network were carried out for the TCGA-HNSC dataset utilizing the "Survival" package, followed by a log-rank test.A Student's t-test was applied to compare the differences between two groups.Multiple comparisons were analyzed by a oneway ANOVA, followed by Tukey's multiple comparisons test.P < 0.05 was considered to be statistically significant.

Analysis of differentially expressed mRNAs in nasopharyngeal carcinoma
The mRNA expression profiles of 12 nasopharyngeal carcinoma and 4 normal samples were extracted from the GSE64634 dataset.The PCA results revealed that nasopharyngeal carcinoma samples differed from normal samples (Fig. 1A).There was a high correlation between the samples at the mRNA level (Fig. 1B).Compared with normal samples, there were 1,078 up-regulated and 1,113 down-regulated mRNAs in the nasopharyngeal carcinoma samples (Fig. 1C−E).The results of the GO enrichment analysis showed that these up-regulated mRNAs were primarily related to cellular proliferation (e.g., cell cycle, mitotic cell cycle, DNA replication, nuclear DNA replication, DNA-dependent DNA replication, cell division, chromosome segregation, mitotic nuclear division and mitotic cell cycle phase transition) (Fig. 1F).The down-regulated mRNAs were primarily correlated with cell motility (Fig. 1G).Similarly, the up-regulated mRNAs were mainly enriched in cell proliferation-related pathways (Fig. 1H) whereas the down-regulated mRNAs were enriched in relation to metabolic pathways (Fig. 1I).

Analysis of abnormally expressed miRNAs in nasopharyngeal carcinoma
We obtained the miRNA expression profiles of 312 nasopharyngeal carcinoma and 18 normal specimens from the GSE32960 dataset.According to the PCA results, no outlier samples were removed (S Fig. 2A).A high correlation between the samples was detected, as shown in S Fig. 2B.There were 91 up-and 94 down-regulated miRNAs in the nasopharyngeal carcinoma specimens in comparison to that of the normal specimens (S Fig. 2C−E).Based on intersections of the target mRNAs of differentially expressed miRNAs from miRTarBase, TargetScan, miRDB, miRanda and miRMap databases, 220 mRNAs were identified to be targeted by these miR-NAs (S Fig. 2F).

Establishment of a PPI network based on mRNAs in the ceRNA network in nasopharyngeal carcinoma
We established a PPI network containing 38 nodes in nasopharyngeal carcinoma based on the mRNAs in the ceRNA network using the STRING database (Fig. 4A).There was a total of 29 up-regulated and nine downregulated mRNAs in the nasopharyngeal carcinoma samples compared to that of the normal specimens.Table 1 lists the degree of nodes in the PPI network, among which, TP53 had the highest degree (n = 13).The functions of the nodes in the network were further analyzed.Figure 4B shows that the cell cycle and metabolic processes were significantly enriched.These genes were found to be involved in key cellular components like the chromosomes, cytoplasm, and cytoskeleton (Fig. 4C).In addition, the mRNAs in the PPI network had the complex molecular functions of adenyl nucleotide, ATP, and enzyme binding (Fig. 4D).The KEGG enrichment results demonstrated that various cancer-associated pathways were distinctly enriched by these genes (Fig. 4E).

HMGA1 is closely related to the prognosis of nasopharyngeal carcinoma
The prognostic values of the mRNAs with a degree ≥ 2 in the PPI network were assessed in depth.As a result, only HMGA1 was found to be related to the prognosis of nasopharyngeal carcinoma.Patients with high HMGA1 expression were associated with undesirable overall survival (HR = 1.434;P = 8.63e − 03 ; Fig. 5A) and disease-free survival (HR = 1.466;P = 2.43e − 02 ; Fig. 5B) compared to those with low expression.

Analysis of differentially methylated genes in nasopharyngeal carcinoma
The GSE62336 dataset was employed to analyze DMGs in nasopharyngeal carcinoma.At the methylation level, distinct differences between nasopharyngeal carcinoma and normal samples were observed (S Fig. 3A).A high correlation between the samples were found, as shown in S Fig. 3B.A total of 12,383 hypermethylated and 853 hypomethylated sites were detected in the nasopharyngeal carcinoma specimens compared to that in the normal samples (S Fig. 3C−E).To obtain the abnormally expressed genes that were affected by DNA methylation, we intersected DEGs with DMGs.As a result, 18 highly expressed and hypomethylated genes (S Fig. 3F), as well Fig. 5 The relationship between HMGA1 expression and the prognosis of nasopharyngeal carcinoma.Kaplan-Meier curves of (A) overall survival and (B) disease-free survival between high and low expression HMGA1 groups as 296 lowly expressed and hypermethylated genes (S Fig. 3G) were found in nasopharyngeal carcinoma.S Fig. 3H shows the PPI network based on lowly expressed and hypermethylated genes.There were 84 nodes in the network, among which, DNAI2 and CCDC39 had the highest degree (n = 6).We also found that highly expressed and hypomethylated genes were distinctly enriched in metabolic and catabolic processes (S Fig. 3I), as well as pyrimidine metabolism and small cell lung cancer pathways (S Fig. 3J).Moreover, lowly expressed and hypermethylated genes were distinctly enriched in the processes governing cilium movement and assembly (S Fig. 3K), as well as various metabolic and cancer-related pathways (S Fig. 3L).

Knockdown of RP11-545G3.1 reduces the migration and invasion of nasopharyngeal carcinoma cells
The results of the wound healing assay revealed that a knockdown of RP11-545G3.1 distinctly reduced the migrated capacity of CNE-2 and NP69 cells (Fig. 7A,  B).Furthermore, we found that the invasive capacities of CNE-2 and NP69 cells were weakened by a RP11-545G3.1 knockdown (Fig. 7C and D).Taken together, these findings indicated that silencing RP11-545G3.1 may decrease the migration and invasion of nasopharyngeal carcinoma cells.

Discussion
Nasopharyngeal carcinoma is the most common malignant tumor of the head and neck particularly in Southern China [3].Because the anatomical components are hidden and difficult to detect early, most patients have already progressed to the middle and late stages of the disease at the time of diagnosis.And the prognosis of patients with advanced stage is not ideal due to the high degree of malignancy associated with nasopharyngeal carcinoma.Thus, clarification of the molecular mechanism associated with nasopharyngeal carcinoma progression will aid in the discovery of novel therapeutic targets and improve patient prognosis.With the continuous deepening of gene-targeted therapy research, the identification of effective target genes is of great significance for improving the prognosis of patients with nasopharyngeal carcinoma.Furthermore, the specific genetic pathways involved in nasopharyngeal carcinoma remain elusive [31].Microarray technology has assisted with probing into the various changes in genes in nasopharyngeal carcinoma, and has been confirmed as convenient method that can be used to screen for promising markers in nasopharyngeal carcinoma [13].In this study, 1,078 up-regulated and 1,113 down-regulated mRNAs were identified in nasopharyngeal carcinoma tissues.These up-regulated mRNAs were primarily related to cellular proliferation whereas the down-regulated mRNAs were mainly correlated with cellular motility and metabolic pathways.Together, these findings indicated that these mRNAs may participate in the progression of nasopharyngeal carcinoma progression.As a ceRNA, lncRNA serves as an miRNA sponge, and may regulate the expression of tumor-related mRNAs.In the present study, a total of 91 up-and 94 down-regulated miRNAs were identified in nasopharyngeal carcinoma specimens.Furthermore, we screened five up-regulated lncRNAs (LOC100144603, RP11-545G3.1, SCARNA22, TAZ, and TIMM8B) and six down-regulated lncRNAs (AC003986.7,AK055386, BC013821, DQ786304, LOC647979 and RP11-179H18.2) in nasopharyngeal carcinoma specimens.It was previously found that SCARNA22 was overexpressed in breast cancer brain metastases [32].However, the functions of other lncRNAs in cancer are poorly understood.We subsequently established a ceRNA network containing 3 lncRNAs (RP11-179H18.2,RP11-545G3.1, and TAZ), 15 miRNAs and 129 mRNAs.In previous studies, Zou et al. constructed a lncRNA-mediated ceRNA network for nasopharyngeal carcinoma using weighted correlation network analyses [33].Moreover, Xu et al. identified the key genes related to a ceRNA network in nasopharyngeal carcinoma by applying bioinformatics analyses [13].
Based on the mRNAs in the ceRNA network, we established a PPI network containing 29 up-and nine downregulated mRNAs in nasopharyngeal carcinoma.Among these mRNAs, only HMGA1 was found to be related to the overall survival and disease-free survival of nasopharyngeal carcinoma patients which will be further validated in a larger population.The carcinogenesis of HMGA1 has been confirmed in numerous malignant cancers [34].For instance, HMGA1 may promote breast cancer aggressiveness and angiogenesis [35].HMGA1 induces tumor growth via regulation of cell cycle as well as facilitates migrated and invasive capacities in cervical cancer [36].Nevertheless, there is no research report on the role of HMGA1 in nasopharyngeal carcinoma.
The advantage of this study is the comprehensive analysis from miRNA-mRNA and lncRNA-miRNA, the ceRNA networks related to nasopharyngeal carcinoma are gradually constructed.On the one hand, the interactions between RNA can be understood at the overall transcription level of cells.On the other hand, it can help to predict the target RNA and protein which may be related to the pathogenesis of nasopharyngeal carcinoma.Other highlights of this study include the first exposure of LncRNA RP11-545GG3.1,which is up-regulated in nasopharyngeal carcinoma and its knockdown decreases nasopharyngeal carcinoma cell viability.However, this study also has certain limitations: For the data used in this study are survival data of head and neck squamous cell carcinoma (and nasopharyngeal carcinoma belongs to head and neck squamous cell carcinoma), which may not represent the most real situation, but has certain clinical reference value.Another limitation of this paper is that batch effects or confounding factors within the data set were not adjusted for factors such as age in the EWAS study.In the future, we will expand the sample size of nasopharyngeal carcinoma and try to remove possible confounding factors for large-scale verification.

Conclusion
Collectively, this study revealed that the ceRNA network and abnormal DNA methylation were associated with gene expression in nasopharyngeal carcinoma, which deepened our understanding of the pathogenesis of this disease.Furthermore, we identified the novel lncRNA RP11-545G3.1 as a promising therapeutic target against nasopharyngeal carcinoma.

Fig. 1
Fig. 1 Analysis of differentially expressed mRNAs for nasopharyngeal carcinoma in the GSE64634 dataset.(A) PCA results of 12 nasopharyngeal carcinoma (green) and 4 normal samples (blue).(B) Heat map for the correlations between samples.The intensity of the color is proportional to the correlation coefficient.(C−E) Scatter plots, volcano diagram, and heat map of the differentially expressed mRNAs between nasopharyngeal carcinoma and normal samples.Blue: down-regulation; red: up-regulation.(F, G) The top 10 GO enrichment results of (F) up-and (G) down-regulated mRNAs.GO database: BP = biological process, CC = cellular component and molecular function, MT = molecular function.(H, I) The top 10 KEGG pathways enriched by (H) upand (I) down-regulated mRNAs.

Fig. 2
Fig. 2 Screening of dysregulated lncRNAs in nasopharyngeal carcinoma in the GSE95166 and GSE126683 datasets.(A, B) PCA results of nasopharyngeal carcinoma and normal specimens.(C, D) Heat map for the correlations between samples.(E−J) Scatter plots, volcano diagram, and heat map for dysregulated lncRNAs between nasopharyngeal carcinoma and normal specimens.Blue: down-regulation; red: up-regulation.(K, L) Venn diagrams for the common (K) up-and (L) down-regulated lncRNAs in the GSE95166 and GSE126683 datasets

Fig. 3
Fig. 3 Establishment of a nasopharyngeal carcinoma-related ceRNA network.(A) A ceRNA network containing lncRNAs (V type), miRNAs (diamond) and mRNAs (circle).The darker the color, the larger the |log 2 FC|.(B) KEGG enrichment results of the mRNAs in the ceRNA network.(C−E) GO enrichment results of the mRNAs in the ceRNA network, including (C) biological process, (D) cellular component, and (E) molecular function

Fig. 4
Fig. 4 Establishment of a PPI network based on the mRNAs in the ceRNA network for nasopharyngeal carcinoma.(A) A PPI network in nasopharyngeal carcinoma.Red: up-regulation; blue: down-regulation.The darker the color, the larger the |log2FC|.(B−D) GO enrichment results of mRNAs in the PPI network containing (B) biological processes, (C) cellular components, and (D) molecular functions.(E) KEGG enrichment pathways of mRNAs in the PPI network

Table 1
The degree of nodes in the PPI network