Research Paper Volume 13, Issue 18 pp 22176—22187

Key genes affecting the progression of nasopharyngeal carcinoma identified by RNA-sequencing and bioinformatic analysis

Yihong Wang1, *, , Manyi Li1, *, , Yan Guo1, *, , Haiping Huang1, , Xuelin Dong1, , Yangguang Sun1, , Jisheng Liu1, ,

  • 1 Department of Otolaryngology, The First Affiliated Hospital of Soochow University, Soochow 215000, Jiangsu Province, China
* Co-first authors

Received: November 13, 2020       Accepted: August 2, 2021       Published: September 20, 2021      

https://doi.org/10.18632/aging.203521
How to Cite

Copyright: © 2021 Wang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Background: The present work was conducted to screen the potential biomarkers affecting nasopharyngeal carcinoma (NPC) progression through RNA-sequencing (RNA-seq), bioinformatic analysis and functional experiments.

Materials and Methods: Six normal samples and five NPC clinical samples were collected for RNA-seq analysis. The expression levels in both groups were determined through student’s t-test. We identified genes of P < 0.01 as the differentially expressed genes (DEGs). In addition, gene set enrichment analysis (GSEA) was conducted. Afterwards, STRING V10 database was employed to extract protein interactions among the DEGs. Later, we established a protein-protein interaction (PPI) network, and used the Cytoscape software for network visualization. qRT-PCR was conducted to verify hub genes from clinical samples. Then, the function of CXCL10 in cell proliferation, apoptosis, invasion and migration was evaluated.

Results: A total of 2024 DEGs were identified, among which, 1449 were down-regulated and 575 were up-regulated. The PPI was constructed, and the hub genes including Insulin Like Growth Factor 1 (IGF1), C-X-C Motif Chemokine Ligand 10 (CXCL10), Interleukin 13 (IL13), Intercellular Adhesion Molecule 1 (ICAM1), G Protein Subunit Gamma Transducin 1 (GNGT1), Matrix Metallopeptidase 1 (MMP1), Neurexin 1 (NRXN1) and Matrix Metallopeptidase 3 (MMP3) were obtained. The expression levels of CXCL10, IGF1, MMP3, MMP1, ICAM1, and IL-13 were significantly up-regulated in tumor tissues. High expression levels of CXCL10, MMP3 and ICAM1 predicted poor prognosis of NPC patients. CXCL10 silencing suppressed NPC cell proliferation and migration.

Conclusions: CXCL10 may serve as a potential key gene affecting NPC genesis and progression.

Introduction

Nasopharyngeal carcinoma (NPC) originates in the upper and lateral nasopharyngeal cavity wall, which shows a high incidence in China and has the highest incidence among otorhinolaryngological cancers. NPC possesses a strong erosive activity and is prone to invasive expansion into surrounding organs and tissues, with frequent occurrence of regional lymph node metastasis (LNM). Due to the lack of specific early clinical manifestations, NPC is easily misdiagnosed as ear, nose and throat (ENT) diseases and is not paid attention to by patients. Therefore, most patients are diagnosed at stages III and IV. Because of the high recurrence rate and distant metastasis (DM) rate, locally advanced patients are mainly treated with radiotherapy supplemented with chemotherapy, biotherapy and other comprehensive treatment to improve the therapeutic effect.

Over the past few decades, RNA-sequencing (RNA-seq) has been applied in the diagnosis and molecular targeted therapy of numerous diseases. It is important to elucidate the functions of novel genes in diagnosing, treating and predicting the prognosis for NPC in clinic. A number of biomarkers in NPC have been reported in recent years. For instance, He et al. established a signature based on the long non-coding RNAs (lncRNAs) in serum to diagnose NPC. According to their results, AL359062, MALAT1 and AFAP1-AS1 were the candidate new serum biomarkers [1]. Luo et al. discovered that FoxM1 was over-expressed in NPC tissues, and FoxM1 induced the progression and cancer stem cell (CSC) features in NPC [2]. According to Fan et al., the CORO1C level was related to NPC metastasis; meanwhile, CORO1C overexpression enhanced the migration and invasion of NPC cells [3].

In the present study, NPC clinical samples and normal samples were collected for RNA-seq. In addition, we selected differentially expressed genes (DEGs), and conducted gene set enrichment analysis (GSEA). Through constructing the PPI network, the hub genes CXCL10, GNGT1, IGF1, MMP3, MMP1, ICAM1, IL13, and NRXN1 were obtained. Later, CXCL10 levels within NPC tissues were analyzed, and the function of CXCL10 silencing in the proliferation, apoptosis, invasion and migration of cells was determined.

Materials and Methods

Ethical statement

The Ethics Committee of The First Affiliated Hospital of Suzhou University approved our study protocol. Each participant provided the informed consent for participation.

Patient samples and the clinicopathological features

