Population and single‑cell transcriptome analyses reveal diverse transcriptional changes associated with radioresistance in esophageal squamous cell carcinoma.

Esophageal squamous cell carcinoma (ESCC) is a tumor composed of heterogeneous cells that easily become radioresistant, which leads to tumor recurrence. The most commonly used treatment for ESCC is fractionated irradiation (FIR) therapy that utilizes ionizing radiation to directly induce cytotoxic cell death. However, this treatment may not be able to eliminate all cancer cells due to high adaptive evolution. To determine whether the transcriptome dynamics during ESCC recurrence formation are associated with FIR response, an in vitro cell culture model for ESCC radioresistance that mimics the common radiotherapy process in patients with ESCC was established in the present study. High‑throughput sequencing analysis of in vitro cultured ESCC cells was performed using different cumulative irradiation doses, as well as tumor samples from FIR‑treated patients with ESCC before and after the development of radioresistance. Radioresistance‑associated genes and signaling pathways that were aberrantly expressed in radioresistant ESCC cells were identified, including autophagy‑related 9B (regulation of autophagy), DNA damage‑inducible transcript 4, myoglobin and plasminogen activator tissue type, which are associated with response to hypoxia, Bcl2‑binding component 3, tumor protein P63 and interferon γ‑inducible protein 16, which are associated with DNA damage response. The heterogeneity and dynamic gene expression of ESCC cells during acquired radioresistance were further studied in primary (41 single cells), 12 Gy FIR‑treated (87 single cells) and 30 Gy FIR‑treated (89 single cells) cancer cells using a single‑cell RNA sequencing approach. The results of the present study comprehensively characterized the transcriptome dynamics during acquired radioresistance in an in vitro model of ESCC and patient tumor samples at the population and single cell level. Single‑cell RNA sequencing revealed the heterogeneity of irradiated ESCC cells and an increase in the radioresistant ESCC cell subpopulation during acquired radioresistance. Overall, these results are of potential clinical relevance as they identify a number of signaling molecules associated with radioresistance, as well as opportunities for the development of novel therapeutic options for the treatment of ESCC.

Introduction

Esophageal squamous cell carcinoma (ESCC) is a common histological subtype of esophageal cancer that occurs primarily in China (1). Standard therapeutic interventions are surgery, radiotherapy and chemotherapy, which result in a poor 5-year survival rate of <20% (2,3). Although recent advances in clinical diagnosis and therapeutics have improved the survival rate, the local recurrence and refractory radioresistance remain high, leading to a limited chance of cure (4). ESCC cells treated by fractionated irradiation (FIR) may survive due to high adaptive evolution of the cancer cells (5). Therefore, the development of promising therapeutic strategies to increase radiosensitivity and reverse irradiation resistance of ESCC cells is crucial. Understanding the dynamic changes in the gene expression regulation network and heterogeneity of ESCC cells before and after radiotherapy is urgently needed for researchers and clinicians.

Numerous studies have identified differentially expressed genes (DEGs) associated with tumor radioresistance by comparing gene expression in parental cell lines with that in radioresistant sub-clones exhibiting different levels of radiation sensitivity (6,7). In radiosensitive ESCC cells, inactive protein tyrosine kinase 7 (PTK7) overexpression has been demonstrated to promote cell survival and inhibit radiation-induced apoptosis, whereas in radioresistant ESCC cells, knockdown of PTK7 decreased the survival of irradiated cells and increased radiation-induced apoptosis; PTK7 regulates radioresistance through nuclear factor-κB in ESCC (8). In addition, microRNA-205 was upregulated in radioresistant ESCC cells compared with their parental cells by enhancing DNA repair, inhibiting apoptosis and activating epithelial-mesenchymal transition (9). Next-generation sequencing (NGS) has been widely used as a tool for cancer transcriptome analysis (10,11), which can provide detailed transcriptome variations between cancer samples. In addition, using single-cell transcriptome sequencing (scRNA-seq) (12), the heterogeneity of cancer cells with various gene signatures and crucial cancer-associated signaling pathways has been revealed (13). In addition, bioinformatics approaches that identify the regulation networks between genes and phenotypes have revealed functional pathways, as well as indicated the key molecules in complex diseases. The Weighted Gene Co-expression Network Analysis (WGCNA) has been extensively utilized to analyze whole-genome gene expression profiles of microarray and RNA-Seq data (14-17).

