Blood Gene Expression Profile Study Revealed the Activation of Apoptosis and p53 Signaling Pathway May Be the Potential Molecular Mechanisms of Ionizing Radiation Damage and Radiation-Induced Bystander Effects

Abstract

Radiotherapy is an effective treatment for local solid tumors, but the mechanism of damage to human body caused by radiation therapy needs further study. In this study, gene expression profiles of human peripheral blood samples exposed to different doses and rates of ionizing radiation (IR) were used for bioinformatics analysis to investigate the mechanism of IR damage and radiation-induced bystander effect (RIBE). Differentially expressed genes analysis, weighted gene correlation network analysis, functional enrichment analysis, hypergeometric test, gene set enrichment analysis, and gene set variation analysis were applied to analyze the data. Moreover, receiver operating characteristic curve analysis was performed to identify core genes of IR damage. Weighted gene correlation network analysis identified 3 modules associated with IR damage, 2 were positively correlated and 1 was negatively correlated. The analysis showed that the positively correlated modules were significantly involved in apoptosis and p53 signaling pathway, and ESR1, ATM, and MYC were potential transcription factors regulating these modules. Thus, the study suggested that apoptosis and p53 signaling pathway may be the potential molecular mechanisms of IR damage and RIBE, which could be driven by ESR1, ATM, and MYC.

Introduction

In the past 20 years, radiotherapy has become the primary cytotoxic therapy for localized solid cancer. ¹ Radiotherapy using ionizing radiation (IR) to destroy cancer cells by irradiating cancerous tumors with high doses of radiation produced by special equipment, thus inhibiting their growth, reproduction, and proliferation. ² However, IR has been shown to cause severe cell damage. ³ Cells exposed to IR showed increased frequency of DNA damage, apoptosis, and chromosomal aberration or mutation. For a long time, it had been generally believed that these biological events are mainly caused by the action of reactive oxygen species formed by ionization and water irradiation of cell structure. Recently, however, attentions have been paid to the bystander effect. ⁴ Radiation-induced bystander effect (RIBE) is induced by reagents and signals emitted by directly irradiated cells, which cause the lowering of survival, cytogenetic damage, apoptosis enhancement, and biochemical changes in neighboring nonirradiated cells. ^3,5 Endocrine molecules of irradiated tumor cells may cause damage to adjacent normal cells. ⁶ Cells exposed to IR and other genotoxic agents (targeted cells) can communicate their DNA damage response (DDR) status to cells that have not been directly irradiated (bystander cells). ⁷ Molecular signals can be transmitted through gap junction, intercellular communication, and mediator transfer mechanisms. However, the biological mechanism of this effect is still not well understood. Many studies have proposed to estimate the dose in irradiated peripheral blood based on gene expressions and showed that peripheral blood is capable of reflecting the degree of radiation exposure. ^{8 -11} In addition, the irradiated tumor tissue contains many capillaries, and factors or biomarkers produced by irradiated tissue can enter the bloodstream through capillaries, though the radiation dosage is relatively small, it could induce RIBE. Therefore, we proposed to use gene expression profiles in peripheral blood to study bystander effects and IR loss, since the gene expression patterns of irradiated peripheral blood may provide information on the molecular mechanism of RIBE.

In our study, radiation-induced gene expression of peripheral blood at the 24-hour time point after irradiation was used to perform weighted correlation network analysis (WGCNA). The analysis showed that apoptosis and p53 signaling pathway (as the major Kyoto Encyclopedia of Genes and Genomes [KEGG] pathway of functional module) may be the potential mechanism of bystander effects and IR. In addition, strong changes were identified in gene expressions in both irradiated and nonirradiated cells by hypergeometric testing. We also found that 3 transcription factors (TFs; ESR1, ATM, and MYC) and their downstream target genes may play a central role in the biological response to IR. Subsequently, a TF-target genes pathway global regulatory network was constructed for the involved genes. Furthermore, our study also identified 13 core molecules that have the ability to distinguish radiation exposure from nonradiation exposure in peripheral blood.

Materials and Methods

Data Selection

Three data sets (GSE65292, GSE55953, and GSE212400) were found in the Gene Expression Omnibus (GEO) database by searching with the keywords “radiation” and “blood gene expression profiles.”

Data Processing

