Investigating Molecular Determinants of Cancer Cell Resistance to Ionizing Radiation Through an Integrative Bioinformatics Approach.

Eradication of cancer cells through exposure to high doses of ionizing radiation (IR) is a widely used therapeutic strategy in the clinical setting. However, in many cases, cancer cells can develop remarkable resistance to radiation. Radioresistance represents a prominent obstacle in the effective treatment of cancer. Therefore, elucidation of the molecular mechanisms and pathways related to radioresistance in cancer cells is of paramount importance. In the present study, an integrative bioinformatics approach was applied to three publicly available RNA sequencing and microarray transcriptome datasets of human cancer cells of different tissue origins treated with ionizing radiation. These data were investigated in order to identify genes with a significantly altered expression between radioresistant and corresponding radiosensitive cancer cells. Through rigorous statistical and biological analyses, 36 genes were identified as potential biomarkers of radioresistance. These genes, which are primarily implicated in DNA damage repair, oxidative stress, cell pro-survival, and apoptotic pathways, could serve as potential diagnostic/prognostic markers cancer cell resistance to radiation treatment, as well as for therapy outcome and cancer patient survival. In addition, our findings could be potentially utilized in the laboratory and clinical setting for enhancing cancer cell susceptibility to radiation therapy protocols.

Introduction

Radiation therapy or radiotherapy (RT) represents one of the optimal, most widely used modalities in the treatment of multiple cancers, either alone or combined with other curative anti-cancer modalities like chemotherapy (Delaney et al., 2005; Begg et al., 2011) or immunotherapy (Tang et al., 2014; Schoenhals et al., 2016). It is estimated that approximately 50% of all cancer patients worldwide undergo radiotherapy throughout their illness trajectory (Baskar et al., 2012).

Advances in radiotherapy contribute greatly to cancer patients' improvement of overall survival and quality of life (Baskar et al., 2012). The aim of radiotherapeutic regimens is to specifically and efficiently sensitize cancer cells to IR in order to eliminate them and prevent cancer recurrence and relapse, minimizing at the same time the adverse effects of radiation on healthy tissue. RT affects cancer cells either directly, by inducing genomic (DNA) lesions, or indirectly, through the generation of DNA damaging intermediates through the interaction with water, like reactive oxygen/nitrogen species (ROS/RNS) and free radicals (e.g., hydrogen ion, hydroxide, etc.) (Mikkelsen and Wardman, 2003; Yamamori et al., 2012).

However, cancer cells have the capacity to develop incredible tolerability and resistance to RT, thereby evading death. Radioresistance represents a major limiting factor in the effective treatment of different types of cancers. The response of tumor cells to radiation depends both on the resistance mechanisms of the cells and also on the accelerated repopulation of the tumor bulk by cells that have developed further radioresistance (Pavlopoulou et al., 2016, 2017). As noted in previous studies, the genes that are differentially expressed (either up- or down-regulated) between radioresistant (RR) and radiosensitive (RR) cancer cells are generally implicated in DNA damage response and repair (DDR/R) pathways, apoptosis, hypoxia, or response to oxidative stress, etc. (Pavlopoulou et al., 2016, 2017). The complexity of radiation resistance mechanisms suggests the involvement of different and diverse biological mechanisms.

During the last decade, the advances in high-throughput (HTP) “omics” technologies (e.g., RNA-Seq and microarrays) enabled the generation of an enormous amount of gene expression data. Data produced with HTP technologies are stored in international public repositories such as NCBI's GEO (Gene Expression Omnibus) (https://www.ncbi.nlm.nih.gov/gds/) (Barrett et al., 2013; Clough and Barrett, 2016). GEO DataSets contains both original records and curated datasets.

Accumulated knowledge over years of research on the biological effects of radiation points toward the development of holistic approaches to “big data” analysis by employing systems biology methodologies (Unger, 2014; Beheshti et al., 2019; Spratt and Speers, 2019; Kanakoglou et al., 2020). Herein, we employed a rigorous systems biology approach to unravel the molecular determinants of resistance of cancer cells to IR, based solely on HTP data. To this end, publicly available transcriptome datasets relevant to cancer cell response to radiation were retrieved from GEO, and specifically, cancer cell lines that displayed enhanced resistance to radiation. Statistical analyses were carried out to identify the differentially expressed genes (DEGs) between radioresistant and radiosensitive tumor cells. Furthermore, functional annotation of genes allowed us to identify specific biological pathways implicated in cancer cell resistance to radiation. Our findings could be applied in the laboratory and clinical setting as biomarkers for the design of targeted and personalized radiotherapy regimens in order to effectively sensitize cancer cells to radiation, enhance tumor control and thereby minimize tumor recurrence and metastasis.

Methods

