Candidate biomarkers and persistent transcriptional responses after low and high dose ionizing radiation at high dose rate

Abstract

Purpose

Development of an integrated time and dose model to explore the dynamics of gene expression alterations and identify biomarkers for biodosimetry following low- and high-dose irradiations at high dose rate.

Material and methods

We utilized multiple transcriptome datasets (GSE8917, GSE43151, and GSE23515) from Gene Expression Omnibus (GEO) for identifying candidate biological dosimeters. A linear mixed-effects model with random intercept was used to explore the dose-time dynamics of transcriptional responses and to functionally characterize the time- and dose-dependent changes in gene expression.

Results

We identified genes that are correlated with dose and time and discovered two clusters of genes that are either positively or negatively correlated with both dose and time based on the parameters of the model. Genes in these two clusters may have persistent transcriptional alterations. Twelve potential transcriptional markers for dosimetry—ARHGEF3, BAX, BBC3, CCDC109B, DCP1B, DDB2, F11R, GADD45A, GSS, PLK3, TNFRSF10B, and XPC were identified. Of these genes, BAX, GSS, and TNFRSF10B are positively associated with both dose and time course, have a persistent transcriptional response, and might be better biological dosimeters.

Conclusions

With the proposed approach, we may identify candidate biomarkers that change monotonically in relation to dose, have a persistent transcriptional response, and are reliable over a wide dose range.

Keywords:

• Dosimetry biomarkers
• linear mixed model
• low- and high-dose irradiation
• persistent transcriptional response
• molecular pathways
• integrated time-dose dynamics

1. Introduction

Human beings are exposed to ionizing radiation (IR) from various sources including the environment, medical diagnostics and treatment, industrial applications, and potential mass-casualty nuclear disasters (UNSCEAR, 2008). There are several types of IR, including α particles, β particles, positrons, neutron particles, X-rays, and γ-rays (Hill 2020). Exposure to IR may cause DNA damage, cell death, and alterations in DNA, RNA, and protein levels, which can lead to a higher risk of cancer and other diseases in the long term (Kamiya et al. 2015; Ozasa et al. 2019; Ghosh and Ghosh 2021; Nagane et al. 2021). To understand the biological effects of IR exposure, it is critical to explore the underlying molecular mechanisms and identify IR-related biomarkers (IAEA, 2010).

With advances in biotechnology, multi-omics data—including gene expression—have been produced at great speed and low cost, allowing us to more efficiently investigate the cellular responses and biological effects of high- and low-dose IR. The major goals of transcriptome studies are to determine the cell- and condition-specific dynamics of transcriptional activity and to investigate regulating mechanisms of mRNA production (Brackmann et al. 2020; ICRP 2007). It is known that acute exposure to low- or high-dose IR leads to profound changes in expression patterns (Sokolov and Neumann 2015). In radiation protection, very low, low, moderate, and high doses of IR exposure are defined as $<10$ mGy, $10-<100$ mGy, $100-1000$ mGy, and $>1000$ mGy, respectively (Little et al., 2021). Transcriptomic responses to different dose ranges might involve different genes. Some genes may only have a short-term response, while others may exhibit changes over a long period of time (Jafer et al. 2020). Despite significant progress made in understanding radiation pathogenesis, IR-related persistent alterations at the molecular level are still not fully explored (Garikipati et al. 2021). However, investigating the persistent transcriptional responses of radiation may help us understand the biology of radiation-induced cellular alterations and can have potential implications for the health-effect studies of A-bomb survivors.

One potential application of radiation-induced expression changes is retrospective estimation of radiation dose (Ostheim et al. 2022). Biodosimetry refers to the measurement of biological markers that are correlated with the magnitude of the dose. Accurately estimating the received radiation dose of accidentally exposed individuals within a population might help to assist medical decision-making in rapid triage and assessing the risk of long-term consequences from radiation exposure (Abend et al. 2022). Even though the gold standard for IR biodosimetry remains the dicentric chromosome assay (DCA) (Sullivan et al., 2013; Sproull et al., 2017; Escalona et al., 2022), gene expression profiles from peripheral blood lymphocytes (PBL) have attracted great interest for biodosimetry development. Gene expression profiles are more time-efficient and less expensive, providing an appealing strategy for biodosimetry (Macaeva et al. 2019). Differentially expressed (DE) genes have been identified from microarrays or RNA-seq in several studies (Lacombe et al. 2018). However, gene expression-based biodosimetry is still an ongoing topic of research. One of the major challenges for expression as a biomarker of IR exposure is that the signal is highly transient and dynamic. Expression changes depend on both time and dose, and the temporal kinetics of expression differ by gene (Manning et al. 2013; Macaeva et al. 2019). To date, biomarkers have been mainly identified with statistical tests or supervised classification, ignoring the time and dose dynamics of the responses, which contributes to poor reproducibility (Serra et al. 2020).

