Irradiation-induced dynamic changes of gene signatures reveal gain of metastatic ability in nasopharyngeal carcinoma

Abstract

Nasopharyngeal carcinoma (NPC), arising from the nasopharynx epithelium, is prevalent among South and East Asia. The radiotherapy is the primary treatment for NPC patients. However, the acquired radioresistance dramatically diminishes the therapeutic effect of radiotherapy. Meanwhile, recurrence and metastasis always occur in line with the radioresistance, but the underlying mechanisms are still unclear. In this study, we established two radioresistant NPC cell lines, CNE1R and SUNE1R, by sequentially irradiated parental CNE1 and SUNE1 cells up to a clinical treatment dose of 72 Gy. A transcriptome profile analysis of CNE1R and CNE1 reveals that activated oncogenic pathways are highly enriched in CNE1R. As the result, CNE1R showed higher proliferation rate but lower apoptosis rate after irradiation, and enhanced metastasis ability in comparison with CNE1. Significantly, a group of metastasis associated genes were increased in CNE1R while the irradiation proceeded, including several matrix metallopeptidase (MMP) members, especially MMP10 and MMP13. With further analysis, we found both MMP10 and MMP13 are highly upregulated in metastatic head and neck cancer specimens compared to non-metastatic ones. More importantly, patients with lower expression of both MMP10 and MMP13 showed a better five-year survival than the double high group. Our findings unveiled the potential mechanisms of radioresistance related metastasis in NPC patients, and the increase of MMP10 and MMP13 may serve as high risk factors for metastasis during radiotherapy.

Keywords: Nasopharyngeal carcinoma; head and neck squamous cell carcinoma; matrix metallopeptidase; metastasis; radioresistance.

Figure 1

Establish a radioresistant NPC cell line CNE1R. (A) Schematic diagram of the establishment of a radioresistant CNE1R cell line. (B) Growth of CNE1R is sustained after IR. CNE1 and CNE1R plated in 24-well culture plates were exposed to IR at 2 Gy and cell growth was monitored daily by counting cell numbers. (C) CNE1R is resistant to IR in a colony-forming assay. CNE1 and CNE1R plated in 6-well culture plates were exposed to IR at 8 or 10 Gy, or without IR. (D) The survival fraction of the colony-forming assay in (C) was calculated and normalized to no IR group. All data are representative of three independent experiments and values are shown in mean ± SEM (B, D), NS means no significant difference, *P < 0.05, **P < 0.01 (Student’s t test).

Figure 2

CNE1R expression profile shows oncogenic pathway activations. A. Volcano plot is used to show different expression genes between CNE1 and CNE1R. Red dots indicate upregulated genes in CNE1R, and green dots indicate downregulated genes in CNE1R. Fold change ≥ 2 or ≤ -2 and q-value < 0.001 were significant. B. 328 upregulated genes in CNE1R was analyzed by KEGG pathway. Gene numbers involved in each pathway was showed as black dots, combined score was showed as histogram in red. C. 513 downregulated genes in CNE1R was analyzed by KEGG pathway. Gene numbers involved in each pathway was showed as black dots, combined score was showed as histogram in blue.

Figure 3

Both CNE1R and CNE1 expression profiles changed in a dynamic manner after irradiation. A-F. Gene expression time course analysis was performed by Mfuzz, and clustered into six groups. Top ranking KEGG pathways in CNE1 and CNE1R were ranked by p-value (-log10). Red dots indicate overlapped pathways in CNE1 and CNE1R. Black dots indicate unique pathways in each cell line.

Figure 4

The irradiation activates signature oncogenic signaling in resistant CNE1R cells. (A) Circos visualization of overlapping different expression genes (DEGs) between CNE1 and CNE1R during different time points after irradiation, extended by informative enriched terms. Gene pairs (from different lists) overlap are linked by purple curves, while genes fall into the same enriched term are linked by blue curves. (B) Heatmap visualization of pathway enrichments. A subset of representative enriched terms based on the GO clusters and ranked by -log₁₀ (p-value). Five different clusters were annotated base on GO items. (C) Hallmark of KRAS activation was enriched in CNE1R compared to CNE1 by gene set enrichment analysis (GSEA). (D) GSEA analysis showed CNE1 transcriptome enriched in negative regulation of response to wounding when compared to CNE1R. (E) Hallmark of apoptosis was enriched in CNE1 compared to CNE1R 48 post-irradiation in gene set enrichment analysis (GSEA). (F) CNE1R is resistant to IR-induced cell death. CNE1 and CNE1R cells were treated with or without 8 Gy IR, and apoptosis was detected by flow cytometer with Annexin V and 7-AAD staining. Cells with PE Annexin V positive/7-AAD negative or PE Annexin V positive/7-AAD positive were defined as apoptosis cells, shown as mean ± SEM. (G) Inflammatory response was enriched in CNE1 compared to CNE1R 48 post-irradiation in gene set enrichment analysis (GSEA). (H) qRT-PCR was used to determine the expression of representative inflammatory response associated genes enriched in (G). All gene expressions were normalized to CNE1. All data are representative of three independent experiments and values are shown in mean ± SEM (F and H). **P < 0.01 (F, Student’s t test; H, ANOVA).