We obtained samples from five NPC cases undergoing surgical treatment at the Department of Otolaryngology, The First Affiliated Hospital of Soochow University (Soochow, China) from May, 2016 to December, 2019 for RNA-seq. Additionally, six healthy controls were collected for RNA-seq, whereas additional 61 NPC patients and 61 healthy controls were collected for clinical verification. The clinicopathological characteristics of patients were recorded. The normal tissues were resected, and cancer tissues were carefully obtained from the middle of NPC. Each collected sample was frozen with liquid nitrogen and preserved under –80°C. At last, we acquired five NPC and six non-carcinoma samples for RNA-seq.

RNA preparation

We utilized Trizol (Invitrogen, Carlsbad, CA, USA) for extracting the total RNA from normal and NPC samples, DEPC-treated water was utilized to re-dissolve total RNA, whereas NanoVue Plus spectrophotometry (GE Healthcare, Fairfield, CT, USA) was carried out to quantify total RNA. Agarose gel electrophoresis was conducted to evaluate the integrity of RNA, and gDNA Eraser (Takara, Tokyo, Japan) was also utilized for eliminating DNA contamination in accordance with specific protocols.

RNA sequencing (RNA-Seq)

Trizol was used to grind the normal and NSP samples, whereas the RNeasy Kit (Qiagen) was employed to isolate total RNA in accordance with specific protocols. Besides, DNaseI was added to prevent tissues from mixing with genomic DNA. In addition, Agilent BioAnalyzer 2100 (Agilent, Santa Clara, CA, USA) was utilized to detect RNA purity. Moreover, we adopted the TruSeq RNA Access Library Prep Kit® (Illumina, CA, USA) to prepare transcriptomic sequencing libraries. Further, we utilized the Illumina HiSeq 2000 instrument to perform 100 bp paired-end sequencing, as supported by Beijing Novogene Biological Information Technology Co., Ltd., (Beijing, China).

Functional enrichment of DEGs

We utilized the clueGO plug-in of Cytoscape to carry out functional enrichment. ClueGO can be used for functional enrichment and classification of functional items with significant enrichment, including Gene Ontology (GO) terms and KEGG pathways [4]. During calculation, we determined the kappa-coefficient to reflect the functional associations of terms or pathways according to gene overlapping between different GO terms or pathways. In this study, the default value of kappa was 0.4. Entries with similar functions were presented by identical color. P < 0.05 was set as the significant enrichment threshold.

PPI network construction

PPI can be obtained from the STRING database based on different data, including fusion, coexpression, text-mining and cooccurrence. This database comprehensively grades every PPI relation pair, with the score ranging from 0 to 1 [5]. A greater score indicates that the PPI relation is more creditable. Generally, 0.4 is used as the threshold of pooled score. In this study, we obtained PPI among DEGs from the STRING V10 database, established the PPI network and visualized it using Cytoscape.

Immunohistochemistry (IHC)

After deparaffinization, fixed sections were rehydrated with gradient xylene and ethanol, and incubated with 0.3% H2O2 solution under ambient temperature for 30 min. Sodium citrate solution (pH 6.0) was used to retrieve the antigen. Thereafter, each section was blocked using the BSA solution for reducing the nonspecific binding. Later, each section was subjected to FGF5 antibody staining, and incubation using secondary antibody (Envision, Gene Technology, Shanghai, China) for visualization. Later, DAB was employed for color developing of each slide, whereas hematoxylin was used for counterstaining.

Cell culture and transfection

The CNE-2 and 5–8F NPC cells were provided by Central Laboratory of Soochow University (Suzhou, China) and cultivated within the RPMI-1640 medium (Gibco, Gaithersburg, MD, USA, A1049101) that contained 10% fetal bovine serum (FBS; Gibco, 10099141) under 37°C and 5% CO2 conditions.

Cell lines and cell culture

CNE-2 and 5–8F cells were obtained from the cell bank at the Chinese Academy of Sciences (Shanghai, China). Afterwards, all cells were cultivated within the RPMI-1640 medium that contained 10% FBS (Gibco, 10099141) under the humidified, 37°C and 5% CO2 conditions.

Cell transfection

The siCXCL10-1 (5′-ACTGCCATTCTGATTTGCTGCCTTA-3′), siCXCL10-1 (5′-GCTGCCTTATCTTTCTGACTCTAAG-3′) and negative control (siNC, 5′-CGGTATGCCATGCCTACGCTATCGAAC-3′) were provided by GenePharma (Shanghai, China). After culture within the antibiotics-free complete medium for 24 h pre-transfection, cells were rinsed by PBS, followed by temporary transfection using siNC or siFGF5 (50 nmol/L) using Lipofectamine 2000 (Invitrogen, USA) in accordance with specific protocols. At 24 h post-transfection, we conducted qRT-PCR for validating the CXCL10 level.

RNA Isolation and qRT-PCR