The human peripheral blood gene expression profiles of GSE65292 ¹² were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), ¹³ which was based on GPL13497. The profiles contain 30 irradiated human whole blood samples and 5 normal controls. The doses (0.56, 2.2, and 4.45 Gy) were delivered by 2 dose rates, acute dose rate of 1.1 Gy/min and low-dose rate of 3.1 mGy/min. In addition, GSE55953 ^14,15 based on GPL14550, including 8 IR samples and 20 normal controls, and GSE21240 based on GPL6480, including 24 IR samples and 24 controls, were used to validate the aberrant expressions of genes of interest. The normalizeBetweenArrays function in the limma package ¹⁶ was used to normalize the gene expression profiles. If a gene responded to multiple probes, the average value of these probes was considered to be the expression value of the corresponding gene. The workflow of the study is shown in Figure 1.

Figure 1.

Flowchart of the present study.

Gene Set Enrichment Analysis and Gene Set Variation Analysis

Gene set enrichment analysis (GSEA) was performed using the normalized gene expression profiles to explore the biological process (BP) and KEGG pathways in relation with different dose- and rate-radiation damage. The Java software of GSEA (version 2-2.2.4) was used in the analysis. The c5.bp.v6.2.symbols.gmt and c2.cp.kegg.v6.2.symbols.gmt data sets in MsigDB V6.2 database ¹⁷ were used as reference gene sets, and GSEA was performed according to default parameters. P < 0.05 was considered significant. In addition, gene set variation analysis (GSVA) package ¹⁸ in R was used to estimate the expression of the gene set in the individual samples.

Differentially Expressed Gene Analysis

Compared to the control samples, the differentially expressed genes (DEGs) in 0.56 Gy dose samples, 2.2 Gy dose samples, 4.45 Gy dose samples, 1.1 Gy/min rate samples, and 3.1-mGy/min rate samples were analyzed using the limma package in R. The genes with P adjusted by the false discovery rate <.01 were considered significant.

Weighted Gene Correlation Network Analysis in GSE65292

All DEGs of 5 comparisons in GSE65292 were extracted to perform WGCNA. ¹⁹ First, hclust function was used for hierarchical clustering analysis. Then, the soft thresholding power value was screened during module construction by the pickSoftThreshold function. Candidate power (1-30) was used to test the average connectivity degrees of different modules and their independence. In the analysis, the power values were automatically estimated by WGCNA. The WGCNA R package was also used to construct coexpression networks (modules), where the minimum module size was set to 30 and each module was assigned a unique color label.

Functional Enrichment Analysis

To further explore the biological significance of the functional modules, Gene Ontology (GO) and KEGG pathway enrichment analyses for the module genes were performed, respectively, using the clusterProfiler package ²⁰ in R. A P < 0.05 was considered significant. In addition, ClueGO ²¹ in Cytoscape ²² was used to perform BPs enrichment analysis for each module.

Hypergeometric Test and Correlation Analysis

In order to predict the upstream TFs of the regulatory modules, hypergeometric test was carried out. Interactions between TFs and their target genes were downloaded from TRRUST v2 database. ²³ Interactions between a regulator and a related functional module were examined using the hypergeometric test in R. Interactions between a regulator and a functional module that showed quantity >2 and P <0.05 were considered significant. Moreover, we analyzed Pearson correlation between a TF and each of its targets to reduce noise and false positives. A P < 0.05 was considered significant. Combining enrichment analysis results, a TF-target genes-pathway network was constructed.

Hub Genes and Receiver Operating Characteristic Curve

In WGCNA, gene significance (GS) is defined as the correlation a gene with phenotype. The module membership (MM) is defined to measure the importance of a gene in the module. In this study, a gene with GS > 0.8 and MM > 0.9 was defined as a hub gene among the candidate modules. Then, the pROC package ²⁴ was used to conduct the receiver operating characteristic (ROC) curve analysis for these hub genes and TF and their target genes. Molecules with area under ROC (AUC) >0.9 were identified as the core genes of IR damage.

Result

Molecular Imbalance in Peripheral Blood Induced by Gradient Dose and Rate of Radiotherapy