The requirements of a good biological dosimeter are that it should be IR sensitive with changes linear in dose, exhibit a stable signal without temporal fluctuations (no or linear correlation with time), be reliable over a wide dose range, and be subject to no influence of gender, age, or health status (Hladik et al. 2022). Moreover, an IR biodosimetry study focuses on the expression of candidate biomarkers between 1– 7 days; therefore, the transcriptional response of a good IR biodosimeter must be stable beyond the first day after exposure. However, most whole genome expression datasets available in the public domain were measured within 24 hours of exposure (Lacombe et al. 2018) and are likely more dynamic, adding to the challenges of dosimetry biomarker identification. Therefore, in this pilot study, we propose a linear mixed effects model to simultaneously explore the effects of both dose and time on transcriptomic changes. The influence of individual characteristics such as sex and age is adjusted explicitly through the random effects; fixed effects are used to explore the IR-related molecular changes and thus identify candidate dosimetry biomarkers.

2. Materials and Methods

2.1. Datasets

As demonstrated in Table 1, three microarray datasets—GSE8917, GSE53151, and GSE23515—were downloaded from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database and utilized for this study.

Table 1. Microarray datasets studied.

	GSE8917	GSE43151	GSE23515
Study & Year	Paul & Amundson, 2008	Nosel et al., 2013	Paul & Amundson, 2011
Platform	Agilent Whole Human Genome Oligo Microarrays (D4112A)	Agilent Whole Human Genome Microarray v2, 4x44K	Agilent Whole Human Genome Oligo Microarrays (D4112A)
Donors	10 healthy donors (5 females/ 5 males)	5 healthy donors (5 males)	24 healthy donors (12 females/ 12 males)
Samples	Whole blood in sodium citrate	Whole blood in heparin diluted 1:10 in IMDM	Whole blood in sodium citrate
RNA Extraction Sources	Whole blood diluted in RPMI-1640/ 10% FBS	CD4+ T lymphocytes	Whole blood diluted in RPMI-1640/ 10% FBS
Dose (Gy)	0.5, 2, 5, 8	0.005, 0.01, 0.025, 0.05, 0.1, 0.5	0.1, 0.5, 2
Dose Rate (Gy/min)	0.82	0.05	0.82
Time after IR (hours)	6, 24	2.5, 5, 7.5, 10	6

IMDM: Iscove’s Modified dulbecco’s Medium; RPMI: Roswell park memorial institute; FBS: fetal bovineserum.

2.2. GSE8917 (high-dose irradiation)

This radiation gene expression dataset was originally generated from an Agilent DNA Microarray with peripheral blood samples obtained from 10 healthy donors including 5 men and 5 women (Paul and Amundson 2008). The dataset was obtained and preprocessed with the GEOquery R package. Fifty microarray samples were generated after irradiation of whole blood at 0.82 Gy/min to total doses of 0, 0.5, 2, 5, or 8 Gy. Irradiated blood was then diluted 1:1 with fetal bovine serum and incubated for 6 or 24 hours at 37 °C in a humidified incubator with 5% CO₂ before total blood RNA extraction, labeling, and hybridization. Each of the five doses was applied to blood samples from the same donor, but data at 6 and 24 hours of post-exposure were from different donors: at 6 hours, 2 men (mean age 46) and 3 women (mean age 30.3); at 24 hours, 3 men (mean age 38.3) and 2 women (mean age 34.5). This study has a high dose (HD) with a maximum of 8 Gy. We get the gene-level expression data by averaging out the expressions at the probe level and excluding genes with missing values. There is a total of 8639 genes after preprocessing. We also normalize both the dose and time into the range of 0–1 to facilitate visualization.

2.3. GSE43151 (low-dose irradiation)

This gene expression dataset was originally collected to study gene-expression changes in CD4+ T lymphocytes after exposure of whole blood (diluted 1:10 in IMDM [Iscove’s Modified Dulbecco’s Medium]) to low and very low doses (LD) of ionizing radiation (Nosel et al. 2013). CD4+ T lymphocytes were extracted by using negative cell sorting (i.e., labeling and removal of unwanted cell types). These data were generated from 5 male donors by using an Agilent-026652 Whole Human Genome Microarray 4x44K v2. Whole blood exposures were performed with a ⁶⁰Co source at a relatively low dose rate (50 mGy min⁻¹) and exposures were performed independently on blood samples of each donor. There are 103 microarray subjects with doses of 0, 5, 10, 25, 50, 100, and 500 mGy and times of 0, 2.5, 5, 7.5, and 10 hours post-irradiation before RNA extraction. Note that the design is not saturated; the reason why some samples at certain combinations of dose and time are not included is unknown. There are 13447 genes after we average out the expression at the probe level and exclude genes with missing values. We also normalize the dose and time into the range of 0 − 1 before model fitting.

2.4. GSE23515 (validation dataset)

This dataset was generated with Agilent’s one-color microarray technology. Peripheral blood cells from each of 24 different donors were exposed to four doses of ionizing radiation (0, 0.1, 0.5, or 2 Gy). The donors comprise 12 women and 12 men, in each sex half of whom were smokers and half nonsmokers. The 6x2x2x4 design is nearly saturated, except that there are only 5 results for female smokers at 2 Gy (total of 95 samples). Peripheral blood from volunteers was first drawn into 0.105 mol/l sodium citrate Vacutainer tubes and exposed at a rate of 0.82 Gy/min to 0, 0.1, 0.5, or 2 Gy γ-rays. Gene expressions were measured 6 hours after irradiation (Paul and Amundson 2011). We made a log transform and quantile normalization and excluded genes with missing values. There are 12152 genes available for downstream analysis. This dataset does not have a time dimension and is mainly used to validate the biomarkers identified from the other datasets.