We utilized Trizol reagent (Takara, Japan) for extracting total tissue RNA from primary tumor and non-carcinoma samples, respectively. Later, the extracted total RNA was used to prepare cDNA by reverse transcription using the PrimeScript RT-PCR kit (Takara, Japan) in accordance with specific protocols. RT-PCR was carried out by the ABI7500 (Applied Biosystem) thermal cycler with the standard SYBR Green PCR kit (Takara). In this study, the primers used were shown below: for CXCL10, 5′-AAGCAGTTAGCAAGGAAAGG-3′ and 5′-GTAGGGAAGTGATGGGAGAG-3′; for GNGT1, 5′-GAAATTGGTTATCGTGGGAT-3′ and 5′-TCACTTCTTTCTTGAGCTGGT-3′; for IGF1, 5′-AGGAGGCTGGAGATGTATTG-3′ and 5′-GTGTTCTTGTTGGTAGATGGG-3′; for MMP3, 5′-GGTCTCTTTCACTCAGCCA-3′ and 5′-GGGTCTCAGGGGAGTCA-3′; for MMP1, 5′-GCTACACGGATACCCCAA-3′ and 5′-CTCAGAAAGAGCAGCATCG-3′; for ICAM1, 5′-AAACACTAGGCCACGCATC-3′ and 5′-CCCACCACTTCCCCTCT-3′; for IL13, 5′-CAGCTCAGGCACACTTCTT-3′ and 5′-CTAGCAGCCACAGTCTTCC-3′; for NRXN1, 5′TCCTCGGGTTAAGAAATGG-3′ and 5′-CCTTGGATGCTTGTGAATG-3′; for CXCR3, 5′-ACAAGCACCAAAGCAGAGG-3′ and 5′-CTGGGCAGCAGCACTTAC-3′; and for GAPDH, 5′-CACCCACTCCTCCACCTTTG-3′ and 5′-CCACCACCCTGTTGCTGTAG-3′. The 2–ΔΔCt approach was adopted for calculating relative gene level, with GAPDH as the endogenous control.

Cell proliferation assay

In the present work, we applied MTT assay for assessing cell proliferation according to specific instructions (Invitrogen, M6494). Briefly, corresponding vectors were transfected into CNE-2 and 5–8F cells, and incubated for 0, 24, 48 and 72 h, separately. After 20 μL MTT solution was added to incubate cells for 4 h, we discarded supernatants. Afterwards, the products were dissolved through the addition of dimethyl sulfoxide (100 μL, DMSO). At last, absorbance (OD) value was detected by using the microplate reader (Bio-Rad, Richmond, CA, USA, 168-1000XC) at 490 nm.

Cell apoptosis assay

The propidium iodide (PI)/Annexin-V Cell Apoptosis Kit (Invitrogen, V13245) was utilized to analyze C666-1 and 5–8F apoptosis rates according to specific protocols. Briefly, vectors were transfected into C666-1 and 5–8F cells and later cultured for 48 h. Then, cells were stained with FITC-Annexin V and PI. At last, we conducted flow cytometric analysis (FACSCantoII, 338960; BD Biosciences, San Jose, CA, USA) for determining cell apoptosis rate.

Transwell invasion assay

The Matrigel-coated (BD, Franklin Lakes, NJ, USA) Transwell chambers (Costar, Manassas, VA, USA) that contained 8-um polycarbonate filters were used for the invasion assay. At 48 h post-transfection, cells cultured within the FBS-free RPMI-1640 (200 μl) were added into the Matrigel (1 mg/ml, 30 μl)-coated upper chamber. Meanwhile, FBS-containing RPMI-1640 was added into the lower chamber, which served as the chemo-attractant. After incubation for 24 h, we scraped cells on membrane surface, then those scraped cells were subjected to PBS rinsing, 100% methanol fixation, and Giemsa dye staining.

Scratched wound healing assay

Cells (1 × 106 cells/well) were inoculated within the 6-well plates in the wound healing assay. When cells reached about 90% confluence, we used the pipette tip (1 ml) to made a wound on the surface of the cell monolayer. Microscopy was conducted to obtain images at 0 and 48 h, respectively, for assessing the wound closure rate. The ImageJ software (NIH, Bethesda, MD, USA) was applied in measuring inter-edge distance. To be specific, the wound closure rate was determined by the following formula: (scratch width at 0 h- scratch width at 24 h/48 h)/scratch width at 0 h × 100%.

Statistical analysis

We used the Fisher’s exact test to assess the clinicopathological features, such as sex, age, tumor stage, tumor size, and differentiation. Animal and functional experiments in vitro were carried out to assess the expression level by student’s t-test. SPSS18.0 (SPSS, Inc., Chicago, IL, USA) was utilized in all statistical analyses (two-sided). P < 0.05 was set as the significance level.

Results

Raw data preprocessing and DEGs screening