Compared to control samples, there are 59 DEGs in 0.56 Gy dose, 3 of which were upregulated, 56 of which were downregulated; 167 DEGs in 2.2 Gy dose, 52 of which were upregulated, 115 of which were downregulated; 281 DEGs in 4.45 Gy dose, 108 of which were upregulated, 173 of which were downregulated; 59 DEGs in 1.1 Gy/min, 5 of which were upregulated, 54 of which were downregulated; 291 DEGs in 3.1 mGy/min, 114 of which were upregulated, 177 of which were downregulated. The 3 genes with the highest significance (ranked by fold-change) in each comparison were visualized (Figure 2A), and these genes may be affected by radiation the most. We found 54 common DEGs in peripheral blood samples that received different radiation doses and rates (Figure 2B). They may be radiation-specific genes. In addition, the expressions of the 54 common DEGs can basically distinguish between different cases and controls (Figure 2C).

Figure 2.

Differentially expressed gene (DEG) analysis. A, Manhattan plot for associations between different degree ionizing radiation damages and control, the first 3 genes with the highest significance were highlighted. B, Common DEGs in dose_0.56 Gy control, dose_2.2 Gy control, dose_4.45 Gy control, rate_1.1 Gy/min control, and rate_3.1 mGy/min control. C, (a) Heatmap of 54 common DEGs in different dose controls. (b) Heatmap of 54 common DEGs in different rate controls.

Synergistic Expressions of Disordered Molecules Were Observed and Significantly Correlated With Doses and Rates

In the WGCNA, power = 0.86 was estimated by the package (Figure 3A), and a total of 4 functional modules were identified (Figure 3B). Each module was given an individual color as its identifier. The colors are brown, blue, turquoise, and gray. The analysis showed that these functional modules were significantly correlated with the dose and the rate of acceptance of IR, of which the brown module is negatively correlated with the dose and rate of IR. The blue and the turguoise modules are positively correlated with the dose and rate of IR (Figure 3C). We also constructed a topological overlap matrix plot and found that the genes in all modules have strong coexpression relationships (Figure 3D).

Figure 3.

Weighted gene correlation analysis. A, Dynamic branch cutting. B, Hierarchical clustering, each color represents a functional module. C, Correlation heatmap of gene modules and phenotypes. Red indicates a positive correlation, and blue indicates a negative correlation. D, Topological overlap matrix (TOM) plot for the gene coexpression network of the intramodules. In the TOM plot, the light color indicates low topological overlap, while the darker color indicates high topological overlap.

Biological Processes and Pathways Involved in Dose- and Rate-Dependent Modules

To study the underlying biological functions of these modules, functional enrichment analysis was carried out. The results of GO analysis (Figure 4A) revealed that genes in the blue module were significantly involved in entry into host cells, entry into hosts, entry into other organism involved in symbiotic interaction, and entry into cell of other organisms involved in symbiotic interaction. Genes in the brown module were significantly involved in B-cell proliferation, B-cell activation, and dendrite extension. Genes in the turquoise module were significantly involved in signal transduction in response to DNA damage, DDR, and signal transduction by p53 class mediator. An analysis of the KEGG pathway (Figure 4B) showed that the blue module was significantly involved in natural killer cell–mediated cytotoxicity and FoxO signaling pathway, the brown module was significantly involved in PI3K-Akt signaling pathway, and the turquoise module was significantly involved in p53 signaling pathway, necroptosis, human papillomavirus infection, apoptosis, and alcoholism. In addition, Clue Go analysis indicates (Figure 4C) that the brown module was involved in deoxyribonucleotide metabolic process, muscle cell cellular homeostasis, and dendrite extension; the blue module was involved in negative regulation of cell-matrix adhesion; and the turquoise module was involved in signal transduction by p53 class mediator, response to UV cellular, response to UV, and DNA cytosine deamination. Moreover, the GSEA results indicated that apoptosis pathway was significantly enriched in the samples of 0.56 Gy dose, 2.2 Gy dose, and 1.1 Gy/min rate, and p53 signaling pathway was significantly enriched in samples of 0.56 Gy dose, 2.2 Gy dose, 4.45 Gy dose, 1.1 Gy/min rate; and 3.1 mGy/min rate (Figure 4D). The GSVA results indicated the gene set variations in both apoptosis and p53 signaling pathway were gradually increased with the increase in dose and rate (Figure 4E).

Figure 4.

