Key radioresistance regulation models and marker genes identified by integrated transcriptome analysis in nasopharyngeal carcinoma.

Nasopharyngeal carcinoma (NPC) is a malignancy that is endemic to China and Southeast Asia. Radiotherapy is the usual treatment, however, radioresistance remains a major reason for failure. This study aimed to find key radioresistance regulation models and marker genes of NPC and clarify the mechanism of NPC radioresistance by RNA sequencing and bioinformatics analysis of the differences in gene expression profiles between radioresistant and radiosensitive NPC tissues. A total of 21 NPC biopsy specimens with different radiosensitivity were analyzed by RNA sequencing. Differentially expressed genes in RNA sequencing data were identified using R software. The differentially expressed gene data derived from RNA sequencing as well as prior knowledge in the form of pathway databases were integrated to find sub-networks of related genes. The data of RNA sequencing with the GSE48501 data from the GEO database were combined to further search for more reliable genes associated with radioresistance of NPC. Survival analyses using the Kaplan-Meier method based on the expression of the genes were conducted to facilitate the understanding of the clinical significance of the differentially expressed genes. RT-qPCR was performed to validate the expression levels of the differentially expressed genes. We identified 1182 differentially expressed genes between radioresistant and radiosensitive NPC tissue samples. Compared to the radiosensitive group, 22 genes were significantly upregulated and 1160 genes were downregulated in the radioresistant group. In addition, 10 major NPC radiation resistance network models were identified through integration analysis with known NPC radiation resistance-associated genes and mechanisms. Furthermore, we identified three core genes, DOCK4, MCM9, and POPDC3 among 12 common downregulated genes in the two datasets, which were validated by RT-qPCR. The findings of this study provide new clues for clarifying the mechanism of NPC radioresistance, and further experimental studies of these core genes are warranted.

INTRODUCTION

Nasopharyngeal carcinoma (NPC) is an epithelial head and neck tumor with marked geographical variation in distribution. It is mostly prevalent in China and Southeast Asia. According to the 2018 global cancer statistics, there were about 129 079 new cases of NPC in 2018, accounting for 0.7% of all new cancer patients. Intensity‐modulated radiotherapy (IMRT) is the standard treatment method for patients with non‐metastatic NPC, which enables good tumor control. , However, 10%–20% of NPC patients develop recurrence after radiotherapy due to radioresistance. Notably, the prognosis of patients with recurrent NPC is still very poor. , Therefore, it is imperative to clarify the molecular mechanism underlying the radioresistance of NPC.

Previous studies have found some molecules and biological processes relevant to the radioresistance of NPC by analyzing radioresistant NPC cells and radiosensitive NPC cells. , , , For example, Chang et al. compared two radioresistant NPC cell lines with their corresponding parental cell lines by cDNA microarray and found that seven genes, including gp96/hsp90b1 and GDF15, were associated with the radioresistance of NPC. In addition, Li et al. identified 15 differentially expressed miRNAs and 372 differential mRNAs associated with NPC radioresistance by comparing radioresistant NPC CNE2‐IR cells with radiosensitive NPC CNE2 cells. However, to date only limited studies have focused on the analysis of the difference between radioresistant and radiosensitive clinical NPC tissues.

Therefore, in this study, we exploited RNA sequencing technology to explore the differences in gene expression profiles between radioresistant and radiosensitive NPC biopsy specimens. Together with known genes and mechanisms involved in the radioresistance of NPC, we set out to define the major radioresistance network models through dataset integration and to make inferences on the key genes. Furthermore, we combined the RNA sequencing results with the existing data in the Gene Expression Omnibus (GEO) database to filter out reliable molecular markers and investigate the impact of these markers on the survival rates of patients. The major regulation models described in this study provide new clues for elucidating the mechanism of radioresistance of NPC and may be used to develop new predictors of NPC radioresponse and patient prognosis.

MATERIALS AND METHODS

Patients and tissue samples