Gene expression was uniformly distributed within samples after preprocessing and the processed data were utilized in later analyses (Figure 1A). We classified 11 samples into 2 groups, including Tumor group (T, n = 5) and Normal group (N, n = 6). ANOVA was performed for analyzing DEGs in both groups, yielding altogether 2024 DEGs (Figure 1B). Among the DEGs, 1449 were down-regulated while 575 were up-regulated. The cluster analysis of DEGs can be observed from Figure 1C. The 10 respective most significantly up-regulated and down-regulated DEGs are presented in Table 1. Of them, the top 10 up-regulated DEGs were RP11-810K23.10, CHI3L1, KREMEN2, NRXN1, CHIT1, IBSP, CLSTN2, MMP12, WNT2, and CNTNAP2. The 10 most significantly down-regulated DEGs included RP11-706O15.3, AZU1, RP11-363J20.2, B3GALT5, ABO, GJB6, RP11-501J20.5, ATP1A2, RP11-646E20.6, and CYP2G1P.

Raw data preprocessing and DEGs screening. (A) The distribution of gene expression levels in the samples was relatively uniform. (B) A total of 2024 DEGs were obtained between Normal group (N, n = 6) and Tumor group (T, n = 5). 1449 DEGs were down-regulated, and 575 DEGs were up-regulated. (C) Heatmap shows the DEGs.

Figure 1. Raw data preprocessing and DEGs screening. (A) The distribution of gene expression levels in the samples was relatively uniform. (B) A total of 2024 DEGs were obtained between Normal group (N, n = 6) and Tumor group (T, n = 5). 1449 DEGs were down-regulated, and 575 DEGs were up-regulated. (C) Heatmap shows the DEGs.

Table 1. The respective most significantly up-regulated and down-regulated DEGs.

GenesFold changeP value
Upregulated
RP11-810K23.103.213.76E-06
CHI3L13.082.08E-09
KREMEN23.022.72E-06
NRXN12.998.70E-06
CHIT12.868.97E-07
IBSP2.852.65E-05
CLSTN22.846.63E-07
MMP122.801.39E-07
WNT22.754.73E-05
CNTNAP22.742.07E-06
Downregulated
RP11-706O15.3–2.942.36E-06
AZU1–2.781.57E-05
RP11-363J20.2–2.637.32E-05
B3GALT5–2.611.17E-05
ABO–2.581.27E-04
GJB6–2.551.97E-04
RP11-501J20.5–2.532.56E-04
ATP1A2–2.502.74E-05
RP11-646E20.6–2.503.35E-04
CYP2G1P–2.472.80E-04

Functional enrichment of NPC related genes

This study also examined the roles of DEGs through GO and KEGG analyses, respectively. Typically, GO enrichment of DEGs directly reflects the enrichment of GO terms including biological process (BP), cell component (CC) and molecular function (MF). We selected 30 most significant GO terms, as shown in Figure 2A. According to Figure 2A, DEGs were enriched in extracellular space (CC), transmembrane transport (BP), and actin binding (MF). We obtained 30 significantly enriched GO terms, including plasma membrane, extracellular region, extracellular space, cornified envelope, integral component of plasma membrane. As observed from the figure of P-value, DEGs were enriched in pathways such as extracellular space, extracellular region, and plasma membrane pathway (Figure 2B). We then performed KEGG pathway enrichment analysis, which revealed that DEGs were related to pathways such as Pancreatic secretion (organismal systems), Alanine, aspartate and glutamate metabolism (metabolism), Breast cancer (human diseases), and Neuroactive ligand-receptor interaction signaling pathways (environmental information processing) (Figure 2C). Meanwhile, as suggested from the KEGG pathway scatter plot, DEGs were mostly related to Neuroactive ligand-receptor interaction, Cytokine receptor interaction, and protein digestion and absorption (Figure 2D).

Functional enrichment of NPC related DEGs. (A–B) DEGs were enriched in GO terms. (C–D) DEGs were enriched in KEGG pathways.

Figure 2. Functional enrichment of NPC related DEGs. (AB) DEGs were enriched in GO terms. (CD) DEGs were enriched in KEGG pathways.

PPI network of NPC related genes

Later, the NPC-related DEGs were extracted from the STRING V10 database to construct a PPI network. Figure 3 presents the established PPI. Typically, the 8 most significant nodes within the PPI network were Insulin Like Growth Factor 1 (IGF1, degree = 27), Interleukin 13 (IL13, degree = 20), Intercellular Adhesion Molecule 1 (ICAM1, degree = 18), G Protein Subunit Gamma Transducin 1 (GNGT1, degree = 17), C-X-C Motif Chemokine Ligand 10 (CXCL10, degree = 17), Matrix Metallopeptidase 1 (MMP1, degree = 17), Matrix Metallopeptidase 3 (MMP3, degree = 16), and Neurexin 1 (NRXN1, degree = 15). These hub nodes had the highest degree and edge values, indicating their irreplaceable roles in NPC. Therefore, further experiments should be performed to investigate the above genes for their clinical meaning.

Protein-protein international (PPI) network of nasopharyngeal carcinoma related genes.

Figure 3. Protein-protein international (PPI) network of nasopharyngeal carcinoma related genes.

Expression levels of the top 8 DEGs in NPC patients

