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.

   Keywords: ionizing radiation damage, radiation-induced bystander
   effect, transcription factor, WGCNA

Introduction

   In the past 20 years, radiotherapy has become the primary cytotoxic
   therapy for localized solid cancer.^[43]1 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.^[44]2
   However, IR has been shown to cause severe cell damage.^[45]3 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.^[46]4 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.^[47]3,[48]5 Endocrine molecules of irradiated
   tumor cells may cause damage to adjacent normal cells.^[49]6 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).^[50]7 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.^[51]8-[52]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 ([53]GSE65292, [54]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
   [55]GSE65292^[56]12 were downloaded from the GEO database
   ([57]https://www.ncbi.nlm.nih.gov/geo/),^[58]13 which was based on
   [59]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, [60]GSE55953^[61]14,[62]15 based on
   [63]GPL14550, including 8 IR samples and 20 normal controls, and
   [64]GSE21240 based on [65]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^[66]16 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 [67]Figure 1.

Figure 1.

   [68]Figure 1.
   [69]Open in a new tab

   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^[70]17 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^[71]18 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 [72]GSE65292

   All DEGs of 5 comparisons in [73]GSE65292 were extracted to perform
   WGCNA.^[74]19 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^[75]20 in R. A P < 0.05 was considered
   significant. In addition, ClueGO^[76]21 in Cytoscape^[77]22 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.^[78]23
   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^[79]24 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 ([80]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 ([81]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 ([82]Figure 2C).

Figure 2.

   [83]Figure 2.
   [84]Open in a new tab

   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 ([85]Figure
   3A), and a total of 4 functional modules were identified ([86]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 ([87]Figure 3C). We also constructed a topological overlap matrix
   plot and found that the genes in all modules have strong coexpression
   relationships ([88]Figure 3D).

Figure 3.

   [89]Figure 3.
   [90]Open in a new tab

   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 ([91]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 ([92]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 ([93]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 ([94]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 ([95]Figure 4E).

Figure 4.

   [96]Figure 4.
   [97]Open in a new tab

   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; [98]Figure 5A). Therefore, a TF-target genes-pathway
   global regulatory network was constructed ([99]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
   ([100]Figure 5C).

Figure 5.

   [101]Figure 5.
   [102]Open in a new tab

   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 ([103]Figure 6A), which were further
   validated in [104]GSE55953 ([105]Figure 6B) and [106]GSE21240
   ([107]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 ([108]Figure 7A), and
   these aberrant expressions were validated by [109]GSE55953 ([110]Figure
   7B) and [111]GSE21240 ([112]Figure 7C).

Figure 6.

   [113]Figure 6.

   [114]Figure 6.
   [115]Open in a new tab

   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 [116]GSE55953. C, ROC curve
   analysis of core molecules in [117]GSE21240.

Figure 7.

   [118]Figure 7.
   [119]Open in a new tab

   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 [120]GSE55953. C, Core
   molecules were upregulated in IR damage of [121]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.^[122]4,[123]25,[124]26 Existing literatures indicated
   that RIBE may be mediated directly by gap junction intercellular
   communication and/or divergent cellular factors increased from infrared
   cells.^[125]27-[126]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.^[127]30 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,^[128]31,[129]32 thus elicited a cellular stress response that
   includes DNA damage recognition and cell-cycle arrest, followed by DNA
   repair or apoptosis.^[130]33 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,^[131]34 and the kinase activity of cells of ATM
   can be activated by IR irradiation.^[132]35,[133]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.^[134]37 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).^[135]38 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.^[136]39-[137]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.^[138]38,[139]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.^[140]44 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
   [141]Click here for additional data file.^ (75.8KB, pdf)

   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

Acknowledgments