A total of 21 fresh frozen NPC biopsy specimens from the Sun Yat‐sen University Cancer Center were retrospectively collected from January 2013 to February 2017. All patients were pathologically diagnosed with non‐metastatic NPC, and tumor tissue samples were obtained before radiotherapy. The samples were immediately frozen in liquid nitrogen and stored. All patients received radiotherapy and platinum‐based chemotherapy. According to the response to radiotherapy, we divided the 21 NPC patients into two groups: a radiosensitive group with 14 patients and a radioresistant group with 7 patients. Here, radiosensitive patients were defined as NPC patients with no residual lesions after 6 weeks of radiotherapy and no recurrence within 5 years after radiotherapy, and radioresistant patients were defined as NPC patients with residual lesions after radiotherapy for more than 6 weeks or NPC patients with recurrence within 1 year after radiotherapy. All patients were re‐staged according to the 8th edition of the American Joint Committee on Cancer (AJCC) Staging Manual. Detailed clinical features of the 21 patients are shown in Table 1. This study was approved by the Ethical Review Committee of Sun Yat‐sen University Cancer Center, and informed consent was obtained from all patients.

TABLE 1

The clinical characteristics of nasopharyngeal carcinoma samples

	Radiosensitive (n=14)	Radioresistant (n=7)	P value
Age (mean ±SD)	45.57±7.76	40.6±8.32	0.190
Sex			0.006
Male	14	3
Female	0	4
T stage			0.374
T1	0	1
T2	2	1
T3	10	3
T4	2	2
N stage			0.677
N0	1	1
N1	5	2
N2	6	4
N3	2	0
TNM stage			0.866
II	1	1
III	9	4
IV	4	2

Total RNA extraction

Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. The purity and concentration of RNA were detected using a NanoDrop ND‐1000 spectrophotometer (NanoDrop, Wilmington, DE, USA). The integrity of RNA was verified using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA) with an RNA integrity number (RIN) >7.0.

Library construction and RNA sequencing

Library construction and RNA sequencing were performed by LC Bio Inc. (Hangzhou, China). Approximately 5 μg of total RNA was extracted, and then ribosomal RNA (rRNA) was removed from the total RNA using the Ribo‐Zero™ rRNA Removal Kit (Illumina, San Diego, USA). The remaining RNA was fragmented and then synthesized into first‐strand cDNA using reverse transcriptase and random primers. Second‐strand cDNA synthesis was then performed using E. coli DNA polymerase I, RNase H, and dUTP. Next, the cDNA strands were end‐repaired and added an “A” base. They were then ligated to the indexed adapters which contained a “T” base overhang. After dUTP strand degradation by the treatment of the UDG enzyme, the cDNA products were amplified by polymerase chain reaction (PCR) for the formation of a library with a fragment size of 300 bp (±50 bp). Finally, paired‐end sequencing was performed using Illumina X Ten (LC Bio, China).

Principal component analysis

To clearly evaluate the similarities and differences between samples and determine whether samples were grouped correctly, principal component analysis (PCA) of two cohorts was conducted using the ellipse package in R software (version 3.5.0). PCA is a dimensionality reduction method that is used to reduce the dimensionality of large datasets, by transforming a large set of variables into a smaller one that still contains most of the information in the large set.

Identification of differentially expressed genes

Clean sequencing reads were aligned with the index built from the human (hg37) genome, and the high‐quality reads were mapped to the reference by HISTA2 v2.1.0, and then the FPKM (fragments per kilobase per million) value of the genes and isoforms were calculated using StringTie v2.1.4 using a combination of Illumina and full‐length transcript‐based annotations.

Differentially expressed genes in RNA sequencing data were identified using the genefilter package in R software (version 3.5.0). Genes with a fold change greater than 1.5 (FC>1.5) and P < 0.05 between the radioresistant and radiosensitive groups were considered significant.

GO and KEGG pathway enrichment analyses