2.5. Methods

2.5.1. Linear Mixed-Effects (LME) model

A linear mixed-effects (LME) model (linear mixed model) is a linear regression model containing both fixed effects and random effects. When the observations are not independent, the random effects in the model allow for an individual-level intercept and/or slope, so that the parameters of a linear regression are no longer fixed but rather can vary around their mean values. In the following model, we account for the variations in baseline gene expression for different persons by including an additional random intercept. Assuming that the baseline expression level of a gene without exposure at time 0 is $E_0{e}^{{\beta }_0}$ and the expression level at time $T$ and dose $D$ is $E,$ a simple dose-response model is defined as: $$E=\mathrm{}{E}_0{e}^{{\mathrm{\beta }}_0+{\mathrm{\beta }}_{1}D+{\mathrm{\beta }}_{2}T+{\mu }_I+\epsilon }\mathrm{}.$$

Taking the natural log, The LME model is specified as: $$\mathit{LFC}=\log \frac{E}{E_0}={\beta }_0+{\beta }_{1}D\mathrm{}+\mathrm{}{\beta }_{2}T\mathrm{}+{\mu }_{I\mathrm{}}+\epsilon ,$$ where LFC is defined as the log-fold change, ${\beta }_0,{\beta }_1,{\beta }_2$ are the fixed-effect parameters, ${\mu }_I$ measures the random effect of the intercept that varies randomly for each person, and ${\mu }_I\sim N(0,{\sigma }_{\mu }^2),$ $\epsilon \sim N(0,{\sigma }^2)\mathit{and}\mathit{cov}({\mu }_I,\epsilon )=0.$

To detect the linear association between dose/time and LFC, it is assumed that individuals only have a random effect on the intercept, i.e., that there are no random effects on dose and time course, so personal characteristics have no influence on the dose and time associations with LFC. Although the assumption is unlikely to be true globally, we are solely interested in identifying genes where this assumption holds for the purpose of biodosimetry. The model was fitted with the fitlme function from the statistical and machine learning toolbox in MATLAB (www.mathworks.com). P-values of the F-test were first examined. The F-test evaluates the amount of variance explained by the proposed model against the constant model with only ${\beta }_0\ne 0.$A large P-value also suggests there is no linear dose-response effect. For genes with a small F-test P-value ($P<1e-05$), differentially expressed (DE) genes were further identified through strong dose-LFC association with a small P-value of ${\beta }_1$ ($P<1e-5).$ The cutoff $P=1e-5$ is chosen to control the family-wise error rate with the conservative Bonferroni correction (Goeman and Solari 2014). Genes with a positive ${\beta }_1$ have increased expression with larger dose and are therefore upregulated, while those with a negative ${\beta }_1$ are downregulated.

2.5.2. Other statistical software