Biological processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analysis of functional modules. A, Biological processes of functional modules, similar biological processes were clustered. B, KEGG pathways of functional modules. C, Biological processes of functional modules in Clue Gene Ontology. D, Apoptosis and P53 signaling pathway were enriched in different samples. In apoptosis, the lines of 0.56 Gy and 1.1 Gy/min were overlapping. E, The gene set variations of both apoptosis and p53 signaling pathway were gradually increased with the increase of dose and rate.

Potential Mechanism of IR Damage

Hypergeometric test result showed that in the turquoise module, 3 TFs (ESR1, ATM, and MYC; Table S1) were significantly correlated with 6 target genes (TNFRSF10B, MDM2, CDKN1A, FAS [also known as TNFRSF6], PCNA, and GADD45A; Figure 5A). Therefore, a TF-target genes-pathway global regulatory network was constructed (Figure 5B). According to gene differential expression analysis, most DEGs regulated by TF in apoptosis and p53 signaling pathway show high expressions, and the expressions of DEGs increased with the increase in radiation dose and rate. These may be relevant to the degree of IR damage, which has been confirmed in GSEA. Thus, apoptosis and p53 signaling pathway were identified as the potential mechanism of RIBE and IR damage (Figure 5C).

Figure 5.

Potential mechanism of IR damage. A, The correlation between TFs and its target genes. B, TF-target genes pathway network. C, Potential mechanism of IR damage. Red genes represent the target genes of TFs (ATM, ESR1, and MYC). IR indicates ionizing radiation; TF, transcription factor.

Core Molecules of Gradient Dose- and Rate-Radiation Damage

Furthermore, according to GS > 0.8 and MM > 0.9, 14 genes were identified as hub genes (PHLDA3, GADD45A, PCNA, TNFSF4, MDM2, E2F7, VWCE, FHL2, DRAM1, FDXR, DDB2, PVT1, ZMAT3, and PAPPA). Therefore, ROC curve analyses were performed on hub genes, TF of participation mechanism, and TF target genes. The result suggested that 13 genes have a strong ability to distinguish whether or not it was exposed to IR and were identified as core genes (Figure 6A), which were further validated in GSE55953 (Figure 6B) and GSE21240 (Figure 6C). In addition, it was found by comparing with the controls, the expressions of these genes were gradually upregulated with the increase in rate and dose in IR sample (Figure 7A), and these aberrant expressions were validated by GSE55953 (Figure 7B) and GSE21240 (Figure 7C).

Figure 6.

Identification of core molecules in ionizing radiation damage. A, (a) Receiver operating characteristic (ROC) curve analysis of core molecules by doses. (b) ROC curve analysis of core molecules by rates. B, ROC curve analysis of core molecules in GSE55953. C, ROC curve analysis of core molecules in GSE21240.

Figure 7.

Trends of gene expressions in different samples. A, Core molecules were gradually upregulated in different doses ionizing radiation (IR) damages and were gradually upregulated in different rate IR damages. B, Core molecules were upregulated in IR damage of GSE55953. C, Core molecules were upregulated in IR damage of GSE21240.

Discussion

Although radiotherapy is a local treatment, the body’s response to radiotherapy is systematic. There are no circulating biomarkers for measuring this systemic response. In our current study, we sought to explore whether this systemic response is reflected in the peripheral blood transcriptome. Radiation-induced bystander effect is one of the systemic responses of the body to radiotherapy. When subjected to radiation therapy, adjacent tissues are inevitably affected by radiation, which may produce RIBE. RIBE refers to the process in which factors released by irradiated cells or tissues affect other parts of nonirradiated parts of the body animals, leading to genomic instability, stress response, and changes in apoptosis or cell proliferation. ^4,25,26 Existing literatures indicated that RIBE may be mediated directly by gap junction intercellular communication and/or divergent cellular factors increased from infrared cells. ^{27 -29} However, it has also been indicated in literatures that gap junction inhibitors or enhancers had no effect on the bystander effect of human lung cancer cell lines or rat tumor cell lines. ³⁰ Therefore, more study on pptential of mechanism of RIBE are necessary. To further explore these important issues in radiotherapy, we carried out a comprehensive analysis of peripheral blood samples irradiated at different doses and rates in this study. We found that the DEGs of samples irradiated by different doses and rates are different. Our WGCNA result identified 4 functional modules, which have similar statistical significances on dose and rate dependence. These modules have different correlations with radiation dose and rate, indicating they may perform different functions in the body's response to radiation. Among them, the coexpression of gray module genes was not significant, and the coexpression of the other 3 modules were significant. Among these 3, the blue module genes were mainly involved in functions such as entry into hosts, the brown module genes were mainly involved in functions such as activation and differentiation of immune cells, and the turquoise module genes are mainly involved in functions such as DNA damage checkpoints. Moreover, our functional enrichment analysis indicated that the functional modules were significantly involved in natural killer cell–mediated cytotoxicity, PI3K-Akt signaling pathway, p53 signaling pathway, and apoptosis. The GSEA showed p53 signaling pathway and apoptosis were also enriched in most radiated samples, suggesting that apoptosis and p53 signaling pathway may be the major potential mechanisms of IR damage. The major consequence of IR exposure is the generation of single- or double-stranded breaks in DNA, which result in DNA damage, ^31,32 thus elicited a cellular stress response that includes DNA damage recognition and cell-cycle arrest, followed by DNA repair or apoptosis. ³³ In addition, the radiated blood in the capillaries of the radiotherapy site circulates throughout the body, which may contribute to the formation of RIBE. These are supportive of our finding that apoptosis and p53 signaling pathway may be the major factor of RIBE.