To characterize the acquired radioresistance in a population of cancer cells, KYSE-180 ESCC cells were continuously exposed to different cumulative irradiation doses (2, 6, 12, 18, 24 and 30 Gy) to induce radioresistance. This FIR exposure procedure mimics the process of radiotherapy in patients with ESCC. KYSE-180 cells with different FIR dosages were collected and their gene expression profiles were analyzed using whole transcriptome sequencing with non-radiated KYSE-180 cells as controls. It was hypothesized that cancer cells with acquired radioresistance accumulated and that the genes associated with radioresistance dynamically changed during FIR treatment. To identify the subgroups of radioresistant cancer cells and to understand the gene expression dynamics in these subgroups, single-cell transcriptome sequencing was performed on the 0, 12 and 30 Gy FIR-treated KYSE-180 cells to characterize the heterogeneity of irradiated cancer cells. By examining the whole genome transcriptome profiles of in vitro cultured ESCC cells and patient tumor samples before and after acquisition of radioresistance at the population and single cell level, the present study aimed to reveal transcriptomic features that corresponded to FIR treatment and to further clarify the potential mechanisms and significant genes responsible for acquired tumor radioresistance.

Materials and methods

Cell culture and irradiation treatment

The human ESCC cell line KYSE-180 was purchased from the Cell Bank of Type Culture Collection of the Chinese Academy of Sciences and passaged for <4 generations. The cells were cultured in RPMI-1640 medium (Gibco; Thermo Fisher Scientific, Inc.) with 500 ng/ml penicillin-streptomycin and 10% fetal bovine serum at 37°C with 5% CO2.

The irradiation treatment was performed as previously described (18). KYSE-180 cells (3×106 cells/flask) were seeded into 75 cm2 culture flasks. When the cells reached ~70% confluence, the medium was replaced and the cells were irradiated with 2 Gy X-rays using a linear accelerator (Elekta Instrument AB) at an average dose rate of 100 cGy/min, a 20×20 cm field and a source-skin distance of 100 cm. Following irradiation, cells were immediately returned to the incubator. The irradiation (2 Gy) was repeated on days 2 and 3, and the cells were cultured for 4 days to recover. When 90% confluence was reached, cells were trypsinized and sub-cultured into new flasks. These procedures were repeated five times to achieve a total dose of 30 Gy (KYSE-180-30 Gy, RA30) (18). Parental cells used as the irradiation control were trypsinized, counted and passaged under the same conditions without irradiation. When the repeated procedures were completed, cells were trypsinized and washed, and single cells were captured by a micromanipulator for scRNA-seq. The remaining cells were collected for bulk cell RNA-seq.

Colony formation assay

KYSE-180 cells (400 cells/well) were seeded in a 6-well plate and cultured for 24 h. The plates were irradiated with 8 Gy FIR; the control group was not irradiated. Following treatment, all plates were returned to the incubator and cultured for one week. The cells were washed with PBS, fixed with methanol for 10 min at room temperature and stained with 0.1% crystal violet (cat. no. C0121; Beyotime Institute of Biotechnology). Following incubation for 30 min at 37°C, the crystal violet staining was removed, and the plates were washed with deionized water and dried at room temperature. The colonies in each plate were counted to calculate the colony formation rate with light microscope (SZX16, OLYMPUS), magnification 4×.

Immunofluorescence assay of phosphorylated H2A histone family member X (γ-H2AX) expression

Primary KYSE180 cells (2,000 cells/well) in 96-well plates at 1daypost-seeding were exposed to 4 Gy radiation for 2 h, fixed with acetone/methanol (1:1) 15 min, and permeabilized with Triton-X 100 (0.1%) in PBS 30 min. Non-specific binding was blocked by 3% BSA (Sigma-Aldrich; Merck KGaA) in PBS at room temperature for 1 h. Cells were incubated with an anti-γ-H2AX antibody (dilution, 1:100; cat. no. 2595; Cell Signaling Technology, Inc.) at room temperature for 2 h in PBS with 0.1% BSA followed by Alexa Fluor 488-conjugated secondary antibody at room temperature for 1 h (dilution, 1:1,000; cat. no. R32723; Invitrogen; Thermo Fisher Scientific, Inc.) to complete the indirect immunofluorescence procedure. Immunofluorescence images were obtained by confocal laser scanning microscopy (magnification, ×40; 4 fields of view).

Bulk RNA-seq for irradiated ESCC cells and tissues

Bulk KYSE-180 (RA0), KYSE-180-2 Gy (RA2), KYSE-180-6 Gy (RA6), KYSE-180-12 Gy (RA12), KYSE-180-18 Gy (RA18), KYSE-180-24 Gy (RA24) and KYSE-180-30 Gy (RA30) cells (4×106 cells) were collected, when the cells reached the respective doses of irradiation during continuous daily irradiation (two replicates per sample).