We collected 61 pairs of NPC and normal samples to verify the expression levels of the top 8 DEGs. Levels of CXCL10, GNGT1, IGF1, MMP3, MMP1, ICAM1, IL13, and NRXN1 were measured through qRT-PCR assay, as presented in Figure 4A. As a result, CXCL10, IGF1, MMP3, MMP1, ICAM1, and IL-13 levels were markedly up-regulated within NPC samples relative to the non-carcinoma samples (P < 0.05). While, difference in the expression of GNGT1 and NRXN1 was not significant in NPC samples compared with non-carcinoma samples. Additionally, the survival analysis of CXCL10, IGF1, MMP3, MMP1, ICAM1, and IL-13 was performed. The results showed that high expression levels of CXCL10, MMP3 and ICAM1 predicted poor prognosis of NPC patients (Figure 4B). The clinical significance of MMP-3 and ICAM1 in NPC is already reported [6], so this study focused on CXCL10.

The expressions of DEGs in NPC patients. (A) The qRT-PCR assay was employed to detect the expressions of CXCL10, GNGT1, IGF1, MMP3, MMP1, ICAM1, IL13, and NRXN1. (B) The survival analysis of IGF1, IL13, ICAM1, CXCL10, GNGT1, MMP1, MMP3, and NRXN1 in NPC patients.

Figure 4. The expressions of DEGs in NPC patients. (A) The qRT-PCR assay was employed to detect the expressions of CXCL10, GNGT1, IGF1, MMP3, MMP1, ICAM1, IL13, and NRXN1. (B) The survival analysis of IGF1, IL13, ICAM1, CXCL10, GNGT1, MMP1, MMP3, and NRXN1 in NPC patients.

CXCL10 knockdown suppressed NPC cell proliferation and invasion

This study detected the CXCR3 (CXCL10 receptor) expression, and the result showed that the CXCR3 was over-expressed in NPC samples (Figure 5A). IHC assay was conducted for identifying CXCL10 levels within NPC samples and non-carcinoma samples. As a result, CXCL10 over-expression was detected within NPC samples (Figure 5A). As presented in Table 2, the high CXCL10 expression showed positive correlation with TNM stage, T stage and N stage. After siNC, siCXCL10-1 or siCXCL10-2 was transfected into 5–8F and CNE-2 cells, RT-PCR was performed to examine the transfection efficiency (Figure 5B). At last, we selected siCXCL10-2 for further experiments. Cell proliferation, apoptosis, invasion and migration were assessed. According to our results, down-regulation of CXCL10 remarkably suppressed NPC cell proliferation, migration and invasion, and promoted their apoptosis relative to cells transfected with siNC (P < 0.01, Figure 5C5F).

Knockdown of CXCL10 inhibited cell proliferation and invasion of NPC cells. (A) The expression of CXCR3 was examined by qRT-PCR assay. The expression of CXCL10 in NPC tumor tissue and normal tissue was examined by IHC analysis. (B) CNE-2 and 5–8F cells were transfected with siNC, siCXCL10-1 or siCXCL10-2, and RT-PCR was performed to examine the CXCL10 expression. (C) Cell proliferation of CNE-2 and 5–8F cells was identified by CCK8 assay. (D) Cell apoptosis of CNE-2 and 5–8F cells was evaluated by FCM. (E) Cell invasion was tested by wound healing assay. (F) Cell migration was detected by Transwell assay. **P

Figure 5. Knockdown of CXCL10 inhibited cell proliferation and invasion of NPC cells. (A) The expression of CXCR3 was examined by qRT-PCR assay. The expression of CXCL10 in NPC tumor tissue and normal tissue was examined by IHC analysis. (B) CNE-2 and 5–8F cells were transfected with siNC, siCXCL10-1 or siCXCL10-2, and RT-PCR was performed to examine the CXCL10 expression. (C) Cell proliferation of CNE-2 and 5–8F cells was identified by CCK8 assay. (D) Cell apoptosis of CNE-2 and 5–8F cells was evaluated by FCM. (E) Cell invasion was tested by wound healing assay. (F) Cell migration was detected by Transwell assay. **P < 0.01 vs siNC group.

Table 2. The relationship between CXCL10 expression and clinicopathologic characteristics of NPC patients.

CharacteristicsNumber of patientsCXCL10CXCL10P value
Low expression
(<median)
High expression
(≥median)
Number612932
Ages0.452
<50311417
≥50301515
Gender0.357
Female291514
Male321418
TNM stage0.048
I-III301812
IV311120
T stage0.015
T1-T2281810
T3-T4331122
N stage0.045
N0-N1321913
N2-N3291019

Discussion