In addition, 3 TFs (ESR1, ATM, and MYC) and their target genes (TNFRSF10B, MDM2, CDKN1A, FAS, PCNA, and GADD45A) were identified, and a TF-target genes pathway global regulatory network was constructed based on the hypergeometric test. Existing studies indicated that MYC was required for activation of the ATM-dependent checkpoints in response to DNA damage, ³⁴ and the kinase activity of cells of ATM can be activated by IR irradiation. ^35,36 Moreover, it was also confirmed that the N-terminal of ATM and ATM’s kinase activity were required for activation of p53’s transcriptional activity and restoration of normal sensitivity to DNA damage. ³⁷ Furthermore, 6 TF target genes are involved in a variety of pathways, including DNA repair (PCNA), cell-cycle progression (CDKN1A and GADD45A), and cell death (TNFRSF10B and TNFRS6). ³⁸ It was also reported in the literature that DNA damage induced by IR involved P53 protein disorders. The increase in p53 protein levels thus can lead to the induction of many genes, including ACTA2, CDKN1A, DDB2, FDXR, GADD45A, PIG3, TNFRSF6, and TNFSF10B. ^{39 -42} The upregulation of these genes leads to a diverse set of events, including cell-cycle arrest, DNA repair, and apoptosis. These results of the existing studies are supportive to our founding.

Moreover, we identified 13 core molecules of IR damage: CDKN1A, DDB2, E2F7, FDXR, GADD45A, MDM2, PAPPA, PCNA, PHLDA3, PVT1, TNFRSF10B, TNFSF4, and ZMAT3. They were also upregulated in IR samples compared to controls, and most of them were involved in the p53 signaling pathway. ^38,43 Among them, CDKN1A, DDB2, FDXR, GADD45A, PCNA, and TNFRSF10B are known radiation-response genes, which have a sustained radiation dose–response and a strong positive correlation with dose. ⁴⁴ Therefore, these genes may be part of the main circulating biomarkers of IR damage. If this is confirmed, it will be able to help us identify and understand IR-related risks associated with radiological examinations and radiotherapy. We plan to further verify our findings using experiments.

Conclusion

This study showed apoptosis and ATM-p53 signaling pathway, and ESR1, ATM, and MYC may be mainly responsible for radiotherapy damage to human body. Thus, understanding their roles in the related transcription regulatory network may help us understand the biological mechanism of radiotherapy damage.

Supplemental Material

Supplemental Material, Table_S1 - Blood Gene Expression Profile Study Revealed the Activation of Apoptosis and p53 Signaling Pathway May Be the Potential Molecular Mechanisms of Ionizing Radiation Damage and Radiation-Induced Bystander Effects

Supplemental Material, Table_S1 for Blood Gene Expression Profile Study Revealed the Activation of Apoptosis and p53 Signaling Pathway May Be the Potential Molecular Mechanisms of Ionizing Radiation Damage and Radiation-Induced Bystander Effects by Guangyao He, Anzhou Tang, Mao Xie, Wei Xia, Pengcheng Zhao, Jianglian Wei, Yongjing Lai, Xianglong Tang, Yi Ming Zou and Heng Liu in Dose-Response