Samples of post-radiotherapy and recurrent tumors were obtained from a patient with ESCC. The recruitment period was from May 2016 to May 2018. A male ESCC patient (48 years) was recruited and the following inclusion criteria were applied: (i) FIR treatment; and (ii) cancer tissues collectable via biopsy. Two samples were collected; in May 2016 cancer tissue after radiotherapy and in March 2017 cancer tissue after recurrence. Patient sample collection, patient consent and recruitment followed protocols approved by the Institutional Review Board of The First Affiliated Hospital of Kunming Medical University (Kunming, China).

Total RNA from cancer tissues (homogenized prior to use) and cells were extracted using the RNeasy Mini kit (Qiagen, Inc.) and the quality of the RNA was evaluated using Agilent Bioanalyzer 2100 (Agilent Technologies, Inc.). Sequence libraries were prepared using a TruSeq Stranded mRNA Library Prep kit (Illumina, Inc.) according to the manufacturer's instructions and sequenced.

Single-cell (sc)RNA-seq for irradiated KYSE-180 cells

FIR-treated KYSE-180 cells (RA0, RA12 and RA30) were trypsinized and resuspended in PBS. Cell numbers were determined using a cell counter and adjusted to 1×103 cells/ml. Single cells were placed in 0.2 ml PCR tubes with a micromanipulation system (cat. no. MPC-285; Sutter Instrument Company) (19). A total of 48 cells were prepared for RA0, 96 cells for RA12 and 96 cells for RA30. Single-cell transcriptomes were amplified with SMART-seq2 (20) and cDNAs were purified using Agencourt Ampure XP beads (cat. no. A63881; Beckman Coulter, Inc.). The yield and quality of cDNAs were determined by fluorometric quantitation. DNA >0.5 ng/µl and fragment distributions of 500-7,000 bp were used to evaluate the quality of the single-cell transcriptome amplification. The cDNA (1 ng/cell) was used for library preparation with the Nextera XT DNA Sample Preparation kit (Illumina, Inc.). Sequencing libraries were prepared using the Illumina HiSeq 2500 platform (Illumina, Inc.). A total of 229 FIR-treated single cells (scRA0, scRA12 and scRA30) were sequenced following 3 days of recovery, yielding ~2 million high quality reads for each cell. ScRNA-seq data were uploaded to the National Center for Biotechnology Information reference sequence database with the accession number GSE81812 (https://www.ncbi.nlm.nih.gov/gds/?term=GSE81812).

ScRNA-seq data selection, normalization and DEG identification

Raw RNA-seq reads were trimmed using Trimmomatic 0.35 (21) and mapped to the GENCODE GRCh38 human reference genome (GENCODE v22) using Tophat 2.1.0 software (22). Gene counts were obtained using the Rsubread 1.34.7 and cleaned by scde 2.12.0 (23). Expression levels were calculated by normalized gene counts and variance-stabilizing transformation for clustering by the DESeq2 1.24.0 in Bioconductor 3.9 (24). Relative standard deviations (RSDs) were calculated as standard deviations of each gene divided by average normalized gene count (25). Standard methods were used to evaluate statistical significance. Normalized counts of each gene were used for the analysis of DEGs by Monocle 2.12.0 (26). The sincell 1.16.0 package (27) was used to draw principal component analysis (PCA) plots and the cluster-Profiler 3.12.0 package (28) was used for gene enrichment analysis. To ensure the fidelity of the data, single cells with <1,800 detected genes were excluded. Genes with <10 counts in each cell were excluded and genes that appeared in ≥5 cells were included. All data from scRNA-seq were analyzed and DEGs were identified by comparing scRA0, scRA12 and scRA30.

WGCNA

WGCNA was used for scale-free network topology analysis of RNA-seq data. The WGCNA 1.51 R package (https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/) was used to cluster highly correlated genes and determine clusters in which gene expression was associated with the examined characteristics. An adjacency matrix based on expression correlation was created using a soft threshold procedure to allow scale-free topology. The clusters created by WGCNA are referred to as modules and the minimum number of genes in a module was set to 30. The functional annotation tool DAVID Bioinformatics Resources 6.7 (https://david.ncifcrf.gov/summary.jsp) was used to determine Gene Ontology (GO) terms enriched by the identified genes. DAVID analyses were performed using genes corresponding to significant WGCNA modules.

DAVID analysis

DAVID (http://david.abcc.ncifcrf.gov) functional annotation bioinformatics microarray analysis was used to identify significantly enriched GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms among the genes that were differentially expressed in control samples. Statistically overrepresented GO and KEGG categories with P≤0.05 were considered significant.

Statistical analysis

Data are presented as the mean ± SEM. All of the experiments were independently performed in triplicate. All graphs were plotted and analyzed using GraphPad Prism 5 (GraphPad Software, Inc.) with one-way ANOVA followed by Dunnett's multiple comparisons tests. Correlation coefficients were calculated by Pearson's correlation analysis using R software 3.6.1 (https://www.r-project.org/). P<0.05 was considered to indicate a statistically significant difference.

Results

Induction of radioresistance in human ESCC cells

To characterize the transcriptome dynamics associated with the FIR response, the radiosensitive human ESCC cell line KYSE-180 was exposed to FIR following the procedure described by Jing et al (18). KYSE-180 cells were isolated from well-differentiated human ESCC (JCRB1083), and carried a cyclin D1 amplification and a p53 mutation that commonly occurs in ESCC in the Asian population (29,30). A preliminary study demonstrated that compared with KYSE150 and KYSE30, KYSE180 was more sensitive to FIR (data no shown). The in vitro radioresistance acquisition model was established to simulate a daily FIR treatment process in patients with ESCC. The results demonstrated that daily low-dose irradiation (2 Gy/day) for 3 days did not influence overall cell viability, and a cumulative dose of 12 Gy for 2 weeks caused instability and apoptosis in KYSE-180 cells (Fig. 1A). Following 4 days of recovery, surviving RA12 cells exhibited dramatic cell morphology changes compared with untreated RA0; ~24.4% of RA12 were abnormal and exhibited fusiform morphology (Fig. 1B). When continuous FIRs reaching a 30 Gy cumulative irradiation dose in a period of 5 weeks was applied to these cells, the morphology of the majority of RA30 cancer cells was aberrantly changed (fusiform cells, 60.4%; Fig. 1B). The levels of radioresistance of RA0, RA12 and RA30 cells were determined using colony formation assays combined with DNA damage analysis using γ-H2AX immunostaining. In the colony formation assay, cancer cells were treated with 4 and 8 Gy FIR, the number of colonies was quantified on day 7. RA12 and RA30 displayed increased resistance to FIR by exhibiting significantly more colonies compared with untreated RA0 controls (Fig. 1C). In addition, γ-H2AX immunostaining demonstrated that RA12 and RA30 cells exhibited significantly fewer γ-H2AX-positive cells compared with RA0, indicating that RA12 and RA30 cancer cells rapidly repaired DNA damage following FIR treatment (Fig. 1D). In summary, these results suggested that an in vitro radioresistance acquisition model for ESCC was successfully established.

Figure 1

Cancer radioresistance in response to FIR in vitro. (A) Regimen for induction of radioresistance in ESCC cells. KYSE-180 cells (yellow) were treated with FIR as indicated: For one cycle, cells were treated FIR 2 Gy/day on days 1, 2 and 3. From days 4-7, cells were cultured. Most stressed cells (red) recovered and few cells died. Five cycles of FIR treatment were performed, and populations in each cycle (total radiation 0, 2, 6, 12, 18, 24 and 30 Gy) were collected and sequenced with bulk RNA-seq. Single cells of RA0, RA12 and RA30 (red box) were picked and sequenced with SMART-seq2. Patient tissues were collected and directly sequenced with bulk RNA-seq. (B) Cell morphology of RA0, RA12 and RA30 cells; images were collected with an inverted microscope (magnification, ×20; scale bar, 100 µm). Type I, original KYSE-180 cells; type II, spindle cells (morphology changed; red arrows). (C) Colony formation assay of RA0, RA12 and RA30 cells; relative colony formation was measured with c-value (colony number with FIR/colony number without treatment). (D) Immunofluorescence analysis regarding the ability of DNA damage repair in RA0, RA12 and RA30 cells; cells were stained with anti-γ-H2AX. Images were recorded using a fluorescence microscope (magnification, ×40; scale bar, 50 µm). γ-H2AX positive cells (damaged cells) were defined with ≥5 foci per nuclei. Data are presented as the mean ± SEM; representative of triplicates. ***P<0.01. FIR, fractionated irradiation; ESCC, esophageal squamous cell carcinoma; seq, sequencing; RA, samples with final radiation dose in Gy; γ-H2AX, phosphorylated H2A histone family member X.

Genes involved in the radioresistance of irradiated ESCC cells

To analyze the response of ESCC cells to FIR, whole transcriptome sequencing of irradiated ESCC cells at cumulative irradiation doses of 0, 2, 6, 12, 18, 24 and 30 Gy were conducted (Table SI). The transcriptional profiles of all irradiated ESCC cells were distinguished from that of the non-irradiated RA0 cancer cells. Unsupervised hierarchical clustering of gene expression in the irradiated ESCC cells demonstrated that transcriptional profiles of RA2+RA6, RA12+RA18, and RA24+RA30 were clustered together (Fig. 2A). These results suggested that RA2/RA6, RA12/RA18 and RA24/RA30 cells were at an early, intermediate and late stage of acquired radioresistance, respectively. To assess the transcriptome dynamics of ESCC cells at different stages of irradiation, WGCNA was performed and multiple gene-network modules associated with early, intermediate and late stages of irradiation were obtained (Fig. 2A; Table SII). At the early stage, genes in the 'brown', 'magenta' and 'blue' modules were highly expressed (Fig. 2B). GO term biological process analysis of these genes revealed significant enrichment in the pathways of 'RNA metabolic signaling' (P<0.001), 'cell cycle' (P<0.001) and 'protein translation' (P<0.001) (Fig. 2C and D; Table SIII). These results indicated that during the early stage of irradiation, ESCC cells undergo a rapid cellular response to environmental stress. During the intermediate stage, genes in the 'grey' and 'red' modules were highly expressed (Fig. 2B). Genes in these modules were significantly enriched in the pathways of 'cell adhesion' (P<0.001), 'cell migration' (P<0.001) and 'cell proliferation' (P<0.001) (Fig. 2E and F; Table SIV), which was consistent with the observed morphological changes in the RA12 cells (Fig. 1B). During the late stage of irradiation, genes in the 'cyan', 'light yellow' and 'green-yellow' modules were significantly upregulated (Fig. 2B). Genes in these modules were enriched in cancer radioresistance-associated signaling pathways, including 'response to hypoxia' (P=0.0015) (31), 'regulation of autophagy' (P=0.0038) (32) and 'DNA damage response' (P=0.016) (33) (Fig. 2G; Table SIII). In addition, KEGG enrichment analysis demonstrated that 'mTOR signaling pathway' (P=0.057) was also activated in the RA24 and RA30 cells (Fig. 2H; Table SIV) (34,35). In summary, transcriptome dynamics revealed that ESCC cells respond to FIR rapidly with increased RNA metabolism at the beginning of radiotherapy, followed by upregulated cell adhesion and finally acquired radioresistance with the activation of multiple known radioresistance-associated signaling pathways during continuous irradiation.

Figure 2

WGCNA revealed gene-network modules responsible for different stages of fractionated irradiation. (A) WGCNA dendrogram indicating expression of gene modules in irradiated esophageal squamous cell carcinoma cells; early stage of irradiationwithRA2 and RA6 clustered, intermediate stage of irradiation with RA12 and RA18 clustered, and late stage of irradiation with RA24 and RA30 clustered. (B) Modules-trait relationships in irradiated ESCC cells; module color codes are preserved; modules expressed in different stages of FIR (red boxes). (C) GO biological process and (D) KEGG analysis of the early stage of irradiation; 'brown', 'magenta' and 'blue' were identified and analyzed with DAVID (P<0.01). WGCNA revealed gene-network modules responsible for different stages of fractionated irradiation. (E) GO biological process and (F) KEGG analysis of the intermediate stage of irradiation; 'grey60' and 'red' were identified and analyzed with DAVID (P<0.01). (G) GO biological process and (H) KEGG analysis of the late stage irradiation; 'cyan', 'lightyellow' and 'greenyellow' were identified and analyzed with DAVID (P<0.01). WGCNA, Weighted Gene Co-expression Network Analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; RA, samples with final radiation dose in Gy.

The modules that reflected unique gene expression dynamics across irradiated ESCC cells were further analyzed. A total of four modules of interest were identified: 'Turquoise', 'brown', 'magenta' and 'cyan' (Fig. 3A and B). In addition, hub-gene-network analysis of these modules revealed highly connected genes in each module, through which key controlling genes in the modular network were identified (Figs. 3C and D and S1). The 'turquoise' module contained a significant number of known cancer-related genes, such as WNT5A (36,37), transmembrane protease serine 2 (TMPRSS2) (38), ATPase copper-transporting α (ATP7A) (39) and others. The 'turquoise' module gene expression significantly decreased in irradiated ESCC cells (Fig. S1), suggesting that these genes were downregulated by radiotherapy. The 'brown' module, which contained genes that were upregulated in the early stage of irradiation and declined during further irradiation treatment, was significantly enriched in genes involved in RNA metabolism (Fig. S1). The expression of genes in the 'magenta' module, which were abundantly expressed in untreated RA0 and early-stage RA2 and RA6 cells, continuously decreased following cumulative irradiation treatment. Thus, the expression of the genes in the 'magenta' module was negatively associated with the acquisition of radioresistance in ESCC and may serve functional roles in the development of cancer cell radioresistance. GO term analysis demonstrated that the genes in the 'magenta' module were significantly enriched in biological processes involved in 'WNT signaling' (P=0.043) and 'chromosome organization' (P=0.0032) (Table SIII). The network of the top hub-genes in the 'magenta' module is presented in Fig. S1. By contrast, genes in the 'cyan' module, which were aberrantly overexpressed in RA30 cells compared with untreated and irradiated ESCC cells at the early and intermediate stages, may be associated with the acquisition of radioresistance in ESCC. GO term analysis demonstrated that the genes in the 'cyan' module were significantly enriched in the biological processes involved in 'response to hypoxia' (P<0.0001) and 'DNA damage' (P=0.006) (Fig. 3C; Table SIII). Previously reported radioresistance-related genes in the 'cyan' module were potassium calcium-activated channel subfamily N member 4 (KCNN4) (40), spermidine/spermine N1-acetyltransferase 1 (SAT1) (41) and metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) (42,43) (Fig. 3D). The hub-gene autophagy-related 9B (ATG9B) in the 'cyan' module may serve important roles in ESCC radioresistance through autophagy. DNA damage-inducible transcript 4 (DDIT4), myoglobin (MB) and plasminogen activator tissue type (PLAT) are associated with hypoxia (44), which may also be associated with radioresistance. Pseudogenes of the WASH protein family that recruit and activate the Arp2/3 complex to induce actin polymerization, which have a key role in the fission of tubules that serve as transport intermediates during endosome sorting (45), WASH2P, WASH3P, WASH5P and WASH6P may be related to radioresistance. Overall, these results demonstrated that transcriptome dynamics were associated with cell adhesion, autophagy, hypoxia and DNA damage during irradiation treatment and revealed a number of genes that may serve crucial roles in the development of radioresistance.

Figure 3

Hub-genes network of cancer radioresistance in ESCC cells and patient. (A) Gene expression per modules in irradiated ESCC cells; 'turquoise', 'magenta' (early stage), 'brown' (intermediate stage) and 'cyan' (late stage) were assessed further. (B) Eigengene connectivity of 'turquoise', 'magenta', 'brown' and 'cyan'. (C) GO analysis of genes in 'cyan' (P<0.01); putative radioresistant signaling pathways in blue. (D) Hub-gene network of 'cyan'; genes associated with radioresistance in red and genes associated with hypoxia in blue, previously reported genes that related with cancer therapeutic resistance were marked with boxes. (E) GO analysis of genes upregulated in relapse cancer tissues compared primary cancer tissues (P<0.01); putative radioresistant signaling pathways in dark blue. (F) Venn diagram of genes in 'cyan' and upregulated genes in ESCC relapse tissue; red highlights genes associated with cancer resistance. ESCC, esophageal squamous cell carcinoma; RA, samples with final radiation dose in Gy; GO, gene ontology.

Radioresistance genes in tumor samples from FIR-treated patients with ESCC

To validate the radioresistant genes identified in the in vitro experiments in the present study, paired tumor tissue samples were collected from one radiotherapy-treated and relapsed patient with ESCC. Comparison between the gene expression profiles of the primary and relapsed tumor samples revealed that 296 genes were upregulated in the relapsed tumor samples (Table SV). GO term biological process analysis of these genes demonstrated significant enrichment in the pathways associated with 'cell response to stress' (P=0.006), 'cell cycle' (P=0.006), 'DNA repair' (P=0.009) and 'hypoxia' (P=0.012). The result obtained using patient samples were consistent with the in vitro results (Fig. 3E; Table SVI). To examine whether the relapsed tumor was derived from residual radioresistant cancer cells, the 296 upregulated genes were compared with the previously identified radioresistance-associated genes in the 'cyan' module. A total of 14 overlapping genes were identified, including certain known resistance-related genes tyrosine kinase 2 (TYK2) (46), growth factor receptor-bound protein 7 (GRB7) (47) and insulin-like growth factor-binding protein 3 (IGFBP3) (48) (Fig. 3F). These results demonstrated a conserved common molecular mechanism in radioresistance between the in vitro ESCC model and the patient samples, indicating that radioresistance genes ('cyan' module) identified in the irradiated cell line model serve crucial roles in acquired radioresistance.

Radioresistance dynamics at the single-cell level

To determine the mechanism of the acquisition of radioresistance in ESCC, whole transcriptome analysis was performed at the single-cell level in RA0, RA12 and RA30 cancer cells. ScRNA-seq reads were mapped to the reference human genome with an average of 65.4% reads mapped within the genome. To test the reproducibility of data, the sequencing was repeated using the same templates; a strong correlation (r=0.96) was obtained, which confirmed the reliability of the sequencing data (Fig. S2A). To obtain accurate estimates of the gene expression levels, the sequencing depth was normalized to ensure that all cells exhibited approximately the same median read depth (Fig. S2B). A total of 217 cells passed this filter, which resulted in 41 qualified RA0 single cells (scRA0), 87 qualified RA12 single cells (scRA12) and 89 qualified RA30 single cells (scRA30). To examine the gene expression levels in individual cells, the RSD of normalized gene counts was used to estimate the dispersion. Genes with higher expression levels exhibited lower degrees of dispersion (Fig. S2C). Genes with <4 normalized read counts (≤2 log2 value) were removed as they provided an unreliable measure of gene expression levels (25). The reproducibility of the Illumina sequencing platform was guaranteed by a high correlation coefficient (r=0.98; Fig. S2D). Overall, the results demonstrated that the quality of the scRNA-seq data was satisfactory for further analysis.

One-way ANOVA was used to analyze the differences between the RNA-seq data of each group. A total of 3,174 DEGs were identified for further analysis. PCA demonstrated thatscRA0 was distinct from scRA12 and scRA30, whereas scRA12 and scRA30 were partially separated, and a number of single cells from scRA12 and scRA30 were clustered together (Figs. 4A, S3A and B). The scRNA-seq data were subjected to WGCNA and multiple gene-network modules were obtained (Table SVII). Unsupervised hierarchical clustering of 217 single cell transcriptome results revealed three major groups (I-III) and three subgroups in group II: (i) Group I, scR0; (ii) group II, scR1, scR2, and scR3; and (iii) Group III, scR4 (Fig. 4B). The distribution of single cells in each group was further examined; the results demonstrated that scR1, scR2, scR3 and scR4 contained different percentages of scRA12 and scRA30. All cells in group I were untreated RA0 cells. From scR1 to scR4, the number of scRA12 gradually decreased, whereas the number of scRA30 gradually increased (Fig. 4C). These results suggested that the cells in different groups obtained various levels of radioresistance, and the ESCC cells in the scR4 subgroup may have developed the highest expression levels of radioresistance-associated genes compared with the cells in the other subgroups (Fig. 4C; Table SVIII).

Figure 4

Single-cell transcriptome sequencing reveals radioresistance. (A) PCA of 217 single cells with differentially expressed genes. (B) WGCNA dendrogram indicating expression of gene modules in single cells; based on the unsupervised hierarchical clustering, with group I including scR0, group II including scR1, scR2, and scR3, and group III including scR4. (C) Distribution of scRA0, scR12 and scR30 in clustered groups. (D) Modules-trait relationships in single-cell groups; color code of modules is preserved. (E) Hub-gene network of scBlue; cancer resistance associated genes highlighted with red box. (F) Hub-gene network of scYellow; cancer resistance associated genes highlighted with red box. WGCNA, Weighted Gene Co-expression Network Analysis; PCA, principal component analysis; RA, samples with final radiation dose in Gy; scR, single radioresistant cell.

The scRA12 and scRA30 were divided into two subgroups (scRA12, scRA12-1 and scRA12-2; and scRA30, scRA30-1 and scRA30-2), and the module-trait relationships were determined (Fig. 4D). The results demonstrated that the genes in the 'yellow', 'blue' and 'brown' modules were highly expressed in single cells from groups II and III, whereas genes in the 'black', 'pink', 'green', 'magenta' and 'turquoise' modules were highly expressed in single cells from group I. Hub-gene network analysis in each module revealed the key controlling genes in the modular network (Figs. 4E, F and S4). The results demonstrated that the hub-genes in the 'blue' module of single-cell data (scBlue) contained multiple known cancer resistant-related genes, including reactive oxygen species modulator 1 (ROMO1) (49), voltage-dependent anion channel 2 (VDAC2) (50), Dickkopf WNT signaling pathway inhibitor (DKK1) (49,51) and paxillin (PXN) (52-54), which have been reported in previous studies. The hub-genes in the 'yellow' module (scYellow) contained neurotrophin 3 (NTF3) (55), FA complementation group E (FANCE) (56) and mutS homolog 3 (MSH3) (57), which have been previously demonstrated to be associated with apoptosis or DNA damage repair processes. Other genes in this module are also reported to associate with cancer resistance, such as ETS transcription factor (ELK1) (58), PBX homeobox 1 (PBX1) (59), GC-rich sequence DNA-binding factor 2 (GCFC2) (60) and mitochondrial leucyl-tRNA synthetase 2 (LARS2) (61). These results revealed the heterogeneity and dynamic gene expression changes of ESCC cells during the acquisition of radioresistance.

Discussion

The present study successfully established an in vitro radioresis-tance acquisition model using radiosensitive ESCC KYSE-180 cells through consistent low dose FIR treatment. The results demonstrated that this in vitro model recapitulated the acquisition process of cancer cell radioresistance following radiotherapy in patients with ESCC. Throughout the irradiation process, multiple ESCC cells receiving different doses of irradiation were collected and, for the first time, the transcriptome profiles of an in vitro ESCC model at the population and single-cell level were analyzed. The results revealed that irradiation-induced dynamics of gene expression and modules were associated with the stage of radioresistance. The genes in the 'cyan' module from the population RNA-seq were specifically enriched in ESCC cells at the late stage of irradiation. These genes were involved in signaling pathways associated with autophagy, hypoxia and DNA damage repair. Considering the crucial roles of autophagy (62) and hypoxia (63) in radioresistance development, the transcriptome analysis results in the present study suggested that MALAT1-ATG9B (autophagy) and DDIT4-MB-PLAT (hypoxia) may induce radioresistance in cancer cells.

The dynamic changes in gene expression profiles were compared between the FIR-induced in vitro radioresistance model and tumor samples from one patient with ESCC undergoing radiotherapy before and after relapse. Although only 5% (14/296) of the genes were consistent between radioresistant ESCC cells and patient-derived samples, two signaling pathways, DNA damage and hypoxia, were conserved in the model and the patient. More patient samples are needed to validate these results; other mechanisms may contribute to radioresistance in other patients in addition to these two signaling pathways. Among the overlapping genes, TYK2 and IGFBP3, which were previously reported as chemoresistance-associated genes (46), were signifi-cantly upregulated in RA30 ESCC cells and in the relapsed tumor sample. GRB7, which has been reported to be associated with Erb-B2 receptor tyrosine kinase 2(HER2) amplification (64), HER2 signaling inhibition (65) and an increased risk of recurrence in triple negative breast cancer (47), may functionally regulate radioresistance in ESCC. These results suggested that the FIR-induced in vitro radioresistance model shared a similar molecular regulation network involving radioresistance; thus, the in vitro model established in the present study may bean efficient tool to study clinical cancer radioresistance.

Single-cell transcriptome analysis of the FIR-induced in vitro radioresistance model in the present study revealed cellular heterogeneity in irradiated ESCC cells. During the process of irradiation, radioresistance is achieved by a continuous dynamic selection process. Different single ESCC cells may activate diverse radioresistance-related signaling pathways and acquire various levels of radioresistance. Among all irradiated single cancer cells (scRA12 and scRA30), four subgroups were identified. Single cells in the scR1 subgroup, in which the expression profile was similar to that in the original single ESCC cells (scRA0 or scR0), may still be sensitive to irradiation. The cells in subgroups scR2, scR3 and scR4, which exhibited significant differences in gene expression profiles compared with the cells in the scR0 subgroup, may acquire varying levels of radioresistance. Analysis of the percentage of scRA12 and scRA30 cells in each subgroup revealed that the cells in scR2, scR3 and scR4 may obtain increased sensitivity to radioresistance, and the cells in scR4 acquired the highest level of radioresistance. The scRA12 cells were distributed among the scR1, scR2, scR3 subgroups, and one cell in the scR4 subgroup, suggesting that cancer radioresistance may occur in the stage of RA12. By contrast, 10% (9/89) of scRA30 cells were in the scR1 subgroup and 18% (16/89) of scRA30 cells were in the scR4 subgroup, which indicated a heterogeneous cell population of scRA30 cells.

In the present study, WGCNA revealed that the scBlue, scYellow (upregulated) and scTurquoise (downregulated) modules significantly represented the molecular features of each subgroup. Gene network analyses of these modules identified the hub-genes involved in cancer radioresistance. The scBlue and scYellow modules exhibited low consistence with the bulk 'cyan' module and samples from the patient with ESCC (Fig. S5), which may be due to a mixed population of radioresistant cancer cells in RA30. However, a number of overlapping genes were identified among these samples, including GRB7 (64,65). Other overlapping genes in subgroup scR4 and relapsed tumor samples may be of importance in acquired radioresistance. Further validation and investigation, including extensive single-cell transcriptome analysis in radioresistant ESCC tissues or ESCC derived PDX models, will be required to determine the precise molecular pathways of acquired radioresistance. The results of the present study supported the hypothesis that single-cell transcriptome analysis of cancer cell lines and patient tumor samples may bring invaluable insight into the development of acquired radioresistance from a specific population of cancer cells. Overall, these results are of potential clinical relevance, as they revealed genes and signaling pathways involved in radioresistance and identified opportunities for developing novel therapeutic options for ESCC.