Data Retrieval

The public repository NCBI GEO DataSets was searched extensively for gene expression datasets using relevant keywords: (“radiation therapy” or “radiotherapy”) and (“cancer” or “tumor”) and (“resistance” or “tolerability”) and (“sensitivity” or “responsive”) and (“human” or “homo sapiens”). A total of three eligible datasets were selected:

The GEO series GSE97543 (Emons et al., 2017) (Supplementary Table 1) includes global gene expression by microarray of both wild-type and radioresistant Dukes' type C colorectal adenocarcinoma (COAD) cell lines that were either non-irradiated or irradiated repeatedly with 2 Gray (Gy) of X-rays in order to acquire a radioresistant phenotype. The Agilent-026652 whole human genome microarray 4x44K, GPL13497 platform was used.

In GSE13280 (Marston et al., 2009) (Supplementary Table 1), genome-wide gene expression by microarray was performed of cell lines derived from pediatric B-precursor acute lymphoblastic leukemia (ALL) after 8 h in vitro exposure to 5 Gy IR. This dataset contains cell lines both resistant and responsive to radiation. The development of resistance and responsiveness to IR was assessed by measuring apoptosis in cells. The Affymetrix human genome U133A array, GPL96 platform was employed.

In GSE120798 (Gray et al., 2019) (Supplementary Table 1), three novel radioresistant breast cancer cell lines were established by exposing the corresponding parental cell lines MDA-MB-231 (metastatic mammary adenocarcinoma), MCF-7 (breast adenocarcinoma), and ZR-751 (luminal breast cancer) to increasing doses of X-rays for 2 and 8 h. Genome-wide gene expression analysis of both the parental and radioresistant cell lines (i.e., MDA/MDAR, MCF/MCFR, ZR/ZRR) was performed using high throughput sequencing. The NextSeq 550 (Homo sapiens), GPL21697 platform was used.

Microarray-Based Transcriptomic Data Analysis

For each microarray study, the gene expression data that represent the gene expression summary for every probe and every sample were recorded. In microarrays, many probes can map to the same Gene Symbol for various reasons, and, conversely, a probe may also map to more than one Gene Symbol if the probe sequence is not specific enough. A simple approach would be to use only the probes with one-to-one mapping for further analysis; however, this approach results in the loss of important information. To conduct an analysis based on genes and not probes, the probe identifiers were firstly converted into gene identifiers, according to Ramasamy et al. (2008) guidelines. To this end, GPL files that contained information about the gene symbols that correspond to probe IDs were used in order to resolve the “many-to-many” relationships between probes and genes by averaging the expression profiles for genes with more than one probe (Ramasamy et al., 2008). We identified the Gene Symbols with the usage of the HUGO Gene Nomenclature Committee (Braschi et al., 2019) and the National Center for Biotechnology Information (NCBI) GENE (Sayers et al., 2020).

The two-sample t-test was employed to identify genes differentially expressed between the case (RR) and control (RS) groups. However, a disadvantage of the t-test in the analysis of microarray data is that if most of the experiments in a given study contain a relatively small number of samples per group, the assumption of normality is untenable. To resolve this, the statistical method t-test with bootstrap was used (Efron and Tibshirani, 1993). Bootstrap provides an ideal method to generate accurate estimates of the standard errors when no formula for the sampling distribution is available or when available formulas make inappropriate assumptions (e.g., small sample size, non-normal distribution). In this study, bootstrap analysis was conducted with 1,000 replicates, a relatively high number, in order to generate accurate estimates of the standard errors.