To understand the functional properties of the differentially expressed genes, Gene Ontology (GO, http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) pathway analyses were performed using the OmicsBean workbench (http://www.omicsbean.cn), an online multiple omics data analysis application. GO analysis was utilized to characterize genes and gene products in terms of cellular component (CC), biological process (BP), and molecular function (MF). KEGG pathway analysis was performed to further identify the pathways in which the differentially expressed genes underwent significant enrichment, thus predicting the potential functions of the differentially expressed genes. In addition, the pathway activation strength (PAS) prediction algorithm implemented in the omicsbean workbench was used to predict the effects of the processes identified by GO enrichment analysis. A positive value of PAS indicates the activation of a signaling pathway, while a negative value indicates the inhibition of a signaling pathway.

Summary of the known genes and mechanisms associated with radioresistance

To summarize the influencing factors and regulatory mechanisms of tumor radiosensitivity and radioresistance that have been known, we retrieved relevant literature from January 1990 to October 2020 using PubMed with the keywords of “radiosensitivity,” “radiation sensitivity,” “radioresistance,” “radiation resistance,” and “cancer.” The representative genes and mechanisms retrieved are listed in Table 2.

TABLE 2

Known genes and mechanisms associated with radioresistance

Radioresistance‐related mechanisms	Radioresistance‐associated genes	Refs
Enhanced DNA damage repair	MRE11, RAD50, NBS1, ATM, ATR, RAD51, BRCA1, BRCA2, DNA‐PK, XRCC4, LIG4, H2AX, MDMX, MDM2, MDC1, 53BP1, TLK1, Rad9, ATF2, and SMC1	48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59
Altered cell cycle	Chk1, Chk2, CDC25A, CDK2, CDC25C, CDK1, p21, p16, GADD45, NF‐kappa‐B, and FANCD2	60, 61, 62, 63, 64, 65, 66, 67, 68
Evasion of apoptosis	TP53, Bcl2, Bax, FAS, TNF, TRAIL, Livin, XIAP, CIAP1, CIAP2, Survivin, Smac, Caspase, RelB, CREB, and SAPK	69, 70, 71, 72, 73, 74
Hypoxia	HIF1	75, 76, 77
Angiogenesis	VEGF	78

Construction of “hub” sub‐network models for the radioresistant system

Through analyzing the relationships between the sequenced differential genes and the reported marker genes related to radiosensitivity, that is, including protein–protein interactions and pathways that genes collectively participate in the “hub” sub‐network models related to radiosensitivity were constructed. Specifically, the sub‐network models were generated by the Cytoscape web application, based on the information obtained from four levels of functional analysis: fold change of genes/proteins, protein–protein interactions, KEGG pathway enrichment, and biological process enrichment. The STRING database (search tool for the retrieval of interacting genes/proteins, http://string‐db.org) was used to analyze protein–protein interactions. Go and KEGG analyses were performed for pathway enrichment analysis.

Combination analyses of differentially expressed genes with GEO data

The microarray dataset GSE48501 based on the Affymetrix Human Genome U133 Plus 2.0 Array platform was downloaded from the GEO database (http://www. ncbi.nlm.nih.gov/geo). The dataset contained the mRNA expression profiles of two samples of radioresistant NPC CNE2‐IR cells and two samples of radiosensitive NPC CNE2 cells. NPC CNE2‐IR cells were derived from the poorly differentiated NPC cell line CNE2 by treating the cells with four rounds of a sublethal dose of radiation. Our RNA sequencing results were combined with the data from GSE48501 to identify overlapping differentially expressed genes. The potential functions of the overlapping differentially expressed genes were further analyzed.

RT‐qPCR validation for the expression of the differentially expressed genes

The expression of the differentially expressed genes was detected by reverse transcription‐quantitative polymerase chain reaction (RT‐qPCR). Specific primer sequences for the genes are shown in Table 3. Total RNA (2 µg) was reverse transcribed into cDNA using a Reverse Transcription Kit according to the manufacturer's protocol. Then, using cDNA as the template, qPCR was performed using qPCR Mix under the following reaction conditions: the initial denaturation step was 95℃ for 10 min, followed by 40 cycles of 95℃ for 10 seconds and 60℃ for 1 min. The internal reference gene was GAPDH, and the relative expression levels of the genes were calculated by the 2‐ΔΔCt method.

TABLE 3

A list of primers used in this study

Gene	Primer Sequence (5’−3’)
DOCK4	F:ATTCCAGAGAGCCAGGAGGT
	R:TGACGTTCTCTCCACCCAGA
MCM9	F:AGGTTCTGGAGTTTGAGCGG
	R:ACAAGCCTGAGAGGCAAGTG
POPDC3	F:TGCACAACCTGGAAGCAAGA
	R:AGAAAACCCAACCCCAGCAA
GAPDH	F:GCATCCTGGGCTACACTGAG
	R:AAAGTGGTCGTTGAGGGCAA

Statistical analysis

Statistical analysis was conducted using R 3.0 (http://www.r‐project.org/). The results are presented as mean ±SD. For differential expression analysis, Student's t‐test between groups was used. The rates of overall survival (OS) and progression‐free survival (PFS) were calculated using the Kaplan–Meier method, and the differences in survival rates between patients with different gene expression levels were compared using the log‐rank test. P values <0.05 were considered statistically significant.

RESULTS

Differentially expressed genes between radiosensitive and radioresistant NPC tissues

In our dataset, the two independent cohorts (radioresistant group and radiosensitive group) can be clearly separated into two clusters with principal component analysis (PCA) (Figure 1A), suggesting that the radio impacts on the transcriptome exhibit expressional closeness within each group.

FIGURE 1

Identification and hierarchical clustering of differentially expressed genes. RR for radioresistant and RS for radiosensitive. (A) Principal component analysis (PCA) of two cohorts. (B) Volcano plot of differentially expressed genes between radioresistant and radiosensitive groups. The cutoff criteria were fold change >1.5 and P < 0.05. The red dots represent the upregulated genes and the blue dots signify the downregulated genes. The black dots indicate the genes with a fold change <1.5 and/or P > 0.05

In total, we identified 1182 differentially expressed genes with filter criteria: fold change >1.5 and P < 0.05. Compared to the radiosensitive group, 22 genes were significantly upregulated and 1160 genes were downregulated in the radioresistant group (Figure 1B).

Gene ontology and KEGG analyses of the differentially expressed genes

Gene ontology analysis of the 1182 differentially expressed genes enriched in a total of 7153 BP, 951 CC, and 1490 MF terms (Figure 2A), among which 4201, 537, and 567 terms were significantly enriched (P < 0.05), respectively. KEGG pathway enrichment analysis revealed a total of 286 pathways (Figure 2A), among which 33 pathways met the P < 0.05 criteria.

FIGURE 2

Enrichment analysis of the differentially expressed genes between the radioresistant and radiosensitive groups. (A) Gene ontology enrichment and KEGG pathway enrichment analysis of the differentially expressed genes. Blue and orange bars indicate enriched total terms and terms exhibiting statistical significance (P < 0.05) in biological process, cell component, molecular function, and KEGG pathway, respectively. (B) Significant (P < 0.05) biological processes enriched by Gene ontology analysis. (C) Predicted top activated and inhibited functional processes based on pathway activation strength (PAS) scores. Brown bars and green bars represent the degree of pathway activation or inhibition, respectively

Among the significantly enriched biological processes, regulation of the metabolic process, regulation of cell communication, apoptotic process, regulation of programmed cell death, and regulation of cell cycle were enriched, playing a leading role in radioresistant events. In addition, these differentially expressed genes were also significantly involved in cell migration, Notch signaling pathway, response to cytokine, mitotic DNA damage checkpoint, and regulation of DNA repair (Figure 2B).

Using the PAS prediction algorithm for each process, it was revealed that in the radioresistant condition, regulation of the metabolic process, apoptotic process, regulation of programmed cell death, and cell cycle regulation were strongly inhibited (Figure 2C).

Based on the KEGG pathway analysis, three metabolic pathways, including fatty acid elongation, and glycosaminoglycan biosynthesis, were significantly enriched. Moreover, three environment‐related pathways including the Wnt signaling pathway, as well as eight cellular processes including endocytosis, apoptosis, lysosome, and focal adhesion, were enriched. Eight organismal system‐related pathways including T‐cell receptor signaling pathway and neurotrophin signaling pathway were also enriched (Figure 3).

FIGURE 3

Significantly (P < 0.05) enriched KEGG pathways and classification.

“Hub” sub‐network models for the radioresistant system

To find sub‐networks of related genes implicated by multiple forms of biological evidence, we integrated the differentially expressed gene data derived from RNA sequencing, as well as prior knowledge in the form of pathway databases. After integrating the sequenced differential gene data with the reported gene data related to radiosensitivity that we retrieved (Table 2), it was found that 10 cancer features, including DNA damage pathway (Figure 4A), cell cycle pathway (Figure 4B), DNA repair pathway (Figure 4C), apoptosis pathway (Figure 4D), stemness pathway (Figure S1A), chromatin pathway (Figure S1B), radiation response pathway (Figure S1C), metal iron response pathway (Figure S1D), epithelial–mesenchymal transition (EMT) pathway (Figure S1E), and cytokine and degranulation pathway (Figure S1F), were enriched. Inspecting each sub‐network through our data revealed several newly discovered common and unique differentially expressed genes related to radioresistance with high connectivity with known genes. It included HIPK2, MCM9 (DNA damage), MAP4K4 (cell cycle), MRNIP (DNA repair), IL2RA (apoptosis), and THBS1 (EMT), implicating several new targets for investigation.

FIGURE 4

“Hub” sub‐network models related to radioresistance. (A‐D) The main NPC radioresistance models including DNA damage pathway (A), cell cycle pathway (B), DNA repair pathway (C), and apoptosis pathway (D) were constructed by integrating the differentially expressed gene data (left half node of the cycle nodes) with the reported genes (right half node of the cycle nodes). Circle nodes indicate genes, with the right half of the circle colored red representing the gene as a marker gene, the left half colored red representing the gene upregulated in differential expression, and the left half colored green representing the gene downregulated in differential expression. Rectangles indicate KEGG pathways or biological processes. Pathways were colored with gradient color from yellow to blue, with smaller p values in yellow and larger p values in blue

Combination analyses of the differentially expressed genes with GEO data

To further search for more reliable genes associated with radioresistance of NPC, we combined the data of RNA sequencing with the GSE48501 data from the GEO database. The data analysis process is shown in Figure 5A. A total of 12 overlapping differentially expressed genes were identified between the two datasets (Figure 5B‐C ), including MAP4K4, DOCK4, NFE2L3, THBS1, EOMES, MCM9, SERPINI1, ARTN, MRPS9, FSIP1, RRP15, and POPDC3. These 12 genes were all downregulated in the radioresistant group (Figure 5d).

FIGURE 5

Identification of the overlapping differentially expressed genes between the data of RNA sequencing and the data of GSE48501. RR for radioresistant and RS for radiosensitive. (A) A flowchart of identifying the overlapping differentially expressed genes. (B) A Venn diagram of the overlapping downregulated expressed genes in both the data of RNA sequencing and the data of GSE48501. (C) A Venn diagram of the overlapping upregulated expressed genes in both the data of RNA sequencing and the data of GSE48501. (D) Heatmap of the 12 overlapping differentially expressed genes between the data of RNA sequencing and the data of GSE48501. The horizontal band at the top: cyan: RR, radioresistant group; pink: RS, radiosensitive group. Each row represents a single gene. Green indicates low expression; red indicates high expression

Survival analysis of the differentially expressed genes

To facilitate the understanding of the clinical significance of the differentially expressed genes, survival analyses using the Kaplan–Meier method based on the expression of the genes were conducted. According to the expression levels of the genes in the NPC specimens, the patients were divided into two levels: high expression and low expression. The results of the Kaplan–Meier analyses revealed that patients with a lower level of DOCK4 exhibited significantly shorter PFS (Figure 6A, P < 0.05). Similar results were shown in gene MCM9 (Figure 6B), POPDC3 (Figure 6C), ARTN (Figure S2A), MRPS9 (Figure S2B), and SERPIN1 (Figure S2C). In addition, the results of the Kaplan–Meier analyses revealed that the lower expression of ARTN, FSIP1, MCM9, MRPS9 was significantly associated with poorer OS (Figure S2D‐G, P < 0.05).

FIGURE 6

Analysis of the three differentially expressed genes. RR for radioresistant and RS for radiosensitive. (A‐C) Progression‐free survival (PFS) for DOCK4, MCM9, and POPDC3, respectively. The expression value of the gene was divided into two parts: low expression (0%‐50%) and high expression (50%‐100%). (D‐E) Validation of the three differentially expressed genes in 35 NPC biopsy specimens with different radiosensitivity by RT‐qPCR (20 radiosensitive samples and 15 radioresistant samples). Mann–Whitney test was performed to calculate significance. *P < 0.05

Validation of the differentially expressed genes by RT‐qPCR

RT‐qPCR was performed to validate the expression levels of the differentially expressed genes between 20 radiosensitive NPC specimens and 15 radioresistant NPC specimens. The results of RT‐qPCR revealed that the expression levels of DOCK4 (P < 0.05), MCM9 (P < 0.05), and POPDC3 (P < 0.05) were all significantly downregulated in the radioresistant NPC specimens when compared to those of the radiosensitive NPC specimens (Figure 6D‐F).

DISCUSSION

At present, radiotherapy resistance has become a major obstacle to the success of NPC treatment. Therefore, increasing attention has been paid to the mechanism of radioresistance in NPC. In recent years, a large number of studies have identified important molecules and biological processes associated with radioresistance by comparing the differences in the expression profiles of genes and other molecules between radiosensitive and radioresistant NPC cell lines. , , , However, the samples for these studies were all from NPC cell lines cultured in vitro. Few studies have focused on the differences between radioresistant NPC and radiosensitive NPC biopsy tissues so far.

In this study, therefore, we exploited the technology of RNA sequencing to compare the differences in gene expression levels between radioresistant NPC tissues and radiosensitive NPC tissues. As is well known, RNA sequencing, which is considered to be a revolutionary tool for transcriptomics, has high levels of accuracy and reproducibility for detecting gene expression levels. The results of this study showed that a total of 22 genes were significantly upregulated and 1160 genes were downregulated in the radioresistant group when compared with the genes in the radiosensitive group. With the GO enrichment analysis, we found that the most enriched pathways were concentrated on regulation of the metabolic process, the apoptotic process, regulation of cell cycle, which were consistent with the radioresistance‐related pathways reported in many previous studies. , , These differentially expressed genes were also found to be involved in Notch signaling pathway, response to cytokine, mitotic DNA damage checkpoint, and regulation of DNA repair. Subsequently, the results of the PAS prediction algorithm for each process revealed that under radioresistant conditions, processes including regulation of the metabolic process, apoptotic regulation, and cell cycle regulation were strongly inhibited which might partially explain the underlying mechanism of radioresistance. According to the KEGG pathway analysis, three environment‐related pathways including the Wnt signaling pathway, as well as eight cellular processes including endocytosis, apoptosis, lysosome, and focal adhesion, were enriched. As reported, the Wnt signaling pathway is one of the important pathways that regulate the proliferation, differentiation, and migration of cells, and dysregulation of the Wnt signaling pathway is closely associated with the development of a variety of tumors such as lung cancer, liver cancer, and breast cancer. Moreover, three immune‐related pathways, including T‐cell receptor signaling pathway were also enriched, suggesting that the immune status of T cells might be associated with radiation resistance.

By integrating the differentially expressed gene data acquired by RNA sequencing with that of the known genes associated with radiosensitivity reported in the previous literature, we established “hub” sub‐networks of genes related to radioresistance. A total of 10 cancer features, including DNA damage pathway, cell cycle pathway, DNA repair pathway, apoptotic pathway, EMT pathway, chromatin organization pathway, cytokine production and degranulation pathway, stem cell differentiation pathway, and metal iron pathway, were finally enriched. After the subsequent detailed analysis of each pathway, several common participants with high connectivity, such as HIPK2, MCM9, and THBS1, were identified. HIPK2 is considered to be a crucial regulator for targeting apoptosis because it phosphorylates the tumor suppressor p53 in response to DNA damage. , It has been shown that HIPK2 knockdown could induce chemoresistance as well as tumor growth in vivo. HIPK2 has been reported to be involved in the hypoxic response as a co‐suppressor of hypoxia‐inducible factor‐1α (HIF‐1α), which is a major factor that regulates the transcription of angiogenesis and invasion‐related genes. THBS1, also known as TSP1, is a member of the thrombospondins (TSPs) family, and its encoded product is a matricellular protein that has the property of limiting angiogenesis by direct effects on endothelial cell migration, proliferation, survival, and apoptosis through CD36, CD47, and integrins. , Given its role in delaying angiogenesis, THBS1 has been shown to suppress tumor growth and has also been found to be positively associated with patient survival in several cancers, such as lung, bladder, gastric, and colon cancers.

Furthermore, we combined our RNA sequencing data of radioresistant and radiosensitive NPC tissues with that of radioresistant and radiosensitive NPC cells in the GEO database to find more reliable core genes. A total of 12 overlapping genes were identified finally and the gene expression in three of them including DOCK4, MCM9, and POPDC3 were validated successfully with RT‐qPCR. DOCK4, a member of the dedicator of cytokinesis (DOCK) family, functions as a guanine nucleotide exchange factor (GEF), converting inactive GDP‐bound small GTPases into their active GTP‐bound form and is involved in the regulation of adherens junctions between cells. DOCK4 has been reported to interact with RAC1, which is associated with chemoresistance, radioresistance, resistance to targeted therapies, and immune evasion, implicating a potential promoting mechanism of radioresistance by downregulated DOCK4. A recent study reported that the overexpression of DOCK4 suppresses the tumorigenicity of glioblastomas (GBM) stem‐like cells, and an increased level of DOCK4 predicts improved patient survival of GBM. Another study showed that DOCK4 expression level is downregulated in paclitaxel‐resistant breast cancers and lncRNA AC073284.4 might sponge miR‐18b‐5p to attenuate the invasion, metastasis, and epithelial–mesenchymal transition of breast cancer cells by upregulating DOCK4 expression. Similarly, DOCK4 was downregulated in radioresistant NPC in our study. More importantly, the results of the Kaplan–Meier analyses revealed that patients with a lower level of DOCK4 exhibited significantly shorter PFS in our study. MCM9, a member of the mini‐chromosome maintenance (MCM) family, has been shown to play a critical role in DNA replication and repair. Several studies have found that the complex consisting of MCM9 and its homolog MCM8 can promote homologous recombination‐mediated DNA repair by facilitating RAD51 recruitment to sites of DNA damage and interacting with the MRE11–RAD50–NBS1 complex. , In addition, MCM9 has been reported to possess a helicase activity which is required for efficient DNA mismatch repair(MMR), and cells with knockdown of MCM9 exhibit microsatellite instability and MMR deficiency. Our analysis showed that MCM9 was downregulated in radioresistant NPC, suggesting that the low level of MCM9 might be associated with the radioresistance of NPC. POPDC3, also known as POP3, encodes a transmembrane protein that can facilitate cyclic adenosine monophosphate (cAMP)‐mediated signaling. Unlike its isoform POPDC1 which acts as a tumor suppressor, POPDC3 has been found to play distinct roles in different cancer types. For example, researchers found that knockdown of POPDC3 significantly increased the migration and invasion of gastric cancer cells. Besides, low expression of POPDC3 was also reported to be associated with metastasis and poor prognosis of gastric cancers. In contrast, a recent study found that high POPDC3 expression was significantly associated with poor prognosis in head and neck squamous cell carcinoma. To date, the function of these three genes in NPC has not been reported yet. Therefore, it is necessary to further study the role of these three genes in the radioresistance of NPC in the future.

In conclusion, we analyzed the differentially expressed genes between radioresistant and radiosensitive NPC tissue samples by RNA sequencing and bioinformatics analysis in this study. In addition, 10 major NPC radiation resistance network models were identified through integration analysis with known NPC radiation resistance‐associated genes and mechanisms. Furthermore, we identified three core genes, DOCK4, MCM9, and POPDC3, that may be involved in the radioresistance of NPC. The findings of this study provide new clues for clarifying the mechanism of NPC radioresistance, and further experimental studies of these core genes are warranted.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTION STATEMENT

F.H performed study concept and design; Z.S, XH.W, and JY.W provided analysis and interpretation of data, and statistical analysis; Z.S, XH.W, JY.W, and J.W performed experimental verification and writing of the paper; RD.H, CY.C, ML.D, and HY.W performed review and revision of the paper; X.L provided technical support. All authors read and approved the final paper.

ETHICS STATEMENT

The current study was approved by the Clinical Research Ethics Committee of Sun Yat‐sen University Cancer Center and was performed in accordance with the Declaration of Helsinki.

FUNDING STATEMENT

The authors received no specific funding for this work.

Figure S1‐S2

Click here for additional data file.