In this work, we performed RNA-seq and bioinformatic analysis to screen the crucial biomarkers affecting the occurrence and development of NPC. The up-regulated and down-regulated DEGs were obtained at first, and then GSEA analysis was processed. Several pathways including extracellular space [7], extracellular region [8], plasma membrane [9], and integral component of plasma membrane [10] were closely associated with tumor genesis, growth and metastasis. Loom et al. investigated that extracellular space was an important compartment for malignant energetic catalysis and therapeutic targeting [7]. It is reported that extracellular region plays a crucial role in oncogenic events of tropomyosin receptor kinase (TRK) in several cancers [8]. In addition, plasma membrane and integral component of plasma membrane pathways are also proved to regulate cancer progression on breast cancer [11]. Some cancer-related KEGG pathways are identified, such as breast cancer, Cytokine receptor interaction [12], and protein digestion and absorption [13]. A lot of studies report that DEGs including mRNAs and non-coding RNAs in many kinds of cancers are enriched in KEGG Cytokine receptor interaction pathway [14, 15]. And, Han et al. identified abnormally methylated DEGs and pathways associated with NPC, and they found that DEGs were also enriched in protein digestion and absorption pathway [16]. Furthermore, we constructed the PPI network, and identified IGF1, IL13, ICAM1, CXCL10, GNGT1, MMP1, MMP3, and NRXN1 as the hub genes. Our clinical data verified that the expression levels of IGF1, IL13, ICAM1, CXCL10, MMP1, and MMP3 were up-regulated in NPC tissues. IGF1 is a polypeptide that is structurally similar to human pro-insulin, which is one of the pathogenetic factors resulting in obesity and other diseases [17]. IGF1 is also a high risk factor for several cancers, such as prostate cancer [18], breast cancer [19], lung cancer [20], gastric cancer [21], and colon cancer [22]. The high serum IGFBP-1-to-IGF-1 ratio is related to the adverse prognostic outcome for NPC [23]. Wang et al. reported that the loss of miR-206 reduced the radiosensitivity of NPC by targeting IGF-1 [24]. IL13 serves as an immunoregulatory cytokine, which is closely related to cancer genesis by impacting the tumor immune surveillance [25]. ICAM1, MMP1, and MMP3 are the well-known biomarkers participating in regulating cancer cell migration and invasion. ICAM1 is involved in tumor cell adhesion to vascular endothelium or neutrophils, and mediates the hematogenous and lymph node metastasis of malignant tumors [26, 27]. The high expression levels of MMP1 and MMP3 enhance cancer cell migration and invasion, and predict poor prognosis for NPC [28, 29]. CXCL10 is over-expressed in several kinds of cancers, and it is reported to regulate cancer progression, such as colorectal cancer (CRC) [30], cervical cancer [31], ovarian cancer (OC) [32], and breast cancer (BC) [33, 34]. However, the function of CXCL10 in NPC is still unknown. As revealed by our results, CXCL10 was over-expressed in NPC tissues, and the up-regulated CXCL10 expression showed positive correlation with TNM stage, T stage and N stage. The above findings displayed that CXCL10 possibly had oncogenic effect in NPC. CXCR3, as a CXCL10 receptor, was also up-regulated in NPC samples. Thereafter, we performed cellular functional assays to identify the role of CXCL10 in tumor cell activities. As a result, the down-regulated CXCL10 expression suppressed CNE-2 and 5–8F proliferation, migration and invasion, but enhanced their apoptosis. These data implied that CXCL10 might serve as a candidate biomarker for molecular diagnosis and targeted therapy of NPC.

Certain limitations should be noted in the present work. First of all, this study did not fit the overall survival time or the gene expression levels. Secondly, the present work involved few clinical indicators and sample size. At last, the detailed molecular mechanism by which CXCL10 affected NPC cell proliferation and invasion should be further studied. Our next study will emphasize the above issues. To sum up, CXCL10 is a candidate key gene affecting NPC genesis and progression. This work sheds new lights on the molecular diagnosis and targeted therapy of NPC.

Author Contributions

YW, ML, YG and JL contributed to the study design. YW, ML and HH were responsible for study implementation. XD, YS and JL analyzed the data. YW, ML, YG and JL wrote the manuscript. All authors have carefully read and approved the final manuscript.

Conflicts of Interest

The authors declare no conflicts of interest related to this study.

Funding

The present work was funded by the Jiangsu Science and Technology Plan Project (BE2019670) and Suzhou people's livelihood science and technology project (SYS2020115).