A typical microarray experiment measures the expression of several thousand genes simultaneously across different conditions. When investigating for potential DEGs between two conditions, each gene is treated independently, and the t-test is performed on each gene separately. The incidence of false positives (i.e., genes falsely declared as DEGs) is proportional to the number of tests performed and the critical significance level (p-value cut-off). In order to account for multiple comparisons, a correction method proposed by Benjamini and Hochberg (1995) which controls False Discovery Rate (FDR) was applied. FDR-controlling procedures have greater power (i.e., they can discover more statistically significant differences), at the expense of increased Type I error rate. Genes with adjusted p-value (or q-value) less or equal to 0.05 were considered as statistically significant in this study. For all statistical analyses, the Stata 13 statistical software package (StataCorp, 2013) was used. For the creation of heatmaps from microarray data, the average linkage clustering with Euclidean distance clustering method implemented in Heatmapper (http://heatmapper.ca/) was utilized.

RNA-Seq-Based Transcriptomic Data Analysis

For the RNA-Seq based transcriptome analysis, the following pipeline was utilized. The FASTQ files were extracted from the respective Sequence Read Archive (SRA) files containing raw RNA-Seq reads by using the SRA Tool Kit v.2.9.0 (Alnasir and Shanahan, 2015) with the “fastq-dump –gzip –skip-technical –readids –dumpbase –clip –split-3” command. The raw RNA-Seq reads in FASTQ files were aligned to the human reference genome GRCh38 (Ensembl version 97) by employing the splice junction aligner HISAT2 v.2.1.0 (Kim et al., 2015) with “hisat2 -p -dta -x {input.index} -U {input.fq} -S {out.sam}” parameters. The generated SAM file was converted to the respective binary BAM file by using SAM Tools v.1.9.0 (Li et al., 2009) with “samtools sort -@ 10 -o {output.bam} {input.sam}” commands. String Tie v.1.3.5 (Pertea et al., 2015) was utilized with “stringtie -e -B -p -G {input.gtf} -A {output.tab} -o {output.gtf} -l {input.label}{input.bam}” parameters for the measurement of gene expression levels. The reconstructed transcripts and transcript abundances were reported in the output GTF file. In order to detect differentially expressed genes between RR and RS samples, we utilized the EdgeR package v3.28.0 (Robinson et al., 2010) of the R statistical computation environment v.3.6.1 (https://www.r-project.org). We firstly applied the trimmed mean of M-values (TMM) normalization (Robinson and Oshlack, 2010) implemented in EdgeR to the count data and we employed generalized linear models with “cell lines,” “resistance,” and “time” as factors. Then, estimating dispersion was computed with the estimateDisp function, and differential expression analysis between the two RNA-Seq groups (RR and RS) was performed using the glmFit and glmLRT functions of the EdgeR package v3.28.0 (Robinson et al., 2010) of the R statistical computation environment v.3.6.1 (https://www.r-project.org). In order to detect statistically significant differentially expressed genes, the threshold for the absolute log2-fold change was set at two (|log2FC≥2|), and for the FDR (Benjamini and Hochberg, 1995)-corrected p-value at 0.05. All statistical calculations for the RNASeq data were performed by using the R software environment. The pheatmap package of R (https://CRAN.R-project.org/package=pheatmap) was utilized to generate a heatmap from RNA-Seq data.

Pathway Enrichment Analysis

To assign biological role(s) to the genes under study that are associated with biological pathways, gene set enrichment analysis (GSEA), or functional enrichment analysis, was conducted. GSEA is a method to identify biological processes or pathways that are over-represented in a large set of genes. To this end, WebGestalt (WEB-based GEne SeT AnaLysis Toolkit) (Zhang et al., 2005; Liao et al., 2019) was employed to identify statistically significant over-represented WikiPathways (Kutmon et al., 2016) cancer terms in the sets of genes; the threshold for the FDR-corrected p-value was set at 10−3, and hypergeometric distribution analysis was used.

Functional Interactions Networks

The associations among the molecules under study were investigated and visualized with the usage of STRING (Search Tool for Retrieval of Interacting Genes/Proteins) v11.0 (Szklarczyk et al., 2019), a database of either known or predicted, direct or indirect, gene/protein associations derived from diverse resources. The highest confidence interaction score (≥0.9) was chosen to display the associations amongst genes/proteins in the generated network.

Survival Analysis

The prognostic potential of the 36 radiogenes was investigated in three types of cancers, namely, breast invasive carcinoma (BRCA), colon adenocarcinoma (COAD), and acute myeloid leukemia (LAML), a type of haematologic cancer like ALL, through the web-based tool GEPIA (Gene Expression Profiling Interactive Analysis) (Tang et al., 2017) version 2 (http://gepia2.cancer-pku.cn/#index), based on data acquired from The Cancer Genome Atlas (TCGA) (https://tcga-data.nci.nih.gov). The cancer patient cohort is divided into the high-risk and low-risk categories; the cut-offs for low and high gene expression level patient cohorts were set at 50%.

Tissue-Wise Differential Gene Expression Analysis

GEPIA2 (Tang et al., 2017), which contains gene expression data from cancer and corresponding normal tissues from TCGA (https://tcga-data.nci.nih.gov) and the Genotype-Tissue Expression (GTEx) (https://gtexportal.org/home/), respectively, was used to investigate the differential expression patterns of the radiogenes under study in BRCA, COAD, and LAML cancer-normal tissue; |logFC| ≥ 2 and FDR-corrected p-value ≤ 0.05.

Results

The overall procedure for data collection and analysis followed in the present study is described illustratively in Figure 1.

Figure 1

Flowchart diagram of the overall methodology employed in this study. Three eligible transcriptome datasets were retrieved from NCBI's GEO. Genes differentially expressed (DEGs) between the radioresistant and radiosensitive cancer cells were identified per dataset. The DEGs of the three analyzed datasets were integrated and compared with an earlier own study on cancer cell radioresistance-related genes/proteins derived from literature research. Functional enrichment analysis of the common genes was performed to obtain the so-called “radiogenes”.

Identification of Differential Expression Patterns in Radioresistance vs. Radiosensitive Cancer Cells

The two techniques for transcriptome profiling, RNA-Seq and microarray, have large inherent differences. RNA-Seq is considered “superior” since it allows the detection of low abundance transcripts and novel transcript isoforms (Marioni et al., 2008). For this reason, we applied different statistical methods for RNA-Seq (Li, 2019) and microarray (Kontou et al., 2018) data processing and analysis.

The number of differentially expressed genes (DEGs) found between radioresistant and radiosensitive cancer cells (Figure 2) for each dataset is 6372 (GSE120798), 782 (GSE13280), and 541 (GSE97543), respectively (Supplementary Table 2). The pathways over-represented in the DEGs of the three datasets are related primarily to DDR/R and cell survival (Figure 3).

Figure 2

Heatmap representing color-coded expression levels of DEGs for the datasets (A) GSE120798 (breast cancer), (B) GSE13280 (ALL), (C) GSE97543 (COAD); columns correspond to samples and rows correspond to genes. RR, Radioresistant; RS, Radiosensitive.

Figure 3

Distribution of the over-represented cancer-related WikiPathways in the DEGs of each transcriptome dataset. The enriched pathways are indicated by gray.

Differential gene expression analysis was performed for the three parental breast cancer cell lines in GSE120798 (Gray et al., 2019) against their radioresistant derivatives, to identify those genes that are significantly dysregulated in response to radiation stress. A total of 11 samples were compared; one biological replicate for each condition (Supplementary Table 1). The number of detected DEGs (Gray et al., 2019) is remarkably higher as compared to the ones derived from microarray gene expression data (Supplementary Table 2). This discrepancy is likely due to the ability of RNA-Seq to detect and quantify, even rare, transcripts without a priori knowledge of a given gene (Metzker, 2010). Accordingly, the number of enriched pathways is also greater in this dataset (Figure 3). Three DDR/R pathways were found to be enriched in the DEGs, including the generic “DNA Damage Response” pathway. Two pathways are specifically implicated in IR-induced DNA damage response, namely “DNA IR-Double Strand Breaks (DSBs) and cellular response via ATM,” and “DNA IR-damage and cellular response via ATR.” Ataxia telangiectasia mutated (ATM) and ataxia telangiectasia and Rad3-related (ATR) proteins are evolutionarily conserved proteins that have a critical regulatory role in DDR and maintenance of genome integrity (Marechal and Zou, 2013; Awasthi et al., 2015). Furthermore, the “PI3K-AKT-mTOR signaling pathway and therapeutic opportunities” was shown to be over-represented in the set of DEGs; increased activity of PI3K in the radioresistant cells of this transcriptome dataset was also observed by Gray et al. (2019). Phosphatidylinositol-3-kinase (PI3K)/AKT/mammalian target of rapamycin (PI3K/AKT/mTOR) signaling is critical to many aspects of tumor cell growth and survival (Porta et al., 2014) and therefore could be likely involved in the survival of irradiated cancer cells.

Differential expression analysis was also carried out of the irradiated COAD cells in the GSE97543 dataset (Emons et al., 2017). A total of six samples were compared; three biological replicates for the wild-type (radiosensitive) cell lines and three for the radioresistant cells (Supplementary Table 1). The notch signaling pathway is significantly over-represented (Figure 3) in the DEGs of this dataset (Supplementary Table 2). Notch signaling is suggested to confer a selective survival advantage on tumors (Capaccione and Pine, 2013). Hence, the Notch network could be implicated in the resistance and survival of the COAD cells to irradiation. Regarding the over-represented DDR/R pathways, one pathway is related to processing IR-induced DNA lesions through ATR signaling and the other to mismatch repair (MMR), which is responsible for detecting and repairing mismatched nucleotides (Iyer et al., 2006; Larrea et al., 2010). MMR reaction is initiated by binding of the MSH2 (MutS homolog 2)/MSH6 heterodimer to the mismatched DNA; both MSH2 and MSH6 were found differentially expressed in the radioresistant COAD cells. The MutS homologs MSH2 and MSH6 form a heterodimer that binds to short insertion/deletion DNA mispairs (Habraken et al., 1996; Edelbrock et al., 2013). The MMR proteins also function in signaling DNA damage (Duckett et al., 1996b; Modrich, 1997). Earlier studies have shown that MSH2 is also involved in the processing of the biologically significant clustered DNA damages, as well as the execution of apoptosis induced by IR (Holt et al., 2009).

In GSE13280 (Marston et al., 2009), the gene expression profiles of radioresistant and radiosensitive ALL cell lines (eleven replicates for each condition) were compared (Supplementary Table 1) for the identification of DEGs (Supplementary Table 2). The enriched transforming growth factor beta (TGFB) signaling pathway (Figure 3), like the Notch network (Capaccione and Pine, 2013), is implicated in several aspects of cancer initiation, promotion, and progression (Syed, 2016). Hence, the Notch- and TGFB-mediated signaling pathways might render ALL cells less vulnerable to IR-induced apoptosis by exerting their cellular pro-survival effect.

Identification of Cancer Cell Radioresistance–Related Genes

The procedure we followed to identify an optimal number of cancer cell radioresistance–related genes is illustrated in Supplementary Figure 1. A list of 175 bio molecules (Supplementary Table 3) was proposed in a previous study by Pavlopoulou et al. (2017) to be implicated in tumor cell radioresistance. These genes/proteins were manually collected through a comprehensive and thorough literature search. The DEGs identified in each of the three transcriptome datasets were merged; a list of 7,185 genes was compiled (Supplementary Table 3). Those genes were compared to the literature-derived molecules, and 88 genes were found in common (Supplementary Table 3). In order to identify an optimal number of genes implicated in radioresistance, we performed functional enrichment analysis of these genes. The pathways over-represented in the 88 genes (Figure 4) are related to DDR/R, similar to those detected in the DEGs of the individual datasets (Figure 3), as well as to apoptosis. Signaling pathways mediated by the cardinal tumor suppressor TP53 are also related to radioresistance (Figure 4). Collectively, 36 genes were found to be implicated in cancer-associated biological pathways listed in Table 1, 26 of those up-regulated and 10 down-regulated in RR cancer cells. The products of the 36 genes also form a highly connected network (Figure 5), with high-confidence interactions, suggesting that these proteins associate, either functionally or physically, to confer cellular resistance to IR. Therefore, we propose 36 interconnected pivotal genes, henceforth referred to as “radiogenes,” which participate in radioresistance-relevant pathways and mechanisms.

Figure 4

Donut chart depicting the over-represented cancer-relevant pathways across the 88 common genes.

Table 1

Gene symbol and expression status of the 36 differentially expressed radiogenes (radioresistant vs. radiosensitive cancer cells).

Gene	Expression status	Cell line
AKT1	Up	BC
ATM	Up	BC
BAX	Down	BC; ALL
BBC3	Down	BC
BCL2	Up	BC
BIRC5	Up	BC
BRCA1	Up	BC
BRCA2	Up	BC
CASP3	Down	BC
CCND1	Up	BC
CHEK1	UP	BC
EGLN1	Up	BC
HIF1A	Up	BC; ALL
JAK1	Down	BC
JUN	Up	BC
MAP2K1	UP	BC
MAP2K2	Up	BC
MCL1	Up	BC
MSH2	Up	BC; COAD
NBN	Up	BC
NFKB1	Up	BC
NFKBIA	Down	BC
PARP1	Up	BC
PLK1	Up	BC
PMAIP1	Down	BC
PRKDC	Up	BC
PTEN	Down	BC
RELA	Up	BC
RNF8	Up	BC
SOD2	Up	BC; ALL
STAT1	Down	BC
STAT3	Up	BC
TERF2	Up	BC
TP53	Down	BC
UBE2D3	Down	ALL
XIAP	Up	BC

ALL, acute lymphoblastic leukemia; BC, breast cancer; COAD, colorectal adenocarcinoma.

Those radiogenes found to be differentially expressed between TCGA-derived cancer tissue and corresponding GTEx matched normal tissue are italicized.

Figure 5

STRING interaction network of the products of the 36 radiogenes. The nodes represent proteins and the edges indicate different modes of interactions with a confidence score ≥ 0.9.

Differentially Expressed Radiogenes in Cancer vs. Normal Tissue

Collectively, 19 radiogenes, found to be differentially expressed in the three radioresistant vs. the radiosensitive cancer cell lines (Table 1), are consistently deregulated in their corresponding cancer-matched normal tissues with the same direction (i.e., up- or down-regulated) (Supplementary Figure 2). This finding could be utilized in personalized tumor treatment for the selective eradication of cancer cells by applying radiotherapy without harming the adjacent healthy tissue at the same time.

Radiogenes Are Potential Cancer Prognostic Markers

The differential expression of the radiogenes can also predict the clinical outcome of cancer patients. In particular, a statistically significant association was found between CHEK1, MAP2K1, and PLK1 overexpression and worse overall survival in breast invasive carcinoma patients, indicated by pooled hazard ratio (HR) values higher than 1 and p-values lower than 0.05. Conversely, a significant relationship was observed between high expression of NFKBIA, which is otherwise underexpressed in radioresistant breast cancer cells, and favorable prognosis, indicated by an HR value <1 (Supplementary Figure 3).

Discussion

Cancer cells confer resistance to irradiation through diverse mechanisms including enhanced DNA damage repair capacity, activation of cell survival signaling pathways, inhibition of apoptotic pathways, and induced autophagy.

DDR/R is a complex entangled process consisting of recognition (or detection, or sensing), signaling, and repair of DNA damage (Rouse and Jackson, 2002). Ionizing radiation usually generates a variety of DNA lesions, including abasic (apurinic/apyrimidinic) sites, oxidized bases, crosslinks between thymine and cytosine bases, DNA single-strand breaks (SSBs), and DNA double-strand breaks (DSBs) (Sutherland et al., 2000). In the case of SSBs, only one strand of the double stranded DNA helix is severed. SSBs are recognized and processed by base excision repair (BER) and nucleotide excision repair (NER) mechanisms (Caldecott, 2008). DSB is the most detrimental type of DNA lesions since both strands of the double helix are broken. DSBs are repaired mainly through homologous recombination (HR) if cells are present in the S/G2 cell-cycle phase or, the less accurate, non-homologous end joining (NHEJ) (Budman and Chu, 2005). Those different types of DNA lesions can be formed either separately or in close vicinity (a few nm), resulting in clustered DNA lesions or locally multiply damaged sites (MDS). Clustered DNA damage, the hallmark of IR, is considered the most severe type of genomic damage because of its complexity. This complex DNA damage includes both DSBs and non-DSBs, usually referred to as oxidatively-clustered DNA lesions (OCDLs), occurring within a DNA region of 15–20 bp. The corresponding DNA damage repair mechanisms are recruited by the cell in response to clustered damaged sites. However, harnessing the corresponding DNA repair machinery to the clustered damaged sites is quite challenging since the presence of a lesion in one strand can delay significantly the simultaneous processing of another lesion on the complementary strand. Furthermore, OCDLs can be rapidly converted into de novo DSBs by a DNA glycosylase during the repair. As a result, unrepaired MDS can lead to increasing levels of genotoxic damage, triggering also systemic responses (Nikitaki et al., 2016). Clustered DNA lesions, if not properly processed, could contribute to increased genomic instability in the form of chromosomal aberrations (e.g., deletions and inversions) and microsatellite instability, leading eventually to carcinogenesis. Therefore, the induction of clustered DNA damage increases the cytotoxic effect of radiation, especially in highly proliferating cancer cells. Radioresistant cancer cells though can counteract this effect through their ability to respond to and repair complex/clustered DNA damage more efficiently without avoiding necessarily increased genomic instability, as compared to their radiosensitive counterparts (Hada and Georgakilas, 2008; Georgakilas et al., 2013; Mavragani et al., 2017, 2019; Bukowska and Karwowski, 2018). The radiogenes identified in the present study are implicated in diverse mechanisms underlying the acquisition of radioresistance in cancer cells.

Among the 36 “radiogenes”, ATM, BRCA1/2, CHEK1, CCND1, MSH2, NBN, PARP1, PLK1, PRKDC, and RNF8, which are consistently up-regulated across the radioresistant cancer cell lines (Table 1) and the corresponding cancer tissue (Supplementary Figure 2), are crucially implicated in different stages of DDR/R. In particular, ATM plays a protagonistic role in the initial stage of DDR/R, that is, DNA damage detection and stress-response signaling (Iliakis et al., 2003; Yang et al., 2004). ATM signaling is activated by a wide variety of DNA lesions, as well as DNA replication stress (Marechal and Zou, 2013; Burger et al., 2019). Mutations in ATM result in the genetic disorder ataxia-telangiectasia (AT), which is characterized by the high sensitivity of AT patients to IR and cancer predisposition (Gatti et al., 1988). Checkpoint kinase 1 (CHEK1), which acts downstream of ATM, is a core regulator of cell cycle checkpoint signaling in DNA damage response (Flaggs et al., 1997; Patil et al., 2013). Cyclin D1 (CCND1) was demonstrated to induce post-DNA damage cell cycle arrest and apoptosis in different types of cancers (Cai et al., 2012; Smith et al., 2016). PRKDC is a serine/threonine DNA-activated protein kinase involved in DSB recognition and DNA damage repair through NHEJ (Soubeyrand et al., 2003). Notably, ATM and PRKDC were found to affect greatly cancer cell response to IR through genome-wide genetic screening in a recent study by Francica et al. (2020). Moreover, BRCA1 and BRCA2 are largely involved in cell cycle control and maintenance of genomic stability in response to DNA damage (Deng, 2006; Gudmundsdottir and Ashworth, 2006). Another radiogene, NBN (nibrin), encodes one of the three components of the MRN complex (MRE11A-RAD50-NBN), which is implicated in the recognition and repair of DSBs (Lamarche et al., 2010). NBN is mutated in patients with Nijmegen breakage syndrome (NBS), and cells from NBS patients are hypersensitive to IR (Taalman et al., 1983). In addition, a functional link between ATM and NBN proteins has been demonstrated by Zhao et al. (2000). Also, the MMR protein MSH2 (MutS homolog 2) is suggested to contribute to radioresistance via SSB processing (Li et al., 2016). Moreover, RNF8 (ring finger protein 8) protein catalyzes the mono-ubiquitination of histones H2A and H2B during DNA damage, thereby facilitating DNA damage repair and activation of cell cycle checkpoint (Kolas et al., 2007; Ma et al., 2011); RNF8 is associated with radioresistance in human nasopharyngeal cancer cells (Wang et al., 2015). The protein encoded by the radiogene PLK1 (Polo-like kinase 1) is involved in cell cycle resumption following DNA damage-induced checkpoint arrest (Hyun et al., 2014; Bruinsma et al., 2017). It has been demonstrated that inhibition of PLK1 renders glioblastoma and non-small cell lung cancer cells sensitive to IR (Pezuk et al., 2013; Van den Bossche et al., 2019). Of particular note, pharmacological inhibitors of the BER enzyme PARP1 [poly(ADP-ribose) polymerase 1], such as niraparib, olaparib, and rucaparib, are widely used for targeted cancer therapy (Sulai and Tan, 2018; Patel et al., 2020). Importantly, CHEK1 and PLK1 were found to be poor prognostic biomarkers for the survival of breast cancer patients (Supplementary Figure 3), further supporting that enhanced DNA damage repair mechanisms in cancer cells play a catastaltic role in efficient radiotherapy (Pavlopoulou et al., 2016, 2017). Therefore, DDR/R might represent a primary “danger” signal, leading eventually to the activation of downstream signaling cascades and pro-survival mechanisms (Nikitaki et al., 2015).

The apoptotic pathway is also affected by the cellular response to radiation-induced genomic damage in cancers, as we suggest in this study. Radioresistant cancer cells have developed the ability to evade apoptosis prompted by their response to extreme and repair-resistant DNA damage, mainly due to deregulation of key pro-apoptotic molecules like TP53 (Fridman and Lowe, 2003; Haupt et al., 2003), PTEN (Lu et al., 2016), PMAIP1 (mediator of damage-induced p53-mediated apoptosis) (Oda et al., 2000), BBC3 (TP53-upregulated modulator of apoptosis) (Han et al., 2001; Nakano and Vousden, 2001), BAX (BCL2-associated X protein) (Pawlowski and Kraft, 2000), as well as anti-apoptotic proteins like MCL1 (Fujise et al., 2000; Glaser et al., 2012), BCL2 (Akl et al., 2014), BIRC5 (Chiou et al., 2003), and XIAP (X-linked inhibitor of apoptosis) (Duckett et al., 1996a; Deveraux and Reed, 1999). In a consistent manner, in this study, the pro-apoptotic genes were shown to be up-regulated, whereas the anti-apoptotic genes are down-regulated in radioresistant cancer cells (Table 1). Of those, BCL2, which is down-regulated in radioresistant breast cancer cells and tissue (Table 1, Supplementary Figure 2), can suppress apoptosis by inhibiting the activity of caspases indispensable for apoptosis, such as caspase-3 (CASP3) (Porter and Janicke, 1999; Swanton et al., 1999). Of note, aberrant activation of the Notch signaling pathway was demonstrated to inhibit TP53-mediated damage response and promote breast carcinogenesis by preventing apoptosis (Stylianou et al., 2006), suggesting a link between pro-survival and apoptotic mechanisms. PTEN was also found to promote apoptosis and cell cycle arrest via PI3K/AKT-dependent and -independent signaling (Weng et al., 2001).

The radiogene NFKB1, which transactivates several pro-inflammatory genes (Liu et al., 2017), was found to be overexpressed in radioresistant breast cancer cells and tissue (Table 1, Supplementary Figure 2). Furthermore, the members of the NFKB1 family, NFKBIA and RELA, can act either as inducers or inhibitors of apoptosis, respectively (Sonenshein, 1997; Barkett and Gilmore, 1999), consistent with their expression status in RR cells (Table 1; over- and under-expressed). The well-known NFKB1 inhibitor alpha (NFKBIA) was shown to be a favorable predictor in breast cancer patients (Supplementary Figure 3). Moreover, NFKB1, together with the inflammatory factors HIF1A (hypoxia-inducible factor 1) and STAT3, both of which were found to be up-regulated in radioresistant breast cancer cells (Table 1) and HIF1A overexpressed in breast cancer tissue (Supplementary Figure 2), are critically implicated in cancer radioresistance and radiation-induced inflammatory responses (Multhoff and Radons, 2012). On the other hand, STAT1, found down-regulated in the same cell lines (Table 1), as opposed to STAT3, elicits pro-apoptotic and anti-proliferative responses and promotes anti-cancer immunity (Avalle et al., 2012).

Suppression of UBE2D3, which is down-regulated in radioresistant ALL cells (Table 1), was demonstrated in a study by Wang et al. (2013) to decrease radiosensitivity of human breast cancer cells by promoting telomere maintenance. In addition, UBE2D3 is negatively correlated with TERF2 (telomeric repeat binding factor 2), the latter of which is primarily involved in telomere maintenance (Kim et al., 2009) and is down-regulated in radioresistance (Table 1).

Radiation can also exert its genotoxic and cytotoxic effects through the indirect and systemic induction of severe oxidative stress and the production of ROS in the organism (Kryston et al., 2011). The radiogene HIF1A plays a pivotal cytoprotective role against oxidative stress (Li et al., 2019) by inhibiting autophagy and cell death (Pouyssegur et al., 2006). Moreover, the superoxide scavenging enzyme SOD2 (superoxide dismutase 2), found overexpressed both in radioresistant breast cancer and ALL cells (Table 1), at normal expression levels provides a cytoprotective effect. Thus, SOD2 can likely exert a protective effect on RR cancer cells by controlling potential ROS-mediated DNA damage via catalyzing the reduction of superoxide into less genotoxic molecules like oxygen (Wang et al., 2018).

Notably, in this study, pro-survival pathways (Figure 3) like Notch signaling, were found to be implicated both in solid and blood cancers, probably by mediating survival of cancer cells to radiation-induced clustered/complex DNA damage (Marston et al., 2009; Capaccione and Pine, 2013). Moreover, the PI3K/AKT/mTOR signaling pathway is suggested to be important in regulating cell survival in response to different types of cellular stress (Hung et al., 2012; Porta et al., 2014), including genotoxic stress. Hence, pro-survival pathways could be considered as potential therapeutic targets in cancer (Porta et al., 2014).

A major limitation of this study, particularly for the two solid cancers datasets, is the lack of patient-derived tumor tissue, as well as cancer stem cells, which constitute a subpopulation in solid tumors that display stem-like characteristics (Pavlopoulou et al., 2016). Instead, the respective experimental studies relied on the use of commercial cancer cell lines; of particular note, the ALL cells were derived from “real” patients. However, the extent to which the individual cancer cell lines can capture the cellular and genomic complexity of tumors is questioned. Further research is needed to determine whether the results derived from the cancer cell lines investigated in this study or other studies could be extrapolated to their corresponding tissues of origin, like breast tumor tissue and colorectal adenocarcinoma. Nevertheless, it is suggested that large panels of cancer cell lines can faithfully capture the genomic heterogeneity of cancers (Sakellaropoulos et al., 2019; Vougas et al., 2019). Beyond the discussed limitations, the over- or under-expressions of the radiogenes in radioresistant phenotypes have been verified, to a great extent, by independent experimental biochemical studies available in the literature. In addition, the clinical implications of those genes are further supported by the patient survival results (Supplementary Figure 3).

Conclusion

In the present study, we employed an integrative bioinformatics approach to analyze transcriptomic data regarding the molecular determinants of cancer cell radioresistance. On the basis of our findings, both solid and hematologic cancer cells likely depend on similar mechanisms to confer resistance to IR (i.e., DDR/R and cell survival). Moreover, we identified 36 functionally associated radiogenes that participate in radioresistance-associated pathways. Most of those radiogenes were also shown to be differentially regulated in the corresponding cancer tissues. Moreover, several of the radiogenes were found to have potential prognostic value for the clinical outcome of cancer patients. However, the availability of clinically derived cancer tissues would provide a more reliable source for conducting research on the response of cancer patients to radiation. The overall data presented herein can be particularly useful for clinicians in selecting suitable targets (e.g., DDR/R inhibitors) for appropriate combination therapy using IR. In conclusion, we suggest that this bioinformatics premise can be harnessed as a first step in the rational design of in vivo experimental studies or in personalized medicine for optimizing tumor response and cancer cell susceptibility to therapeutic ionizing radiation and reduction of the total effective radiation dose administered to the patient.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author/s.

Author Contributions

AP conceived the study. AP and AGG designed and supervised the study. HIT, GK, PK, and AP performed the experiments and analyzed the data. HIT, GK, PK, HA, AGG, and AP wrote the manuscript. All authors reviewed and approved of the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.