Figure 5

Radioresistant CNE1R cells show enrichments in the metastasis associated gene sets. A. Venn diagram showed overlapped upregulated genes in CNE1 and CNE1R without IR or 6 h post-IR, KEGG pathway and GO analysis were used for analyzing all overlapped genes. B. Venn diagram showed overlapped upregulated genes in CNE1 and CNE1R without IR or 48 h post-irradiation, KEGG pathway and GO analysis were used for analyzing all overlapped genes. C. Venn diagram showed overlapped upregulated genes in CNE1 and CNE1R without IR or 6 h and 48 h post-irradiation, KEGG pathway and Wiki pathway analysis were used for analyzing all overlapped genes. D. Heatmap showed pro-resistance genes in CNE1R: upregulated in CNE1R; no significant change at 6 h post-irradiation; upregulated at 48 h after IR. E. Heatmap showed anti-resistance genes in CNE1R: downregulated in CNE1R; no significant change at 6 h post-irradiation; downregulated at 48 h after IR. F. String 10 was used to analysis interaction proteins of radioresistance associated genes. Two main groups of genes involved in metastasis associated pathways: Regulation of migration (in Red), and extracellular matrix organization (in Blue).

Figure 6

Radioresistant CNE1R cells gain metastasis ability. (A) CNE1R and SUNE1R cells show higher wound healing associated genes expression. qRT-PCR was used for detecting gene expression, all gene expressions were normalized to which of CNE1 or SUNE1. (B) Metastasis associated genes MMP10, MMP13, LMO1 and ELMO1 increased gradually during the process of sequential irradiation in CNE1R cells. (C) CNE1R cells have increased migration ability in the wound healing assay. Wounding area was calculated and normalized to day 0. (D) SUNE1R cells have increased migration ability in the wound healing assay. Wounding area was calculated and normalized to day 0. (E) The loss of MMP10 and MMP13 attenuates CNE1R wound healing ability. Wounding area was calculated and normalized to day 0. (F) The loss of MMP10 and MMP13 attenuates SUNE1R wound healing ability. Wounding area was calculated and normalized to day 0. (G) CNE1R cells gain metastatic ability in the nude mice lung metastasis model, IR enhances the metastatic tumor incidence. H&E staining of mouse lungs are shown. Metastatic node is shown by rectangle frame and magnification is shown in the right corner. Scale bar, 500 um. The lung metastatic tumor incidence of different group is shown in the table. All data are representative of three independent experiments and values are shown in mean ± SEM (A-F). **P < 0.01 (A, B, E and F, ANOVA; C and D, Student’s t test; G, Chi-square test).

Figure 7

Gene signatures indicated gain of metastasis ability after the irradiation. (A) MMP10 and MMP13 showed significant predication ability of overall survivals in HNSCC patients by ROC analysis. (B) MMP10 expression was positively correlated with MMP13 in HNSCC patients (TCGA, n = 499). (C) MMP10 were upregulated in metastatic (n = 11) compared to non-metastatic (n = 11) HNSCC tumors (GSE9349). (D) MMP13 were upregulated in metastatic (n = 11) compared to non-metastatic (n = 11) HNSCC tumors (GSE9349). (E) MMP10 expression was positively correlated with MMP13 in HNSCC tumors (GSE9349). (F) Kaplan Meier-plotter showed MMP10^LowMMP13^Low group patients had a better five-year survival than MMP10^highMMP13^high group. *P < 0.05, **P < 0.01 (C and D, Student’s t test). Correlation analysis was calculated by Pearson correlation coefficient (B and E).

Figure 8

The model of radioresistance related functional alterations after irradiation. After irradiation, radioresistant cells have altered gene expression profiles and acquire higher metastasis ability.