References

  • 1. He B, Zeng J, Chao W, Chen X, Huang Y, Deng K, Huang Z, Li J, Dai M, Chen S, Huang H, Dai S. Serum long non-coding RNAs MALAT1, AFAP1-AS1 and AL359062 as diagnostic and prognostic biomarkers for nasopharyngeal carcinoma. Oncotarget. 2017; 8:41166–77. https://doi.org/10.18632/oncotarget.17083 [PubMed]
  • 2. Luo W, Gao F, Li S, Liu L. FoxM1 Promotes Cell Proliferation, Invasion, and Stem Cell Properties in Nasopharyngeal Carcinoma. Front Oncol. 2018; 8:483. https://doi.org/10.3389/fonc.2018.00483 [PubMed]
  • 3. Fan L, Wei Y, Ding X, Li B. Coronin3 Promotes Nasopharyngeal Carcinoma Migration and Invasion by Induction of Epithelial-to-Mesenchymal Transition [Retraction]. Onco Targets Ther. 2020; 13:7355–56. https://doi.org/10.2147/OTT.S273162 [PubMed]
  • 4. Dalmer TRA, Clugston RD. Gene ontology enrichment analysis of congenital diaphragmatic hernia-associated genes. Pediatr Res. 2019; 85:13–19. https://doi.org/10.1038/s41390-018-0192-8 [PubMed]
  • 5. Cook HV, Doncheva NT, Szklarczyk D, von Mering C, Jensen LJ. Viruses.STRING: A Virus-Host Protein-Protein Interaction Database. Viruses. 2018; 10:519. https://doi.org/10.3390/v10100519 [PubMed]
  • 6. Busson P, Zhang Q, Guillon JM, Gregory CD, Young LS, Clausse B, Lipinski M, Rickinson AB, Tursz T. Elevated expression of ICAM1 (CD54) and minimal expression of LFA3 (CD58) in Epstein-Barr-virus-positive nasopharyngeal carcinoma cells. Int J Cancer. 1992; 50:863–67. https://doi.org/10.1002/ijc.2910500605 [PubMed]
  • 7. Loo JM, Scherl A, Nguyen A, Man FY, Weinberg E, Zeng Z, Saltz L, Paty PB, Tavazoie SF. Extracellular metabolic energetics can promote cancer progression. Cell. 2015; 160:393–406. https://doi.org/10.1016/j.cell.2014.12.018 [PubMed]
  • 8. Lange AM, Lo HW. Inhibiting TRK Proteins in Clinical Cancer Therapy. Cancers (Basel). 2018; 10:105. https://doi.org/10.3390/cancers10040105 [PubMed]
  • 9. Schmit K, Michiels C. TMEM Proteins in Cancer: A Review. Front Pharmacol. 2018; 9:1345. https://doi.org/10.3389/fphar.2018.01345 [PubMed]
  • 10. Vastrad B, Vastrad C, Tengli A, Iliger S. Identification of differentially expressed genes regulated by molecular signature in breast cancer-associated fibroblasts by bioinformatics analysis. Arch Gynecol Obstet. 2018; 297:161–83. https://doi.org/10.1007/s00404-017-4562-y [PubMed]
  • 11. Kim I, Choi S, Kim S. BRCA-Pathway: a structural integration and visualization system of TCGA breast cancer data on KEGG pathways. BMC Bioinformatics. 2018 (Suppl 1); 19:42. https://doi.org/10.1186/s12859-018-2016-6 [PubMed]
  • 12. Yang W, Ma J, Zhou W, Li Z, Zhou X, Cao B, Zhang Y, Liu J, Yang Z, Zhang H, Zhao Q, Hong L, Fan D. Identification of hub genes and outcome in colon cancer based on bioinformatics analysis. Cancer Manag Res. 2019; 11:323–38. https://doi.org/10.2147/CMAR.S173240 [PubMed]
  • 13. Tang D, Zhao X, Zhang L, Wang Z, Wang C. Identification of hub genes to regulate breast cancer metastasis to brain by bioinformatics analyses. J Cell Biochem. 2019; 120:9522–31. https://doi.org/10.1002/jcb.28228 [PubMed]
  • 14. Chen X, Cai S, Wang L, Zhang X, Li W, Cao X. Analysis of the function of MAGE-A in esophageal carcinoma by bioinformatics. Medicine (Baltimore). 2019; 98:e15774. https://doi.org/10.1097/MD.0000000000015774 [PubMed]
  • 15. Hu Y, Guo G, Li J, Chen J, Tan P. Screening key lncRNAs with diagnostic and prognostic value for head and neck squamous cell carcinoma based on machine learning and mRNA-lncRNA co-expression network analysis. Cancer Biomark. 2020; 27:195–206. https://doi.org/10.3233/CBM-190694 [PubMed]
  • 16. Han B, Yang X, Zhang P, Zhang Y, Tu Y, He Z, Li Y, Yuan J, Dong Y, Hosseini DK, Zhou T, Sun H. DNA methylation biomarkers for nasopharyngeal carcinoma. PLoS One. 2020; 15:e0230524. https://doi.org/10.1371/journal.pone.0230524 [PubMed]
  • 17. AsghariHanjani N, Vafa M. The role of IGF-1 in obesity, cardiovascular disease, and cancer. Med J Islam Repub Iran. 2019; 33:56. [PubMed]
  • 18. Cao Y, Nimptsch K, Shui IM, Platz EA, Wu K, Pollak MN, Kenfield SA, Stampfer MJ, Giovannucci EL. Prediagnostic plasma IGFBP-1, IGF-1 and risk of prostate cancer. Int J Cancer. 2015; 136:2418–26. https://doi.org/10.1002/ijc.29295 [PubMed]
  • 19. Al-Delaimy WK, Flatt SW, Natarajan L, Laughlin GA, Rock CL, Gold EB, Caan BJ, Parker BA, Pierce JP. IGF1 and risk of additional breast cancer in the WHEL study. Endocr Relat Cancer. 2011; 18:235–44. https://doi.org/10.1530/ERC-10-0121 [PubMed]
  • 20. Pal S, Yadav P, Sainis KB, Shankar BS. TNF-α and IGF-1 differentially modulate ionizing radiation responses of lung cancer cell lines. Cytokine. 2018; 101:89–98. https://doi.org/10.1016/j.cyto.2016.06.015 [PubMed]
  • 21. Xu L, Zhou R, Yuan L, Wang S, Li X, Ma H, Zhou M, Pan C, Zhang J, Huang N, Shi M, Bin J, Liao Y, Liao W. IGF1/IGF1R/STAT3 signaling-inducible IFITM2 promotes gastric cancer growth and metastasis. Cancer Lett. 2017; 393:76–85. https://doi.org/10.1016/j.canlet.2017.02.014 [PubMed]
  • 22. Dhifallah H, Aissi S, Njima M, Zakhama A, Kenani A. IGF1 polymorphisms and colon cancer risk in Tunisian population. Tunis Med. 2019; 97:1407–14. [PubMed]
  • 23. Feng X, Lin J, Xing S, Liu W, Zhang G. Higher IGFBP-1 to IGF-1 serum ratio predicts unfavourable survival in patients with nasopharyngeal carcinoma. BMC Cancer. 2017; 17:90. https://doi.org/10.1186/s12885-017-3068-0 [PubMed]
  • 24. Wang T, Dong XM, Zhang FL, Zhang JR. miR-206 enhances nasopharyngeal carcinoma radiosensitivity by targeting IGF1. Kaohsiung J Med Sci. 2017; 33:427–32. https://doi.org/10.1016/j.kjms.2017.05.015 [PubMed]
  • 25. Su T, Mi Y, Zhang L, Wang S, Lu H, Shi L, Sun H, Wu X, Zhang W, Zuo L, Zou J. Association between IL13 gene polymorphisms and susceptibility to cancer: a meta-analysis. Gene. 2013; 515:56–61. https://doi.org/10.1016/j.gene.2012.11.035 [PubMed]
  • 26. Tsai ST, Wang PJ, Liou NJ, Lin PS, Chen CH, Chang WC. ICAM1 Is a Potential Cancer Stem Cell Marker of Esophageal Squamous Cell Carcinoma. PLoS One. 2015; 10:e0142834. https://doi.org/10.1371/journal.pone.0142834 [PubMed]
  • 27. Reina M, Espel E. Role of LFA-1 and ICAM-1 in Cancer. Cancers (Basel). 2017; 9:153. https://doi.org/10.3390/cancers9110153 [PubMed]
  • 28. Tao L, Li Z, Lin L, Lei Y, Hongyuan Y, Hongwei J, Yang L, Chuize K. MMP1, 2, 3, 7, and 9 gene polymorphisms and urinary cancer risk: a meta-analysis. Genet Test Mol Biomarkers. 2015; 19:548–55. https://doi.org/10.1089/gtmb.2015.0123 [PubMed]
  • 29. Li Z, Ge H, Xie YG, Xie GY, Lv C. Matrix metalloproteinase-1 (MMP1) polymorphism is associated with lowered risk of nasopharyngeal carcinoma in Asian population. Cell Biochem Biophys. 2015; 71:999–1004. https://doi.org/10.1007/s12013-014-0299-4 [PubMed]
  • 30. Bai M, Chen X, Ba YI. CXCL10/CXCR3 overexpression as a biomarker of poor prognosis in patients with stage II colorectal cancer. Mol Clin Oncol. 2016; 4:23–30. https://doi.org/10.3892/mco.2015.665 [PubMed]
  • 31. Zhao M, Ma Q, Xu J, Fu S, Chen L, Wang B, Wu J, Yang L. Combining CXCL10 gene therapy and radiotherapy improved therapeutic efficacy in cervical cancer HeLa cell xenograft tumor models. Oncol Lett. 2015; 10:768–72. https://doi.org/10.3892/ol.2015.3281 [PubMed]
  • 32. Bronger H, Singer J, Windmüller C, Reuning U, Zech D, Delbridge C, Dorn J, Kiechle M, Schmalfeldt B, Schmitt M, Avril S. CXCL9 and CXCL10 predict survival and are regulated by cyclooxygenase inhibition in advanced serous ovarian cancer. Br J Cancer. 2016; 115:553–63. https://doi.org/10.1038/bjc.2016.172 [PubMed]
  • 33. Zhou H, Wu J, Wang T, Zhang X, Liu D. CXCL10/CXCR3 axis promotes the invasion of gastric cancer via PI3K/AKT pathway-dependent MMPs production. Biomed Pharmacother. 2016; 82:479–88. https://doi.org/10.1016/j.biopha.2016.04.069 [PubMed]
  • 34. Jafarzadeh A, Fooladseresht H, Nemati M, Assadollahi Z, Sheikhi A, Ghaderi A. Higher circulating levels of chemokine CXCL10 in patients with breast cancer: Evaluation of the influences of tumor stage and chemokine gene polymorphism. Cancer Biomark. 2016; 16:545–54. https://doi.org/10.3233/CBM-160596 [PubMed]