R packages including VennDiagram and ggplot2 are used to generate the Venn diagram and pathway enrichment figures. Dose-response surface plots are generated with our own script in MATLAB. Enrichment analysis was performed with Enrichr (https://maayanlab.cloud/Enrichr/) and STRING (https://string-db.org/).

3. Results and discussion

3.1. Different genes but similar functional pathways respond to low- and high-dose IR exposure

We identified 354 and 266 DE genes for low (GSE43151) and high (GSE8917) dose IR exposures, respectively, with the P-value of ${\beta }_1$ ($P<1e-5)$ (Supplementary Table 1).

As demonstrated in Figure 1(A), only 25 DE genes are common to both datasets, indicating that for the most part, different genes are involved after low- and high-dose exposures. Note that the low-dose dataset (GSE43151) was based on CD4+ lymphocytes whereas the other datasets were based on whole blood. Transcriptional responses may depend on both the dose of radiation and the type of cells.

Figure 1. Differentially expressed (DE) genes after low- and high-dose IR exposures. A: Venn diagram of the DE genes. B: the percentages of identified DE genes that Overlapped with known functional databases. C: top 5 bioprocesses (biological process), KEGG (Kyoto Encyclopdia of genes and Genomes), and Wiki (Wikipedia) pathways enriched in DE genes after low-dose IR exposure. D: top 5 biological process, KEGG and Wiki pathways enriched in DE genes after high-dose IR exposure.

Figure 1(B) shows percentages of identified genes involved in various known functional databases, including cell cycle (https://cyclebase.org/) (Santos et al. 2015; Riba et al. 2022), immune-related genes (https://www.innatedb.com/), and DNA damage repair (https://repairtoire.genesilico.pl/) (Milanowska et al. 2011). Respectively, 31.2%, 10.9%, and 3.0% of high-dose DE genes are represented in these functional databases, whereas 23.7%, 13.3%, and 3.1%, of the low-dose DE genes overlap with these functional databases. Therefore, different percentages of DE genes are involved in immune functions and cell cycle after low- and high-dose exposures. Besides the dose variations, distinct cell types of the low- and high-dose datasets may also contribute to the differences in DE genes.

Figures 1(C and D) are the functional enrichment analysis for DE genes after low- and high-dose IR exposures, respectively. We report the top 5 pathways, from the BioProcess (biological process), KEGG (Kyoto Encyclopedia of Genes and Genomes), and Wiki (Wikipedia: https://www.wikipathways.org/) databases, that are enriched within the DE gene sets. Interestingly, although only 25 DE genes are common to both gene sets, both groups of responding genes are involved in similar biological functions. For instance, both low- and high-dose DE genes are involved in the biological processes of regulation of programmed cell death and regulation in the apoptotic process. Moreover, both are enriched in 3 KEGG pathways: p53 signalling pathways, Epstein-Bar virus infection, and Chronic myeloid leukemia. They are also involved in 4 common Wiki pathways: p53 transcriptional gene network, genotoxicity pathway, DNA damage response, and miRNA regulation of DNA damage response. KEGG and Wiki are alternative representations of the same biological pathways, so it is not surprising that they overlap substantially. These pathways were also reported in previous investigations (Khan and Wang, 2022; Lumniczky et al., 2021; Huang & Zhou, 2020).

3.2. Distinct DE genes are involved in different phases of the cell cycle

The cell cycle is a process fundamental to cell growth and replication. It is typically divided into 4 phases: G1 (Gap1), S (Synthesis), G2 (Gap 2), and M (Mitosis). IR exposure can cause DNA damage and affect cell cycle progression (Lonati et al. 2021). Cell-cycle genes with known phases were downloaded from Cyclebase (https://cyclebase.org/). We further report DE genes with known phases and checkpoint functions in Figure 2.

Figure 2. DE genes with known cell-cycle phases and checkpoint functions, where genes in red are upregulated with a positive ${\beta }_1,$ and genes in blue are down-Regulated with a negative ${\beta }_1.$ A: Cell cycle related DE genes after high-dose IR exposure (GSE8917). B: Cell cycle associated DE genes after low-dose IR exposure (GSE43151).

Genes in red indicate upregulation and a positive association (${\beta }_1>0)$ between dose and the LFC of expression, whereas genes in blue denote downregulation and a negative correlation between dose and LFC with ${\beta }_1<0.$ One observation from Figure 2 is that there are fewer DE genes for cell cycle with known phases after high-dose (Figure 2(A)) than low-dose (Figure 2(B)) irradiation. For instance, 9 and 4 DE genes are upregulated at the G2M and S1 phases after low-dose IR exposure (Figure 2(B)), but only 3 and 2 DE genes are positively associated with dose during G2M and S1 phases after high-dose exposure (Figure 2(A)). Furthermore, most DE genes during these phases are different and more DE genes are involved during the G2M phase than other phases after both low- and high-dose exposures. Both the dose of radiation and the type of cells may induce distinct DE genes in the two datasets.

Interestingly, although most DE genes at different phases are distinct for low- and high-dose irradiation, there are several DE checkpoint genes common to both. There are 16 checkpoint DE genes (13 upregulated and 3 downregulated) after low-dose IR exposure and 8 checkpoint DE genes (6 upregulated and 2 downregulated) after high-dose exposure. Of these DE genes, 5 out of 6 checkpoint genes upregulated after high-dose IR exposure (BAX, GADD45A, PCNA, PLK3, XPC) are also upregulated after low-dose irradiation. The G1 phase gene BBC3 is also upregulated after both low- and high-dose IR exposures. Gene RNF113A, on the other hand, has a negative association with LFC after a high-dose IR exposure but is positively associated with LFC after a low-dose exposure.

3.3. Persistent transcriptional responses of low- and high-dose IR exposure

DE genes may or may not be persistent over time, depending on the sign of ${\beta }_2.$ To study persistent transcriptional responses after irradiation, we further divide the DE genes into 4 distinct clusters with positive and negative ${\beta }_1$ and ${\beta }_2$ values as demonstrated in Figure 3.

Figures 3(A and B) demonstrate the 4 clusters for high-dose (HD) and low-dose (LD) IR exposures, respectively. The average LFCs in each cluster against both dose and time are plotted in columns 1 and 3 of Figure 3, while the marginal 2D plots of average LFCs vs either dose or time are presented in columns 2 and 4 of Figure 3. Cluster 1 (C1) includes the DE genes with both ${\beta }_1$ and ${\beta }_2$ having positive values ($++).$ Positive ${\beta }_1$ and ${\beta }_2$ indicate that LFCs in C1 are positively associated with both dose and time. Cluster 2 (C2) consists of the DE genes with positive ${\beta }_1$ but negative ${\beta }_2$ ($+-).$ The LFCs in C2 increase with dose but decrease with time course. Cluster 3 (C3) contains the DE genes with negative ${\beta }_1$ but positive ${\beta }_2$ ($-+).$ Genes in C3 are negatively associated with dose but positively associated with time. Cluster 4 (C4) consists of DE genes with both ${\beta }_1$ and ${\beta }_2$ being negative $(--),$ indicating that DE genes in C4 are downregulated over both dose and time. The kinetics of genes in each cluster with normalized dose and time (0 − 1) are reported in Figure 3. As demonstrated in Figure 3, DE genes in C1 and C4 are likely to have a persistent effect, because LFCs are up or down, respectively, over both time course and dose. Genes in C2 and C3, on the other hand, might represent the process of DNA damage and cell recovery over time and only have a short-term transcriptional response, because gene expression is either up over time but down with dose or vice versa. In these clusters, the effects of dose and time on LFCs are in opposite directions, so that LFC might return to zero within a short amount of time. Therefore, DE genes in C1 and C4 merit more attention. Some of these genes might be useful as radiation biomarkers.

Figure 3. Average log-fold changes (LFC) of the DE genes in each cluster over time and dose DE genes are divided into 4 clusters with the sign of ${\beta }_1$ and ${\beta }_2,$ where cluster I (C1) contains DE genes with ${\beta }_1$ and ${\beta }_2$ oth being positive (+ +); cluster 2 (C2) consists of DE genes with positive ${\beta }_1$ but negative ${\beta }_2$ (+ -); cluster 3 (C3) comprises DE genes with negative ${\beta }_1$ but positive ${\beta }_2$ (- +); and cluster 4 (C4) includes DE genes with ${\beta }_1$ and ${\beta }_2$ both being negative (- -). A: the 4 clusters after high-dose (HD) IR exposure. B: the 4 clusters after low-dose (LD) IR exposure. From left to right, columns 1 and 3 are the 3D plots for LFC vs both dose and time course, while columns 2 and 4 are the marginal 2D plots for LFC against dose and time, separately.

The results of further analysis are reported in Figure 4. As demonstrated in Figures 4(A and B), the components of DE genes following high- and low-dose IR exposures are different. Overall, there are 36% upregulated genes (C1 + C2) and 64% downregulated genes (C3 + C4) after high-dose exposure, but 55% upregulated genes (C1 + C2) and 45% downregulated genes (C3 + C4) after low-dose exposure. For genes in clusters C1 and C4 with potentially persistent effect, there are 9% and 35% of DE genes in C1 and C4, respectively, after high-dose irradiation, but 38% and 17% of DE genes in C1 and C4, respectively, after low-dose irradiation. Although transcriptional responses may vary due to different cell types in the two datasets, fewer upregulated and more downregulated genes after high-dose IR exposure might indicate reduced biological activities due to DNA damage and cell death.

Figure 4. Pie charts showing the percentage of DE genes in different clusters and the KEGG pathway enrichment analysis for genes in clusters C1 (${\beta }_1$ > 0 & ${\beta }_2$ > 0) and C4 (${\beta }_1$ < 0 & ${\beta }_2$ < 0) following low- and high-dose IR exposures. LD-C1: Top 5 KEGG (Kyoto Encyclopedia of genes and Genomes) pathways for upregulated genes dose IR exposure; and HD-C1: Top 5 KEGG pathways for upregulated genes in C1 after a high-dose IR exposure. There is no KEGG pathway enriched for downregulated genes in LD-C4.

DE genes with cluster information are reported in Supplementary Tables 1 and 2. The top 5 enriched KEGG pathways for DE genes in C1 and C4 after low- and high-dose exposures are reported in Figure 4(C). The 135 (38%) upregulated genes in cluster 1 (C1) after low-dose irradiation is enriched primarily in the p53 signalling pathway, natural killer (NK) cell-mediated cytotoxicity, and glutathione metabolism, whereas the 25 (9%) upregulated genes in C1 after high-dose irradiation are enriched primarily in the p53 signalling pathway and apoptosis. The 93 (35%) downregulated genes in C4 after high-dose IR exposure are enriched primarily in the NK cell-mediated cytotoxicity, the B cell receptor signalling pathway, and the T cell receptor signalling pathway. There are no enriched KEGG pathways for the 61 (17%) downregulated genes after low-dose exposure.

Table 2. Common DE genes over high- and low-dose IR exposure.

	High-Dose Irradiation (GSE8917)				Low-Dose Irradiation (GSE43151)
Gene ID	${\mathit{\beta }}_1$	${\mathit{\beta }}_2$	${\mathit{P}}_2$	${\mathit{C}}_{\mathit{H}}$	${\mathit{\beta }}_1$	${\mathit{\beta }}_2$	${\mathit{P}}_2$	${\mathit{C}}_{\mathit{L}}$
ARHGEF3	0.449	−0.428	1.72E-06	C2	0.512	−0.089	2.46E-01	C2
BAX	0.317	0.062	4.60E-01	C1	1.287	0.282	3.40E-02	C1
BBC3	0.672	−0.026	8.18E-01	C2	2.187	−0.547	2.55E-03	C2
CCDC109B	−0.644	0.230	3.56E-03	C3	−0.254	0.105	6.18E-02	C3
DCP1B	0.613	−0.670	2.15E-09	C2	0.648	−0.091	1.15E-01	C2
DDB2	1.308	−0.531	2.43E-04	C2	1.394	−0.128	2.44E-01	C2
F11R	0.586	0.211	3.22E-02	C1	0.175	0.138	4.39E-04	C1
GADD45A	1.399	−0.000	9.99E-01	C2	1.086	−0.703	3.25E-11	C2
GSS	0.428	0.060	3.57E-01	C1	0.260	0.240	1.61E-04	C1
PLK3	0.688	−0.079	3.06E-01	C2	0.453	−0.819	3.96E-20	C2
TNFRSF10B	1.085	0.095	5.30E-01	C1	1.444	0.510	6.33E-08	C1
XPC	0.964	−0.360	3.30E-03	C2	0.901	0.042	5.68E-01	C1
CBX3	−0.348	0.054	3.52E-01	C3	0.279	−0.040	5.26E-01	C2
PPP3CC	−0.427	−0.127	1.03E-01	C4	0.176	0.451	2.01E-16	C1
RNF113A	−0.359	−0.139	2.19E-02	C4	0.201	−0.032	4.68E-01	C2
ASCC3	1.176	−0.452	1.22E-02	C2	1.093	0.811	7.58E-10	C1
FBXO22	0.864	−0.542	9.87E-04	C2	0.522	0.253	1.19E-06	C1
FBXW7	0.455	−0.447	5.40E-06	C2	0.306	0.686	1.87E-18	C1
FDXR	1.971	−2.028	1.01E-12	C2	2.826	1.668	3.71E-14	C1
PCNA	1.744	−0.749	2.81E-05	C2	1.338	0.636	5.13E-10	C1
PHPT1	0.908	−1.203	8.07E-14	C2	1.634	0.794	3.89E-08	C1
PPM1D	0.893	−0.340	1.06E-02	C2	0.471	0.187	2.94E-04	C1
REV3L	0.428	−0.254	9.34E-04	C2	0.350	0.302	2.78E-08	C1
SESN1	0.799	−0.650	2.64E-07	C2	0.615	0.377	2.21E-04	C1
SRA1	0.582	0.266	1.37E-02	C1	0.298	−0.204	4.70E-03	C2

The DE genes were selected with $P_1<1e-05$ (not reported in the Table). We report the estimated coefficients ${\mathrm{\beta }}_1$ and ${\beta }_2$ for dose and time, the corresponding P value ($P_2)$ of ${\beta }_2,$ and the clusters each gene belongs to. Genes in yellow are the potential transcriptional markers for dosimetry. The middle 3 genes in orange have opposite dose effects after low- and high-dose exposures. The unshaded genes at the bottom have opposite time effects after low- and high-dose irradiations. The full description of the identified genes can be found at https://www.genecards.org/.

3.4. Biomarker identification

We identified 25 genes common to both datasets (Table 1) with the P value for ${\beta }_1$ ($P_1<1e-05$). Full descriptions of the identified genes can be found at https://www.genecards.org/. The fixed-effect coefficients ${\beta }_1,{\beta }_2,P_2,$ and cluster information ${(C}_H,C_L)$ after high and low-dose IR exposures are reported in Table 2.

These 25 genes have a linear correlation with both time and dose over a wide range of doses and are potential candidate biomarkers for dosimetry. Ideally, a candidate biodosimeter should belong to the same cluster and have the same positive or negative linear association with both dose and time in the two datasets. When a gene has different (positive and negative) associations with time in different datasets, the P-value of ${\beta }_2$ should be insignificant ($P_2>0.05$) in at least one dataset. We identified 12 genes as candidate biodosimeters—ARHGEF3, BAX, BBC3, CCDC109B, DCP1B, DDB2, F11R, GADD45A, GSS, PLK3, TNFRSF10B, and XPC—as reported on the top panel of Table 1 (shaded in yellow) and Supplementary Figure 1. Of these, BAX, F11R, GSS, and TNFRSF10B, which belong to cluster C1, have a persistent effect and might be better candidates for biodosimetry. Another 3 genes (CBX3, PPP3CC, and RNF113A) on the middle panel of Table 1 (shaded in Orange) and in Supplementary Figure 2 have opposite associations with dose after low- and high-dose IR exposures. They are activated with a positive ${\beta }_1$ after low-dose exposure but have reduced activities with a negative ${\beta }_1$ after high-dose exposure. Although these 3 genes are not appropriate for use as dosimetry markers, they might be important for understanding the effects on expression of radiation dose and cell type. The remaining ten genes (ASCC3, FBXO22, FBXW7, FDXR, PCNA, PHPT1, PPM1D, REV3L, SESN1, and SRA1) have opposite associations with time after low- and high-dose irradiation (Table 1 and Supplementary Figure 3). These genes might have nonlinear dose-time interactions and are therefore excluded from the candidate biodosimeters.

3.5. Validation with an independent dataset

We further validate the 12 candidate biomarkers with an independent dataset (GSE23515). There are 10 out of 12 of our candidate genes (BAX, BBC3, DCP1B, DDB2, F11R, GADD45A, GSS, PLK3, and TNFRSF10B) present in this validation dataset; the two other candidate genes (CCDC109B and XPC) were not measured by this microarray platform. Boxplots for the LFC at 3 different doses of irradiation (0.1 Gy, 0.5 Gy, and 2 Gy) are presented in Figure 5.

Figure 5. Boxplots of the LFC of 10 candidate biomarkers following different doses of IR exposure (0.1 Gy, 0.5 Gy, and 2 Gy) in the validation dataset (GSE23515). the y-Axis represents LFC. The red dot and the horizontal line inside each box represent the median.

As shown in Figure 5, average LFCs of all candidate biodosimeters increase with larger doses of IR exposures in the validation dataset. Apart from F11R, these genes have statistically significant monotonic relationships with dose (P < 1e-4, ANOVA with Random Effects) and their average (and median) gene expression increases as the dose becomes larger, indicating that the candidate biodosimeters are reproducible in an independent dataset.

3.6. Advantages of the LME model

Although a considerable number of studies have been conducted on the direct effects of radiation exposure on gene expression in peripheral blood samples (Ostheim et al. 2022), most of them are based on fold changes, statistical tests, and supervised classifications within discrete time points, ignoring the time-dose dynamics of the transcriptional response. However, both the dose and time dependences of gene expression can be important for understanding the molecular mechanisms of radiation response and play a major role in accurate dose reconstruction. More importantly, persisitent effects of high- and low-dose irradiation on gene expression are still not fully understood (Jafer et al. 2020). Therefore, in this pilot study we proposed a LME model to investigate the joint dose-time dynamics of transcriptional responses and to functionally characterize the time- and dose-dependent changes in gene expression. Coefficients of the proposed model aid understanding of the effects of irradiation on gene expression beyond 24 hours after exposure. One advantage of the proposed approach is that DE genes are identified with a strong linear association between dose and LFC of expression, adjusting for time course and eliminating effects of a specific time.

Transcriptional responses are highly dynamic, particularly within 24 hours of IR exposure. We are more interested in genes with persistent changes in expression profiles because such genes would be more useful for dose reconstruction. Such modifications in expression may also be the result of upstream signaling switches or epigenetic variation after IR exposure. We simply divide the DE genes into 4 clusters with the signs of ${\beta }_1$ and ${\beta }_2.$ Genes in clusters C1 and C4 with ${\beta }_1$ and ${\beta }_2$ being either both positive or both negative have persistent transcriptional responses. In particular, reduced expressions of genes in C4 over both dose and time course might indicate dying cells and be the result of irreparable DNA damage caused by IR exposure. Previous studies indicate that most DNA double-strand breaks (DSBs) are repairable, many stress-responsive genes may return to normal over 24 hours, and only a small proportion of DE genes associated with unrepaired DSBs remain differentially expressed beyond 24 hours (Noda 2018; Noda et al. 2012). The present approach helps us to identify genes with transcriptional response persisting beyond one day, which are applicable even if the gene expressions are measured within 24 hours of IR exposure.

3.7. Pathways enriched after low and high dose irradiation

Even though peripheral blood cells are typically considered to be dormant in the body, dysregulations in the cell cycle and checkpoint genes may still have significant implications and can affect cell proliferation, differentiation, immune response, and other processes. In particular, the P53 signalling pathway is activated after both low- and high-dose radiation exposure and plays a key role in cell cycle arrest, DNA damage repair, senescence, and apoptosis. Persistent upregulation of the P53 pathway may indicate unrepaired DNA damage. Interestingly, NK cell-mediated cytotoxicity functions are enriched among the upregulated genes (PPP3CC, ICAM2, TNFRSF10B, FAS, KIR3DL3) after low-dose exposure, but are enriched among downregulated genes (PPP3CC, SH2D1B, GZMB, PIK3R1, CD2) after high-dose exposure, and gene PPPC3CC (Protein Phosphatase 3 Catalytic Subunit Gamma) on this pathway has associations with both time and dose after low-dose irradiations that are in the opposite direction from those after high-dose irradiations (Supplementary Figure 2). NK cells are innate immune lymphocytes that can recognize and destroy abnormal cells, such as tumor or DNA-damaged cells (Paul et al. 2013; Chen et al. 2020). Therefore, the upregulated NK-mediated cytotoxicity pathway indicates that NK cell functions are activated after low-dose IR exposure. On the other hand, enrichment of downregulated genes in NK cell-mediated cytotoxicity, T-cell receptor signalling, and B-cell receptor signalling pathways might indicate that cells undergo permanent arrest of DNA replication or apoptosis after high-dose irradiation and the immune functions of the cells are partially impaired. Our finding is consistent with the fact that lymphoid cells (T, B, and NK cells) are generally more vulnerable to damage from radiation because they divide and proliferate more frequently than other types of blood cells (Heylmann et al. 2021)

3.8. Candidate biomarkers for dosimetry

As mentioned earlier, a good biodosimeter must be IR sensitive, change linearly with dose, have a stable signal without temporal fluctuations (no linear correlation with time), be reliable over a wide dose range, and have no influence of gender, age, or health status. We, therefore, built an LME model and identified 12 genes common to both datasets and having similar correlations with dose and time after high- and low-dose exposures. The 12 genes (ARHGEF3, BAX, BBC3, CCDC109B, DCP1B, DDB2, F11R, GADD45A, GSS, PLK3, TNFRSF10B, and XPC) are mainly involved in cell cycle, DNA repair, and immune functions. Of these, BAX, BBC3, DDB2, GADD45A, TNFRSF10B, and XPC were also reported as candidate biomarkers in previous studies (Lacombe et al. 2018). The 12 candidate biomarkers were also reproducible when validated with an independent dataset. Interestingly, PLK3 belongs to the highly conserved polo-like kinase family, is a key regulator of the cell cycle, and has the typical function of chromosome segregation at G2 phase. Genes in C1 and C4, including BAX, GSS, and TNFRSF10B, might be better biodosimeters as they have persistent transcriptional responses. We precluded 13 other genes—ASCC3, FBXO22, FBXW7, FDXR, PCNA, PHPT1, PPM1D, REV3L, SESN1, SRA1, CBX3, PPP3CC, and RNF113A—from consideration as candidate dosimeters, although several of them (ASCC3, FBXO22, FDXR, PCNA, PHPT1, and SESN1) were reported as potential dosimeters in previous studies that did not consider the joint dose and time dynamics (Lacombe et al. 2018). For instance, FDXR was eliminated from our biomarker list because it has different correlations with time after low- and high-dose exposures, indicating nonlinear dose-time interactions and temporal fluctuations, although it was the gene with the greatest significance in the statistical test. FDXR has been thoroughly studied over the last several years, but a discrepancy in FDXR expression between mice and human cell lines was found and it failed to predict the dose accurately (Abend et al. 2022; Macaeva et al. 2019).

3.9. Limitations of the study

It is challenging to mine gene expression datasets in the public domain. First, besides different cell types, microarray data were collected with different platforms. There are 8639 genes measured after high-dose irradiation, but 13,447 genes are available that were measured after low-dose irradiation; this difference might create certain biases in gene identification. Second, measured time points are different. Maximum measured times are 10 and 24 hours over low- and high-dose irradiations, respectively. This time discrepancy might introduce bias and systematic errors in the analysis. Third, observed gene expression can change when blood is put in a culture or if the RNA sources are different (Paul et al. 2011; Paul et al. 2013). We would also fail to detect marker genes that are associated with both dose and time but not measured in the datasets. Furthermore, other factors such as cell cycle, peri- and post-irradiation conditions, and RNA preparation methods could influence the observed gene expression. Nevertheless, the monotonic association between dose and LFC was confirmed in an independent validation dataset for 9 out of 10 identified candidate biomarkers.

4. Conclusion

We propose an LME model to detect genes that have persistent transcriptional responses or may serve as candidate biological dosimeters with integrated transcriptome analysis. We identified two subsets of genes that are either positively or negatively associated with both dose and time course after exposure and are likely to have persistent expression changes. One especially interesting finding is that the NK cell-mediated cytotoxicity pathway is downregulated after high-dose IR exposure but is upregulated after low-dose IR exposure, and there are more upregulated genes after low-dose irradiation than after high-dose irradiation. Although different types of cells in the low- versus high-dose datasets might impact the finding, we demonstrated that generally cell cycle genes respond to low doses whereas immune suppression mechanisms are induced at high doses. We also identified 12 potential biological dosimeters (ARHGEF3, BAX, BBC3, CCDC109B, DCP1B, DDB2, F11R, GADD45A, GSS, PLK3, TNFRSF10B, and XPC), of which 3 genes—BAX, GSS, and TNFRSF10B—might be the better biomarkers as they have persistent transcriptional responses. We also eliminated 13 other genes that have an opposite association with either time or dose after low- and high-dose irradiations. The findings were validated with an independent dataset. Further validation in a wet lab should be performed to confirm their biological and clinical relevancies.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The Radiation Effects Research Foundation (RERF), Hiroshima and Nagasaki, Japan is a public interest foundation funded by the Japanese Ministry of Health, Labor and Welfare (MHLW) and the US Department of Energy (DOE). SAA received support from the Center for High-Throughput Minimally-Invasive Radiation Biodosimetry, National Institute of Allergy and Infectious Diseases Grant Number U19 AI067773, USA. This research was also funded in part through DOE award DE-HS0000031 to the National Academy of Sciences. The views of the authors do not necessarily reflect those of the two governments.

Notes on contributors

Zhenqiu Liu

Zhenqiu Liu, PhD, is a Senior Scientist at the Department of Statistics, Radiation Effects Research Foundation (RERF), Hiroshima, Japan.

John Cologne

John Cologne, PhD, ELS, is a Senior Scientist at the Department of Statistics, Radiation Effects Research Foundation (RERF), Hiroshima, Japan.

Sally A. Amundson

Sally A. Amundson, PhD, is an Associate Professor of Radiation Oncology at the Center for Radiological Research, Columbia University Irving Medical Center, New York City, NY, USA.

Asao Noda

Asao Noda, PhD, is the Chief of the Department of Molecular Biosciences, Radiation Effects Research Foundation (RERF), Hiroshima, Japan.