Gene expression-based biodosimetry for radiological incidents: assessment of dose and time after radiation exposure

Abstract

Purpose: In order to ensure efficient use of medical resources following a radiological incident, there is an urgent need for high-throughput time-efficient biodosimetry tools. In the present study, we tested the applicability of a gene expression signature for the prediction of exposure dose as well as the time elapsed since irradiation.

Materials and methods: We used whole blood samples from seven healthy volunteers as reference samples (X-ray doses: 0, 25, 50, 100, 500, 1000, and 2000 mGy; time points: 8, 12, 24, 36 and 48 h) and samples from seven other individuals as ‘blind samples’ (20 samples in total).

Results: Gene expression values normalized to the reference gene without normalization to the unexposed controls were sufficient to predict doses with a correlation coefficient between the true and the predicted doses of 0.86. Importantly, we could also classify the samples according to the time since exposure with a correlation coefficient between the true and the predicted time point of 0.96. Because of the dynamic nature of radiation-induced gene expression, this feature will be of critical importance for adequate gene expression-based dose prediction in a real emergency situation. In addition, in this study we also compared different methodologies for RNA extraction available on the market and suggested the one most suitable for emergency situation which does not require on-spot availability of any specific reagents or equipment.

Conclusions: Our results represent an important advancement in the application of gene expression for biodosimetry purposes.

Keywords:

• Biodosimetry
• gene expression
• radiation
• time point prediction
• dose prediction

Introduction

Biodosimetry is the dose estimation after exposure to ionizing radiation by means of changes in biological endpoints, or biomarkers. In the case of large-scale radiological accidents, when physical dosimetry is not available for all the individuals at risk of exposure, these biomarkers could be used to identify individual exposure cases. In such situations triage decisions have to be undertaken as soon as possible in order to split the exposed subjects into different categories, depending on their exposure dose and radiosensitivity (Etherington et al. 2011). This kind of triage will allow focusing the medical staff and facilities only on those subjects in need of urgent medical assistance (Etherington et al. 2011).

As of yet, the gold standard method in biodosimetry is the detection of dicentric chromosomes in peripheral blood lymphocytes (dicentric chromosome assay or DCA) (International Atomic Energy Agency 2011). Besides its sensitivity to doses down to 20 mGy when scoring a few thousands metaphases (Lloyd et al. 1992), this method has many other advantages, such as a high specificity to ionizing radiation, the possibility to detect partial body exposure and the possibility of exposure assessment even months after irradiation (International Atomic Energy Agency 2011). However, DCA is low in throughput: it is time-consuming, laborious and requires well-trained personnel for scoring. In 2010, the total capacity of European Union laboratories specialized in biodosimetry for DCA was estimated to be 1493 samples in the triage mode and 187 samples in the full mode per week, excluding the 48 h needed for lymphocyte culturing (Wojcik et al. 2010). This would be insufficient in case of a large-scale emergency with thousands of potentially irradiated subjects. Several automated systems allowing for faster dicentric scoring have been developed so far (Gruel et al. 2013; De Amicis et al. 2014; Garty et al. 2015), and several studies demonstrated the possibility of scoring a limited number of cells (20–50) for fast triage purposes (Beinke et al. 2013; Oestreicher et al. 2017). A simplification of dicentric scoring can also be achieved with the application of telomere and centromere fluorescence in situ hybridization – the technique which allows to simultaneously stain telomeres and centromeres (M'Kacher et al. 2015). FISH-based translocation analysis is also the only validated biodosimetry method allowing dose assessment years after actual exposure (Hall J et al. 2017), but this does not offer an advantage for emergency triage.

A promising new approach for biodosimetry that offers superior time-efficiency (Abend et al. 2016), is the analysis of changes in gene expression levels. Several genes, which respond to radiation exposure have been studied using different methodological approaches, such as whole genome microarray methods (Dressman et al. 2007; Fachin et al. 2007; Paul and Amundson 2008; Boldt et al. 2012; Knops et al. 2012; Macaeva et al. 2016) or quantitative PCR (Joiner et al. 2011; Kabacik et al. 2011; Manning et al. 2013; Tucker et al. 2014; Brzoska and Kruszewski 2015). Most of the identified genes are known to be regulated by p53 (e.g. MDM2 (Levine et al. 1997), DDB2 (Hwang BJ et al. 1999), FDXR (Hwang PM et al. 2001), PCNA (Xu and Morris 1999), GADD45A (Hollander et al. 1993), RPS27L (He and Sun 2007), SESN1 (Velasco-Miguel et al. 1999)), and are involved in canonical p53-mediated pathways such as cell cycle regulation, DNA damage repair and apoptosis. Because some of them are induced in isolated peripheral blood cells after exposure to doses as low as 5–25 mGy (Knops et al. 2012; Riecke et al. 2012; Manning et al. 2013; Nosel et al. 2013), it can be stated that in terms of sensitivity to low doses gene expression equals the DCA method. As a result, the development of a biodosimetry gene signature and associated assays which can be configured as devices suitable for low-cost, ‘point-of-care’ measurements make an appealing strategy (Brengues et al. 2010; Huang et al. 2011; Joiner et al. 2011).

Recently, we used microarrays to analyze the transcriptional response of human peripheral blood mononuclear cells (PBMCs) exposed ex vivo to radiation doses of 0.0, 0.1, and 1.0 Gy and observed that many differentially expressed genes were also alternatively transcribed/spliced in response to radiation (Macaeva et al. 2016). We thereby identified a signature of genes and exons that showed high performance in dose prediction (Macaeva et al. 2016).

One of the main difficulties in using changes in gene expression as a biomarker of exposure is the highly dynamic and transient nature of the signal. The expression of every single gene following radiation exposure is affected not only by the dose, but also by the time, and the kinetics of expression is different for every gene (Manning et al. 2013). This means that knowledge of the time span between exposure and measurement is pivotal for correct dose prediction. Using a signature of genes rather than one single gene may therefore, allow assessing the time after exposure based on the combination of their expression profiles.

In addition, finding the right methodological approach to monitor gene expression as early as possible following exposure is also important because the signal is lost within days (Hall J et al. 2017). Possible solutions to this include immediate snap freezing of blood in liquid nitrogen or dry ice, which might be challenging in field conditions, or addition of special whole blood preservation buffers (Schwochow et al. 2012), which would also solve the problem of effective preservation of easily-degradable RNA. Another challenge of using whole blood for gene expression studies is the heterogeneity of blood cells. About 99% of blood cells are red blood cells, including immature reticulocytes, which contain high levels of globin mRNA accounting for ∼70% of all mRNA in whole blood. This can compromise the detection of other specific mRNAs from white blood cells (Field et al. 2007). Although qPCR is less affected by globin mRNA contamination, this parameter is highly important for such techniques as microarrays (Liu et al. 2006) and next-generation sequencing (Schwochow et al. 2012).

In the present study we firstly compared RNA quality from blood samples processed using different RNA extraction methods. Next, the potential of a gene expression assay using a customized qRT-PCR array to predict both dose and time after exposure was investigated. Genes were chosen based on our previous study, taking into consideration the different sensitivities of single exons. Besides, we included three genes (PF4, GNG11 and CCR4) which were up-regulated in response to low-dose exposure (0.05 Gy) and down-regulated after a higher dose (1 Gy) (El-Saghire et al. 2013).

Materials and methods

Comparison of RNA extraction methods

As a first step of the present study, the performance of two RNA extraction kits specifically designed for RNA extraction from whole blood – QIAamp RNA Blood Mini Kit (Qiagen, Hilden, Germany) and PAXgene Blood RNA Kit (PreAnalytiX, Hombrechtikon, Switzerland) – was tested under different conditions (Figure 1).

Figure 1. Schematic representation of the experimental procedures.

Peripheral blood samples were collected from five healthy donors with informed consent and ethical approval from the local SCK•CEN Ethics Committee. The procedure was carried out in accordance with the ethical standards of the Helsinki Declaration of 1975, as revised in 2000.

RNA extractions were performed following the manufacturer’s instructions, unless mentioned otherwise (Figure 1). RNA concentration was measured on a Trinean Xpose instrument (Trinean, Gent-Brugge, Belgium) and the quality of total RNA samples was assessed using Agilent 2100 Bioanalyser (Agilent Technologies, Santa Clara, CA, USA) by calculating the RNA integrity number (RIN). RIN for a sample is computed using several characteristics of an RNA electropherogram trace, the most important being the ratio of the area under the 18S and 28S rRNA peaks. Based on the results each sample is assigned a value of 1 to 10, with 10 being the least degraded. Globin mRNA contamination was assessed using qRT-PCR with primers specific to HBA1 and HBB genes using isolated PBMCs as a reference. PBMCs isolation was performed as previously described in Macaeva et al. (2016). RNA extraction from isolated PBMCs was performed using the RNeasy Mini Kit (Qiagen) following the manufacturer’s instructions.cDNA synthesis was performed using the GoScript™ Reverse Transcription System (Promega, Leiden, The Netherlands) with random hexamer primers. For each gene, qRT-PCR reactions were run in duplicate using the MESA GREEN^® qRT-PCR kit (Eurogentec, Seraing, Belgium) on an Applied Biosystems^® 7500 Real-Time PCR instrument following the manufacturer’s instructions. qRT-PCR data were analyzed by 7500 Software v2.0.6 and Microsoft Excel using the Pfaffl method (Pfaffl 2001). The relative amount of transcript of the selected genes was normalized to PGK1 and HPRT1 reference genes using the geometric mean of the threshold cycle (Ct) values of these reference genes (Vandesompele et al. 2002).

Blood collection and in vitro irradiation

Peripheral blood samples were collected in EDTA coated tubes from healthy donors with informed consent and ethical approval from the local SCK•CEN Ethics Committee. All procedures were carried out in accordance with the ethical standards of the Helsinki Declaration of 1975, as revised in 2000. Half of the samples were used as reference samples, the other half was used as blind samples. Donor information can be found in Table 1.

Table 1. Donor information.

Donor	Gender	Age	Smoking status	Recent viral/bacterial infection	Reference/blind sample
D1	Male	24	No	No	Reference
D2	Female	24	No	No	Reference
D3	Male	29	No	No	Reference
D4	Male	23	No	Yes	Reference
D5	Female	28	No	No	Reference
D6	Female	27	No	No	Reference
D7	Female	46	No	No	Reference
D8	Male	56	Recently stopped	No	Blind
D9	Male	58	No	No	Blind
D10	Male	54	No	No	Blind
D11	Male	49	No	No	Blind
D12	Male	40	No	No	Blind
D13	Male	56	Yes	No	Blind
D14	Male	30	No	Yes	Blind

Each blood sample was collected in an EDTA-coated tube and aliquoted in either 5 ml (reference samples) or 2 ml (blind samples) tubes for irradiation. Reference samples were irradiated with 0, 25, 50, 100, 500, 1000, and 2000 mGy, after which the tubes were placed on a rocking platform in an incubator at 37 °C without CO₂ supply. The irradiations were performed at room temperature with an Xstrahl machine (250 kV, 1.4 mm Cu +3.8 mm Al filtration) at a dose rate of 0.14 Gy/min. At 8, 12, 24, 36 and 48 h after irradiation, a 1 ml aliquot was taken from each sample and used for RNA extraction. The doses and fixation time points for blind samples were assigned randomly and are given in Table 2.

Table 2. Information on doses and time points assigned to blind samples.

Donor	Code	Dose, mGy	Time point, h
D8	B	0	8
D8	R	200	8
D8	P	900	12
D9	I	700	8
D9	C	60	12
D9	L	0	24
D10	S	1200	12
D10	A	30	24
D11	W	100	8
D11	F	1600	36
D11	M	0	48
D12	Q	60	36
D12	T	200	36
D12	G	2000	36
D13	Y	400	24
D13	U	400	48
D13	D	1200	48
D14	E	100	12
D14	J	0	24
D14	O	700	48

RNA extraction, quantification and quality control

The QIAamp RNA Blood Mini Kit was used to extract RNA from the blood samples used in the biodosimetry part of the present study. All procedures were performed following the manufacturer’s instructions, unless stated otherwise. RNA concentration was measured on a Trinean Xpose instrument and the quality of total RNA samples was assessed using Agilent 2100 Bioanalyser (Agilent Technologies, Santa Clara, CA, USA). All samples had a RIN >8 and were therefore considered as suitable for further processing.

Reverse transcription and qPCR

cDNA synthesis on the samples used for biodosimetry part of this study was performed using the RT² First Strand Kit (Qiagen) following the manufacturer’s instructions. qPCRs were run using custom RT² Profiler PCR Arrays and the RT² SYBR Green Mastermix (Qiagen) on an ABI7500 Fast instrument (Applied Biosystems). Standard 96-well plate-format arrays were used, including 25 genes of interest, four reference genes (HPRT1, PGK1, GAPDH and B2M), a positive PCR control, a human genomic DNA contamination control and a reverse transcription control. Each 96-well array could therefore be used to run three samples. The list of genes present on the arrays is provided in Supplementary Table 1. Based on our previous results (Macaeva et al. 2016), for the genes alternatively transcribed/spliced in response to irradiation, primers were designed to target the most responsive exons (Supplementary Table 1). qRT-PCR data were analyzed using the dedicated software available at:

http://www.qiagen.com/be/shop/genes-and-pathways/data-analysis-center-overview-page/custom-rt2-pcr-arrays-data-analysis-center/

After comparing the variability of expression (standard deviations of Ct values in 245 reference samples) of the four reference genes, it was decided to use only B2M for normalization. The obtained differences in threshold cycle (ΔCt) values compared to that of B2M were used to train the prediction models. Relative expression levels were tested for statistical significance using a 2-way ANOVA with Bonferroni post-hoc test.

Classification

Two models were built to predict both the radiation dose and time point of each blind sample based on the gene expression profiles. Two separate datasets were used for training and testing. For training the models, the data acquired from seven donors (reference samples, Table 1) were used, while for testing we used the data of seven different donors (blind samples, Table 1). Training both models consisted of five consecutive steps, (1) identifying the genes most likely responsible for the prediction of the radiation doses and time points (i.e. feature identification), (2) selecting the most informative features (i.e. feature selection), (3) building both models using the classifiers available (i.e. model training), (4) testing the models with a database unexposed to the training step (i.e. model testing), (5) selecting the best performing model for both problems (i.e. model nomination).

Feature (gene) identification was performed in our previous study (Macaeva et al. 2016). Both dose and time point were treated as a regression variable ranging from 0 mGy to 2000 mGy and 0 h to 48 h, respectively. Importantly, the testing dataset was not only constructed from new subjects but also new intermediate doses such as 30, 60, 400, 700, 900, 1200 and 1600 mGy not used for training the model were included. This is particularly important to illustrate the independency of the model from the set of training instances.

Secondly, we attempted to perform a feature selection step, important to reduce the complexity and to increase the accuracy of the model. To investigate if all features were necessary for the predictive model, we used the WEKA software implementation to assess the importance of all attributes. Interestingly, however, performing the principal component analysis illustrated the necessity for all features to describe 95% of the variance within the training dataset (data not shown). Consequently, all features were included for training both models using a wide range of classifiers available in WEKA (Hall et al. 2009).

Different classification and regression models were considered including Support Vector Machines (SVM), linear regression, Multi-Layer Perceptron (a neural network), nearest neighbour, and three decision tree models: M5Base (implementing base routines for generating M5 Model trees and rules), decision stump, Fast decision tree learner, as implemented in WEKA (Hall et al. 2009). Next, each of the trained classifiers was tested using the test dataset, during the fourth step (model testing). To select the optimal classifier, we trained and tested various models. We reported for each classifier Pearson correlation coefficient and the relative absolute error (RAE) calculated for dose prediction as following: $${\mathit{R}\mathit{A}\mathit{E}}_{\mathit{D}\mathit{o}\mathit{s}\mathit{e}}={\sum }_{\mathit{i}=1}^{\mathit{n}}|{\mathit{x}}_{\mathit{i}}-{\mathit{y}}_{\mathit{i}}|/{\sum }_{\mathit{i}=1}^{\mathit{n}}|\mathit{z}-{\mathit{y}}_{\mathit{i}}|$$ where $x_i\mathrm{ }$is the predicted dose, $y_i$ is the true dose, $z$ is the average of x values, for every n cases.

For time point prediction RAE was calculated as following: $${\mathit{R}\mathit{A}\mathit{E}}_{\mathit{T}\mathit{i}\mathit{m}\mathit{e}}={\sum }_{\mathit{i}=1}^{\mathit{n}}|{\mathit{x}}_{\mathit{i}}-{\mathit{y}}_{\mathit{i}}|/{\sum }_{\mathit{i}=1}^{\mathit{n}}|\mathit{z}-{\mathit{y}}_{\mathit{i}}|$$ where $x_i\mathrm{ }$is the predicted time point, $y_i$ is the true time point, $z$ is the average of x values, for every n cases.

Finally, we nominated the best performing model for each case (predicting the radiation doses and time points) achieving the lowest mean absolute error whilst having the highest correlation coefficient.

Results

Optimization of RNA extraction protocols for emergency situations

To optimize RNA extraction for emergency situations, different protocols were compared and the total RNA content, integrity and contamination were quantified (Figure 2).

Figure 2. Comparison of RNA extraction protocols. (A) Quantity of RNA extracted following each protocol. Statistical comparison performed using two-tailed paired t-test. In case two bars are mark with the same letter(s) the difference is not statistically significant (p > .05), if the letters are different the difference between two conditions is statistically significant. (B) Quality of RNA extracted following each protocol. Statistical comparison performed using two-tailed paired t-test. In case two bars are mark with the same letter(s) the difference is not statistically significant (p > .05), if the letters are different the difference between two conditions is statistically significant. (C) Levels of contamination of RNA extracted from whole blood samples with PAXgene and QIAmp kits with excessive globin mRNA assessed as expression of α- and β-globin coding genes. Statistical comparison performed using unpaired t-test. SP: standard protocol; MP: modified protocol. Bars represent the mean of five measurement, error bars represent SD. **p < .005, ***p < .001.

In this part of the study a few modifications of the standard protocols of two commercially available kits were tested: PAXgene tubes, which are specifically designed for direct collection of blood and preservation of RNA for up to 3 days at room temperature and up to 50 months at –20 °C and the QIAmp kit intended for molecular biology applications. The PAXgene system showed consistently high yields of high-quality RNA, in case blood was directly collected in PAXgene tubes (Figure 2, Standard PAXgene protocol, Modified protocols 1 and 2). Alternatively, the option of collecting the blood in standard EDTA-coated tubes followed by transfer into PAXgene tubes was also explored (Figure 2, Modified protocol 3). This approach resulted in extraction of high-quality RNA, although at lower yields compared to the standard PAXgene procedure. Importantly, sufficient RNA of acceptable quality could also be extracted from the blood samples which were frozen at –20 °C and transferred into PAXgene tubes after thawing (Figure 2, Modified protocol 4). Since the QIAmp protocol does not allow processing of frozen blood, the options to modify the protocol were limited (Modified protocol 5). Nevertheless, both tested protocols involving this kit resulted in high yields of RNA of high quality (Figure 2, Modified protocol 5 and QIAmp standard protocol). Also the possible contamination of RNA samples with excessive globin mRNA was tested (Figure 2(C)). In this respect, the QIAmp kit clearly outperformed PAXgene, showing low globin mRNA contamination comparable to those obtained for isolated PBMCs, used as a reference in which no contamination was to be expected.

Robust dose- and time-dependent response of the biodosimetric gene panel

In order to identify the best endogenous controls consistently expressed across the sample population, we compared the variability of expression of PGK1, HPRT1, GAPDH and B2M as reference genes. The standard deviations of Ct values for B2M were lower than those for PGK1, HPRT1 and GAPDH (Supplementary Table 2). Also the differences between the minimal and maximal registered Ct values for B2M were the lowest (∼5 cycles difference compared to ∼10 cycles for the other genes). Therefore, in the subsequent analyses, the Ct values of target genes were normalized to B2M.

Blood samples from seven healthy volunteers exposed to X-ray doses ranging from 25 to 2000 mGy were used as reference samples at five different time points following exposure: 8, 12, 24, 36 and 48 hours. In total, expression levels of 25 genes were assessed (Figure 3). The Ct values for all the samples are given in Supplementary Table 3.

Figure 3. (A–E) Heatmaps showing the average of expression levels in seven donors (relative to control samples per time point, log₂ transformed) for all the 25 genes (shown in rows) included in the qPCR arrays in response to exposure to different doses (shown in columns) at five time points. (F) Heatmap showing the relative expression levels in sham-irradiated samples from seven donors relative to sham-irradiated samples at 8 h for all the 25 genes (shown in rows) at five time points (shown in columns).

Table 3. Dose and time point prediction results for blind samples.

Sample code (with remark)	True time point, h	Predicted time point, h	Error, h	True dose, mGy	Predicted dose, mGy	Error, mGy
A	24	23.3	–0.7	30	20	–10
B (ex-smoker)	8	9.5	+1.5	0	20	+20
C	12	10.2	–1.8	60	20	–40
D (smoker)	48	44.0	–4.0	1200	900	–300
E (infection)	12	8.9	–3.1	100	30	–70
F	36	39.6	+3.4	1600	1370	–230
G	36	43.2	+7.2	2000	1370	–630
I	8	4.5	–3.5	700	410	–290
J (infection)	24	22.1	–1.9	0	20	+20
L	24	15.5	–8.5	0	20	+20
M	48	42.7	–5.3	0	120	+120
O (infection)	48	43.9	–4.1	700	900	+200
P (ex-smoker)	12	14.0	+2.0	900	1370	+470
Q	36	39.0	+3.0	60	120	+60
R (ex-smoker)	8	10.5	+2.5	200	420	+220
S	12	18.4	+6.4	1200	420	–780
T	36	40.1	+4.1	200	120	–80
U (smoker)	48	41.0	–7.0	400	900	+500
W	8	10.0	+2.0	100	230	+130
Y (smoker)	24	21.4	–2.6	400	410	+10

Bold: time point predictions with more than 4 h difference from the true time point and dose predictions with more than 500 mGy difference from the true dose.

Overall, we observed dose-dependent increases in gene expression (up to 22-fold in the case of FDXR) for the majority of the genes, at all the different time points (Figure 3). It was clear, however, that different genes had very different kinetics. For instance, induction of FDXR peaked after 8 h (Figure 4(A)), while that of DDB2 was stable up to 36 h after irradiation (Figure 4(B)). TNFRSF10B was stably induced up to 24 h, less so after 36 h and no longer induced after 48 h (Figure 4(C)). Importantly, for many genes (FDXR, DDB2, TNFRSF10B, AEN, XPC, BAX, ASTN2, NDUFAF6, MAMDC4, PHPT1, ASCC3, TRIAP1, LR5) these changes were significant at least at one-time point after exposure to the lowest dose of 25 mGy. We also found, in agreement with previous findings (Manning et al. 2013; Abend et al. 2016), that most genes showed a plateau at doses above 1000 mGy at which no further induction could be observed. However, some genes did not show any significant changes at all, e.g. SESN1, and the genes previously identified as having opposite changes in response to low- and high-dose radiation (PF4, GNG11, CCR4) (El-Saghire et al. 2013) (Figure 3).

Figure 4. (A–C) Changes in expression of FDXR, DDB2 and TNFRSF10B relative to sham-irradiated samples at 8 h time point. Bars represent the average of samples from seven donors irradiated in vitro with doses from 0 to 2000 mGy and fixed at 5 time points after exposure. (D–F) Changes in expression of FDXR, DDB2 and TNFRSF10B relative to sham-irradiated samples per time point. Bars represent the average of samples from seven donors irradiated in vitro with doses from 0 to 2000 mGy and fixed at 5 time points after exposure.

Importantly, a time-dependent up-regulation in gene expression was seen for several genes (Figure 3(F), Figures 4(D–F)), particularly those involved in apoptosis (AEN, BAX, FDXR, TNFRSF10B) and DNA damage repair (ASCC3, DDB2, XPC). This was most evident after 48 h and suggested that even in unirradiated samples, cells were becoming apoptotic over time. Several genes were induced in control samples already after 12 h (Figure 3(F)), which affected the radiation response for these genes at later time points. Nevertheless, we could still observe dose-dependent changes in gene expression, even after 48 h although radiation-induced changes were in general less outspoken at this time point. This may be, at least in part due to the fact that gene expression was already induced in control samples.

Dose and time after exposure can be predicted from radiation-induced changes in gene expression

The ΔCt values of all the reference samples were used as the training dataset to build the models for dose and time point prediction. Next, the selected model was used to predict the dose and the time after exposure of the blind samples (Table 3 and Figures 5 and 6).

Figure 5. The decision tree for the REPTree model built for dose prediction, ending with 11 leaves (for each the dose prediction is calculated).

Figure 6. The decision tree for the M5Base model built for time point prediction, ending with 23 leaves (for each a linear model named LM1-LM23 has been calculated).

For the dose prediction, the fast decision tree learner (implemented as REPTree in WEKA) could achieve the highest Pearson correlation coefficient of 0.86 between the true and the predicted doses with the lowest RAE of 42%. The merit of this classifier is that it builds the decision based on the information gain and prunes the tree branches using reduced-error pruning (with backfitting). Dose prediction was more accurate for doses below 0.5 Gy (Pearson correlation 0.85) compared to higher doses (Pearson correlation 0.55), probably resulting from the previously described plateau effect of the gene expression response at doses above 1 Gy (Manning et al. 2013; Abend et al. 2016; Manning et al. 2017). Although most of the genes were required to explain the variability within the data, a set of four genes (MDM2, FDXR, ASCC3 and CDKN1A) was considered particularly important for dose prediction. For the time point prediction, the M5Base decision tree classifier (implemented in WEKA as M5P) was able to achieve a correlation coefficient between the true and the predicted time points of 0.96 (with lowest RAE of 28%). M5Base implements base routines for the generation of the M5 model (Quinlan 1992; Wang and Witten 1997). It functions as a decision tree by splitting the data into branches and leaves based on few parameters, but instead of assigning the ‘time points’ for each leave, a linear regression function is calculated for each leave (enabling a continuous numeric prediction).

The 48-h time point was the most difficult to predict with four predictions out of four outside of the ±4 h interval (Table 3). The most discriminative genes for time point prediction were AEN, ASCC3, CDKN1A, GNG11 and CCR4, although all genes were necessary to describe the dissimilarity within the data.

In case of both dose and time point prediction using only the most important genes did not result in significantly better model performance (data not shown). All three samples misclassified according to the true dose were also misclassified according to the true time point (Table 3).

Among the blood donors, one was a smoker, one recently stopped smoking and one individual reported a recent infection. Prediction of the doses or time points for these individuals was not particularly more problematic (Figure 7). Also, the median age of donors whose blood was used for reference samples (27 years) was lower than that of the blind sample donors (54 years). Nevertheless dose prediction was still efficient suggesting that age does not affect dose prediction. It is also important to mention that all the blind samples were collected from male donors completely at random. To demonstrate the model’s performance when testing female records, a cross-validation approach for training and testing the model was applied. However, in this case calibration and blind samples were pooled, thereby violating the independence of the two datasets. The performance of this model was very close to the initial one (data not shown).

Figure 7. Dose (A) and time point (B) prediction results for the blind samples. CC: correlation coefficient.

Discussion

Gene expression is emerging as a highly powerful readout for biodosimetry. Several studies demonstrated the applicability of microarray technology for analyzing large sets of transcripts for dose prediction for exposure to acute doses (Dressman et al. 2007; Meadows et al. 2008; Paul and Amundson 2008; Macaeva et al. 2016). However, as was demonstrated in a recent biodosimetry exercise, the analysis of a short list of genes, or even one gene by means of qPCR technology, is far more straightforward and cheaper compared to microarray analysis, while the accuracy of dose prediction is essentially similar (Abend et al. 2016). Many of the genes included in our biodosimetry panel were also previously tested by other groups (Riecke et al. 2012; Manning et al. 2013; Tucker et al. 2014; Brzoska and Kruszewski 2015), suggesting their robustness as radiation biomarkers. However, our panel composition is unique because we selected the most radiation-sensitive exons of the respective genes for qRT-PCR analysis based on our previous results (Macaeva et al. 2016). We used seven X-ray doses ranging from 25 mGy to 2 Gy to generate reference samples. Moderately high doses of 1–2 Gy may result in acute radiation syndrome, therefore the individuals receiving doses above this threshold will benefit greatly from prompt treatment (Sullivan et al. 2013). However, in case of a radiological accident, most victims would be exposed to much lower doses. Although lower doses down to 100 mGy are not associated with immediate health effects they might require medical follow-up due to the risk of long-term effects, particularly cancer (UNSCEAR 2006). To add to that, the fears about the effects of radiation and the uncertainty about the exposure doses may lead to further stress and mental health effects, as has already happened after the Chernobyl and Fukushima accidents (WHO 2006; Brumfiel 2013; Bromet 2014). In this respect, we believe that assessment and communication of any radiation dose is important. As it is unlikely that any single biodosimetry assay can be used as a stand-alone tool to meet the requirements of throughput, time-efficiency and accuracy needed, it was suggested to perform stepwise triage starting with assessment of less than or greater than a 2-Gy dose, followed by a secondary high-throughput screening to further define individual doses (Sullivan et al. 2013). Thus, a combination of high-throughput methods with high accuracy in the low (e.g. gene expression), and the high dose ranges could prove to be an ideal match for practical biodosimetry in emergency cases.

An important aspect that we have addressed in this study is the dynamic nature of the gene expression response to radiation exposure which can have a major impact on accurate dose prediction in cases when the exact time point of irradiation is not known. Taking the example of FDXR, a gene often considered to be the most sensitive and accurate gene biodosimeter (Abend et al. 2016; Hall et al. 2017; Manning et al. 2017; O’Brien et al. 2018). In our experimental set-up, a ∼5-fold up-regulation of FDXR could indicate (i) the early (8–24 h), maximal response to a low dose of 100 mGy, (ii) the very early (<8 h) sub-maximal response to a higher dose, or (iii) the late (>36 h) residual response to a higher dose. This example shows that, no matter how sensitive, a single gene cannot sufficiently predict the dose when there is no indication of the time elapsed since exposure. In such a case, the combination of a panel of genes, each with their own kinetics and amplitudes is needed. The same example also highlights that is pivotal to take the time elapsed after exposure into consideration. The classification method used in our study allowed to classify the samples according to the time elapsed since exposure with high precision (errors ≤4 h). The importance of the time factor is also highlighted by the fact that all three samples misclassified according to the true dose were also misclassified according to the true time point. It is important to note, however, that the gene expression kinetics might be different following in vivo exposure. Thus, adaptation of our methodology for in vivo exposure situation may be required. The possibility of using specific weighting factors for adaptation of in vitro generated calibration curves for dose prediction of in vivo exposed blood samples from radiotherapy patients was successfully explored in a recent biodosimetry exercise (Abend et al. 2016). The potential of FDXR as radiation biomarker of in vivo exposures was thoroughly investigated in a recent study by O’Brien and co-authors (O'Brien et al. 2018). This study demonstrated that the expression of this gene was significantly up-regulated already at 2 h following a diagnostic CT-scan (low-dose exposure to 2.9–20.9 mGy dose to the blood) (O'Brien et al. 2018). The expression of FDXR was also significantly up-regulated 24 h following the first fraction of radiotherapy treatment in most patients involved in this study and its expression remained up-regulated during the fractionated treatment (O'Brien et al. 2018). Importantly, no significant differences in FDXR expression were found between ex vivo and in vivo irradiated samples from the same patients (O'Brien et al. 2018).

In many studies published so far, fold changes in gene expression between irradiated and control samples from the same subject were used (Kabacik et al. 2011; Riecke et al. 2012; Manning et al. 2013). In a real large-scale emergency, such data will not be available. Therefore, in the present study, we opted to use ΔCt values compared to the reference gene. This approach has been successfully applied in several exploratory studies (Tucker et al. 2014; Brzoska and Kruszewski 2015), as well as in two recent biodosimetry exercises (Abend et al. 2016; Manning et al. 2017).

For most of the studied genes we found a plateau in expression for high doses (1000 and 2000 mGy), which contributed to less precise dose prediction in the higher dose range as exemplified by the lower correlation coefficient (CC = 0.55) compared to that of the lower dose range (CC = 0.85). This may reflect a saturation of the response, but it might also be the result of in vitro culture conditions. Indeed, similar findings have been previously observed in other studies using ex vivo irradiated whole blood or cultured blood cells (Manning et al. 2013; Ghandhi et al. 2015; Abend et al. 2016; Manning et al. 2017). In contrast, in experiments using ex vivo irradiated blood which was diluted tenfold in RPMI culture medium, this effect was only observed at doses above 6 Gy (Tucker et al. 2014), comparable to what was seen in whole body irradiated mice in a study from the same group (Tucker et al. 2013). Similarly, in radiotherapy patients undergoing total body irradiation several genes included in our panel (e.g. CDKN1A, FDXR, and PHPT1) showed further up-regulation after 3.75 Gy compared to 1.25 Gy (Paul et al. 2011).

In the present study, we opted for using undiluted blood for our experiments, to interfere as little as possible with the samples as gene expression changes might be affected by culture conditions, and more specifically by addition of fetal bovine serum (Shahdadfar et al. 2005; Bieback et al. 2010). A similar approach as ours was used by Manning and co-authors, showing 99% cell viability in blood samples incubated for 24 h at 37 °C (Manning et al. 2013). In our study, an increase in expression with time, especially at 36 and 48 h was seen for a number of genes involved in apoptosis, suggesting an increase in the number of dying cells with time. This may have also influenced the radiation-induced response of these genes, and therefore the dose prediction. Nevertheless, we were able to accurately (absolute error ≤300 mGy) predict the radiation doses of 6 out of 8 blind samples from the 36-h and 48-h time points.

In vivo irradiation experiments performed by Tucker and co-authors in mice confirmed the possibility of dose prediction based on gene expression for at least 7 days after exposure (Tucker et al. 2013). Whether the same is true in humans remains to be investigated. Studies involving radiotherapy patients showed that prediction of in vivo radiation dose using gene signatures was possible for at least 24 h following exposure (Paul et al. 2011; Abend et al. 2016). Further in vivo validation of biodosimetric gene signatures at longer time intervals is therefore warranted.

In our study, we did not address the effect of dose rate on gene expression. Several studies addressing this issue were previously performed in total body irradiated mice (Paul et al. 2015) and ex vivo irradiated human blood (3.1 mGy/min vs 1.03 Gy/min) (Ghandhi et al. 2015). Overall, these studies showed that a significant number of genes responded similarly to low dose rate and acute exposures. This was especially seen for typical P53-regulated genes, many of which were also included in our gene panel, although the amplitude of the response was sometimes higher after high dose rate exposure. Interestingly, in both studies classification models could very efficiently discriminate between low and high dose rate irradiated samples, showing that gene signatures may be used to discriminate between acute and protracted exposures.

An important step for biodosimetric triage of casualties of a radiological accident is to appropriately collect and preserve (blood) samples for further analysis. For gene expression measurement the time elapsed from exposure to blood collection and sample fixation is critical for correct dose estimation, as gene expression in non-frozen and non-stabilized blood would still be subjected to change not only as a result of radiation exposure, but also due to ex vivo incubation itself (Baechler et al. 2004). Several approaches to solve this issue were previously tested. Brzoska and Kruszewski extracted RNA for gene expression analysis from whole blood which was frozen and stored at –75 °C (Brzoska and Kruszewski 2015). In a real-life situation this approach would, however, not be practical due to the lack of very low-temperature freezers. An alternative approach, tested during a recent RENEB biodosimetry exercise, includes the addition of RNA stabilization reagents which protect RNA from degradation allowing for sample storage and transportation at room temperature for several hours or even days. However, this methodology still requires the availability of significant quantities of such reagents at the accident site while the further processing of the samples might be affected by transport conditions (Manning et al. 2017). Therefore, we tested the freezing of whole blood at –20 °C, which would allow blood preservation without any specific equipment or reagents. Upon delivery to the biodosimetric laboratory, blood samples would be thawed and transferred to PAXgene tubes, allowing for recovery of sufficient RNA of acceptable quality. A recent biodosimetry exercise showed that dose prediction accuracy using our qRT-PCR arrays is not affected by the differences in RNA extraction methods, including those which result in high globin mRNA contamination (Manning et al. 2017). Thus, in case of significant damage to infrastructure after a large-scale accident, the suggested methodology of blood freezing at –20 °C could be substituted with RNA stabilization reagents as described above.

Taken together, in the present study, we proved the usefulness of a biodosimetric panel of genes in predicting both dose and time after in vitro exposure. Our results confirm that the analysis of expression of these genes, which can be carried out in virtually any laboratory possessing a qPCR instrument, can certainly provide sufficient information for triage purposes in comparatively short amount of time. The validity of our customized biodosimetric qPCR arrays was recently confirmed in an interlaboratory comparison exercise organized by RENEB (Manning et al. 2017). Although in this exercise different experimental procedures, including different dose assessment approach were used, we could still achieve good dose prediction (Manning et al. 2017), confirming the flexibility and versatility of our method. We also suggest a method allowing for further RNA extraction from blood samples stored at –20 °C, which is the most realistic for field conditions. Further research is needed to confirm the time- and dose-wise validity and applicability of our signature for in vivo situation, as well as to reveal the potential impact of confounding factors on the reliability of the gene expression-based biodosimetry.

Acknowledgments

The authors wish to thank the SCK•CEN medical service staff for the help with blood drawing, as well as all the blood donors involved in this study. Ellina Macaeva is a recipient of a joint doctoral SCK•CEN/UGent grant.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This study was supported by the EU-FP7 project DoReMi [grant 249689]. Ellina Macaeva was a recipient of a SCKCEN/UGent PhD grant.

Notes on contributors

Ellina Macaeva

Ellina Macaeva is a PhD student at the Radiobiology Unit of SCK•CEN working on transcriptional radiation biomarkers.

Mohamed Mysara

Mohamed Mysara is a bioinformatician at the Belgian Nuclear Research Centre.

Winnok H. De Vos

Winnok De Vos (PhD Bioscience Engineering) is research professor at Antwerp University and guest professor at Ghent University. He specialises in cell-centered research, with particular focus on the role of the cell nucleus in accelerated aging diseases. To this end, he relies on a multidisciplinary cocktail of techniques, including genome editing and deep coverage microscopy.

Sarah Baatout

Sarah Baatout is head of the Radiobiology Unit at the Belgian Nuclear Research Centre, guest professor at Gent University and maître de conférence at University of Namur.

Roel Quintens

Roel Quintens (PhD) is a senior scientist in the Radiobiology Unit at the Belgian Nuclear Research Centre. Dr. Quintens coordinates two research lines: (i) effects of prenatal irradiation on embryonic brain development, and (ii) identification of novel transcriptional radiation biomarkers.

Gene expression-based biodosimetry for radiological incidents: assessment of dose and time after radiation exposure.