To better predict clinical outcome after radiation exposure, it is very important to know the absorbed dose and body areas exposed. Previously we found that 22 miRNAs appeared to predict total- and partial-body irradiation (TBI and PBI, respectively) patterns and were suggestive of the percentage of the body exposed in a baboon model. Motivated by these results, we performed a similar analysis on the transcriptional level (mRNAs) using whole genome microarrays. From 17 irradiated baboons, blood samples were taken before, and at 1, 2, 7, 28 and 75–106 days postirradiation to an equivalent TBI dose of 2.5 or 5 Gy applied either to the total body or to different parts of the body such as the upper body (UBE) or left hemibody (LHB). We compared quantile normalized log2-transformed gene expression values with three exposure pattern comparisons, namely TBI vs. PBI, TBI vs. LHB and UBE vs. LHB using Kruskal-Wallis and logistic regression analysis for receiver-operator characteristic (ROC) calculation. We found several hundred significantly (P < 0.05) and ≥2-fold deregulated mRNAs per exposure pattern comparison with a peak of 163–860 mRNAs at day 28. Lower numbers on day 2 (60 mRNAs) and day 7 (91–162 mRNAs) were observed, with the lowest number of deregulated mRNAs at day 75–106 (22–58 mRNAs). The 14 most promising mRNAs (e.g., LTF, DEFA3, OLFM4) appeared 10.1–46.2-fold upregulated and the exposure groups were completely or almost completely discriminated (ROC between 0.8–1.0). Several of the mRNA gene expression changes were significantly associated with the percentage of the body exposed. The numbers of overlapping genes used for diagnosis on consecutive days postirradiation were mostly 0 or less than 10. Bioinformatic analysis confirmed that at each time point different biological processes predominated. Our results suggest mRNA changes over time may be used to retrospectively determine radiation exposure patterns as partial or total body. mRNA gene expression changes likely could be applied over a longer time frame (2–75 days postirradiation) than miRNA, but due to the transient gene expression changes a different set of candidate mRNAs appears to be required at each day after irradiation.
After 2 Gy total-body irradiation (TBI) and higher in humans, life-threatening acute health effects develop, referred to as acute radiation syndrome/sickness (ARS) (1 ). The pathomechanism is rapid cell death (reproductive, interphase, apoptotic or necrotic) leading to a functional organ and organ system deficit and a systemic whole-body inflammation (1, 2 ). Several features of radiation such as absorbed dose, dose rate, radiation quality, fractionation, and whole- or partial-body exposure affect cell survival and, ultimately, clinical outcome (3–5 ). For instance, a local arm exposure with 4 Gy will reduce the number of bone marrow stem cells in the irradiated limb dramatically, but hematopoietic syndrome occurrence is unlikely, since bone marrow stem cells in other body areas will compensate for any potential reduction in the number of peripheral blood cells. In contrast, after 4 Gy whole-body exposure, the total number of bone marrow stem cells will decrease, with subsequent reduction in lymphocytes, granulocytes and thrombocytes, resulting in life-threatening conditions including immune suppression and hemorrhage. Therefore, clinical outcome after exposure hinges on both total absorbed dose and exposed body area. While possible radiation exposure scenarios include total-body irradiation, a large percentage of irradiation victims in a mass casualty scenario will have heterogenous exposure patterns of partial-body irradiation (PBI) (6 ). For instance, shielding by building structures resulted in several partial-body exposure patterns after the atomic bomb in Hiroshima and Nagasaki (7–9 ), whereas total-body exposures are predominantly reported from the Chernobyl nuclear power plant accident (10 ). With regard to biodosimetry in a scenario with heterogenous exposure patterns, a previously published murine study suggested that the use of multiple gene expression signatures may help to distinguish the anatomic site of radiation exposure (11 ). Further published findings indicate that a novel combination of radiation-responsive biomarker proteins also shows the potential utility of applying dose prediction models of TBI to PBI (12 ).
In collaboration with the French Army Biomedical Research Institute, we sought to elucidate the physiologic effects of absorbed dose and body exposure pattern. Blood samples from irradiated baboons were taken before (day 0) and at 1, 2, 7, 28 and 75–106 days after either total-body, left hemibody (LHB), upper-body (UBE) irradiations, or irradiation to other parts of the body, of 2.5 or 5 Gy. The doses were not primarily lethal, but sufficient to induce hematological ARS by simultaneously avoiding severe gastrointestinal (GI) ARS. In general, an observed hematological ARS is known to start from 2–2.5 Gy TBI, and GI ARS from 8–10 Gy TBI or 12 Gy PBI delivered to GI tract (4 ). The clinical observations regarding these baboons were reported previously (4, 13 ). In previously published work concerning gene expression, we reported that several miRNAs (n = 22) provided indications of exposure patterns and a suggestion of the percentage of total-body area exposed (14 ). Motivated by these results, we performed a corresponding analysis on the transcriptional level (mRNAs) employing whole genome microarrays in a nonhuman primate model. The goal of this work is to enable health status assessment of an individual receiving partial-body and total-body irradiation.
MATERIALS AND METHODS
All applicable international, national and/or institutional guidelines for the care and humane treatment of animals were followed. All procedures performed in studies involving animals were in accordance with the ethical standards of the institution or practice where the studies were conducted.
Eighteen male baboons were bred by the Centre National de la Recherche Scientifique (Rousset sur Arc, France) for the purpose of biomedical research. In the nonhuman primate facility of the French Army Biomedical Research Institute, the baboons were placed in individual cages at 21°C, with a relative humidity of 55% and a 12:12 h light-dark schedule. The animals received fresh fruit and solid food twice a day, and had access to water ad libitum. The baboons had an average age of 8.1 years (±3.3 years) and weighed 23.7 (±5.2 kg). The experiment was approved by the French Army Animal Ethics Committee (no. 2010/12.0). All baboons were treated in compliance with the European legislation related to humane animal care and protection. Nevertheless, due to unusual blood cell counts before irradiation and his sudden death after irradiation, one baboon was excluded, leaving 17 baboons eligible for analysis.
Irradiation and Exposure Pattern
The animals were anesthetized with a combination of tiletamine and zolazepam (6 mg.kg–1 intramuscularly, Zoletil 100; Virbac, Carros, France) before irradiation. Then, the baboons were placed in restraint chairs, sitting orthogonally, front to a horizontal and homogeneous field of gamma rays delivered by a 60Co source (IRDI 4000; Alsthom, Levallois, France) to perform either TBI or PBI. To attain different patterns of PBI, a 20 cm thick lead screen was used to shield different parts of the body as detailed in Table 1. Two baboons received 5 Gy TBI and two others received 2.5 Gy TBI. Nine different exposure patterns were simulated and up to two baboons were exposed per pattern, summing to 12 baboons receiving PBI [details are provided elsewhere (13 ) and in Table 1] corresponding to an equivalent TBI dose of 2.5 or 5 Gy. Two dose rates were used (8 cGy/min for 5 Gy TBI and 5 Gy 50% PBI, and 32 cGy/min for all other situations), because the 60Co source was changed during this study. Moreover, to achieve the same homogeneous radiation field whatever the dose rate, all baboons were irradiated at the same distance from the source. Consequently, radiation exposures lasted between 8 min and 62 min. The mid-line tissue (right anterior iliac crest) dose in air was measured with an ionization chamber. Delivered doses were controlled by alumina powder thermoluminescent dosimeters placed on different cutaneous areas [thorax, thoracic and lumbar vertebrae, head, tibia, femur, femoral head; for details see (13 )]. A mathematical reconstruction for dosimetry of organ doses was not performed.
Based on the exposure pattern and number of animals, five groups were defined (Table 1):
total-body exposure with either 2.5, 5 Gy or 7.5/2.5 Gy (per hemibody, n = 5);
upper-body exposure with 5.55 or 6.25 Gy with one or both legs shielded (n = 5);
total-body exposure, but head and neck shielded, with 6.25 Gy (n = 1);
left hemibody exposure with 5 or 10 Gy only (n = 4);
only head and arms exposure with 15 Gy (n = 2).
Partial-body exposure is a combination of groups 2, 3, 4 and 5 (n = 12). In addition to these exposure pattern categories, we constructed another variable to reflect the percentage of the body area that was irradiated (Table 1, right-side column).
RNA Extraction and Quality Control
Whole blood samples (2.5 ml) were processed following the PAXgene™ Blood RNA system (BD Diagnostics, PreAnalytiX GmbH, Hombrechtikon, Switzerland). In brief, blood was drawn into a PAXgene Blood RNA tube at the French Army Biomedical Research Institute. The tube was gently inverted (10 times), stored at room temperature overnight then at –20°C. After all samples were collected, the PaxGene tubes were sent to Germany for further processing. After thawing, washing and centrifugation, cells in the supernatant were lysed (Proteinase K) followed by addition of lysis/binding solution taken from the mirVana™ Kit (Life Technologies, Darmstadt, Germany). With the mirVana kit total RNA, including small RNA species, was isolated by combining a phenol-chloroform RNA precipitation with further processing using a silica membrane. After several washing procedures, DNA residuals were digested on the membrane (RNAse free DNAse Set, Qiagen, Hilden, Germany). RNA was eluted in a collection tube and frozen at –20°C. Quality and quantity of isolated total RNA were measured spectrophotometrically (NanoDrop, PeqLab Biotechnology, Erlangen, Germany). RNA integrity was assessed using the 2100 Agilent Bioanalyzer (Life Science Group, Penzberg, Germany) and DNA contamination was controlled by conventional PCR using an actin primer. We used only RNA specimens with a ratio of A260/A280 ≥ 2.0 (Nanodrop) and RNA integrity number (RIN) ≥ 7.3 for qRT-PCR analysis. From 2.5 ml whole blood we isolated between 6–12 µg total RNA. RNA integrity (RIN) with a mean value of 8.6 [standard deviation (SD) ± 0.6, min 7.3, max 9.5] suggested high quality RNA.
Whole Genome Microarrays
Whole genome analysis for differentially expressed genes (protein coding mRNAs) was performed on 102 RNA samples originating from 17 baboons examined before, and at 1, 2, 7, 28 and 75–106 days postirradiation. We used the Agilent oligo microarray GE 8x60K (Agilent Technologies, Waldbronn, Germany) combined with a one-color based hybridization protocol of GeneSpring GX12 software for data analysis as described in detail elsewhere (15 ). Gene expression was analyzed using quantile normalized log2-transformed probe signals as an outcome. The non-parametric Kruskal-Wallis (KW) test was used to compare gene expression across different TBI and PBI exposure patterns as described below (in Statistical Analysis). Only those gene transcripts that had a call “present” in at least 50% of RNA specimens were included in the analysis of gene expression and only genes with KW P values ≤ 0.05 and a ≥2-fold gene expression difference among compared groups were considered. Due to the explorative nature of this study and the low sample size, we did not correct for multiple comparisons, but considered this within our bioinformatic approach and the linear regression analysis (see below). Of note, in previously published work, we successfully validated our microarray data by using qRT-PCR in an independent previously unused sample set (16 ), lending credibility to our approach.
All genes associated with P values ≤ 0.05 and a ≥2-fold gene expression difference (up or down relative to the reference) underwent gene set enrichment analysis using PANTHER pathway software version 15.0 (http://www.pantherdb.org/). PANTHER groups genes with similar biological function based on their annotation (reference list was the current Homo sapiens GO database). For these P values, we corrected for multiple testing by employing the false discovery rate algorithm (17, 18 ).
During analysis of gene expression changes between TBI and PBI exposure patterns over time, in most instances the PBI group was used as the reference. Following the exposure groups as shown in Table 1, we performed three group comparisons and a percentage of body area irradiated was used as a continuous variable:
TBI (n = 5) vs. PBI group I (using all remaining samples and combining them into one group, n = 12). We refer to this comparison as “TBI vs. PBI”.
TBI (n = 5) vs. left hemibody exposure group (n = 4). We refer to this comparison as “TBI vs. LHB”.
Upper body (n = 5) vs. left hemibody exposure group (n = 4). We refer to this comparison as “UBE vs. LHB”.
Gene expression changes over the percentage of the irradiated body area (Table 1, right-side column, n = 17).
Comparisons for binary groups 1–3 were performed using a nonparametric (KW) test and logistic regression for receiver-operator characteristic (ROC) calculations. Exposure groups represented the outcome variables and quantile normalized gene expression values the explanatory variable. Using logistic regression models we calculated odds ratios (OR), 95% confidence intervals (95% CI) and corresponding P values (Wald chi-square) as well as the area under a ROC curve to assess diagnostic accuracy. ROC areas of 1.0 indicate a complete distinction between group comparisons. The fourth comparison was performed using univariate linear regression analysis and employed forward and backward selection procedures on multivariate linear regression models. We assessed the assumptions of normality (Kolmogorov-Smirnov) and, if necessary, we log- or boxcoxtransformed mRNA expression values. Finally, for these RNA species, we used the Bonferroni method (19 ) to correct P values for multiple comparisons. Descriptive statistics (n, mean, standard deviation, minimum, maximum) and P values (non-parametric KW) were calculated for each of the mRNA species and per time point. All calculations were performed using SAS (release 9.4; Cary NC). Graphical presentations were performed using either SAS or SigmaPlot™ 14 (Jandel Scientific, Erkrath, Germany).
Data Availability Statement
The data sets generated or analyzed during the current study are available from the corresponding author on reasonable request.
Radiation Exposure Patterns and Number of Deregulated mRNA Species Over Time
Irrespective of the three group comparisons, we observed a wave-like increase in the number of deregulated mRNA species with a peak at 28 days postirradiation and a decline thereafter (Fig. 1). For instance, for the comparison TBI vs. PBI, the number of deregulated mRNA species increased from 28, 60 and 162 to 253 at days 1, 2, 7 and 28 postirradiation and decreased to 22 mRNA species at 75–106 days postirradiation (Fig. 1A and Table 2). A similar pattern was observed for the TBI vs. LHB comparison, but the number of deregulated mRNA species was >50 as early as the first day postirradiation (Fig. 1B and Table 2). Also, most of the deregulated mRNAs in the TBI vs. PBI comparison were upregulated, but for TBI vs. LHB, the number of up- and downregulated genes appeared similar. Regarding the UBE vs. LHB comparison, we observed a wave-like increase and decrease of up- and downregulated mRNAs (approximately the same proportion) with a peak of 860 mRNA species at day 28 postirradiation (Fig. 1C, Table 2). For intercomparison, we juxtaposed corresponding post-transcriptional miRNA data from previous studies that also showed a wave-like pattern, but the peak occurred at day 7 postirradiation (Fig. 1D).
Exposure Patterns and Overlap in Number of Deregulated mRNA Species Over Time
TBI vs. PBI comparisons at the five different time points revealed virtually no overlap of mRNA species deregulated over more than one time point (Fig. 2A). For a few time points an overlap, in which the same mRNA species were deregulated at multiple time points, was found, e.g., between days 2 and 7 (n = 2) or between days 7 and 28 (n = 10 and Fig. 2B and C). These features were observed for the TBI vs. LHB and UBE vs. LHB group comparisons as well (Fig. 2C). Compression of data over time revealed a high number of overlapping deregulated genes (n = 112) between TBI vs. PBI (n = 331) and TBI vs. LHB comparisons (n = 257; Fig. 2D). This is of minor significance regarding the UBE vs. LHB comparison with the other two comparisons.
Radiation Exposure Pattern and Magnitude of Gene Expression Changes Over Time
The majority (78–80%) of mRNA species appeared 2-to-3-fold up- or downregulated in all group comparisons (Table 3). The number of genes being 3–5 or >5-fold deregulated (13–17%) reached 12–107 at days 7 and 28 postirradiation and, in general, a low number of genes below 10 was found for the other time points. Examination of the most promising 14 mRNA species, based on the fold change and P value, revealed mRNAs that were 10.1-to-46.2-fold upregulated [corresponding to threshold cycle (Ct), differences ranging between 3.3–4.4] and the binary exposure group comparisons could be either completely or almost completely discriminated (ROC between 0.8–1.0; Fig. 3 and Supplementary Table S1; https://doi.org/10.1667/RADE-20-00121.1.S1).
We assessed whether the percentage of the exposed body area (Table 1) was also associated with gene expression changes. We found several mRNA species with gene expression changes significantly associated with the percentage of the body area exposed to radiation (Supplementary Table S1; https://doi.org/10.1667/RADE-20-00121.1.S1). No combination of examined genes resulted in a statically significant model and increased explained variance. Most of these mRNAs (80%) remained significant after Bonferroni correction for multiple comparisons (see also P values in Supplementary Table S1). Strong associations of gene expression (e.g., P < 0.001 for LTF) with the percentage of the exposed body area were significant at multiple days after irradiation as well as day 28 (data not shown).
Bioinformatic Analysis Over Time
The comparison between TBI and PBI groups showed no significant enrichment of genes coding for biological processes at days 1 and 2 postirradiation (Supplementary Table S2; https://doi.org/10.1667/RADE-20-00121.1.S2). At days 7 and 28 postirradiation, reactome-related processes regarding “neutrophil degranulation” and the “immune system” appeared significantly enriched. The enrichment from day 7–28 was at times increased twofold and corresponding lower P values were observed at day 28 relative to day 7 postirradiation (Supplementary Table S2). Except for these processes in general, no further overlapping biological processes could be identified at the time points of 7, 28 and 106 days postirradiation.
Many different characteristics related to radiation exposure must be considered for a meaningful prediction of acute radiation health effects (3 ). In addition to absorbed dose, the exposure pattern of the body is also a critical determinant of ARS. For instance, during the Georgian accident, local skin exposures of up to 150 Gy led to acute injuries (ulcerations), but individuals survived (20 ). Such would not be the case with a whole-body dose of the same magnitude.
Measuring gene expression changes in peripheral blood represents a fast, early and high-throughput methodology in easily accessible biosamples (21–23 ). Taking advantage of a non-human primate model with high genetic homology to humans, we sought to discriminate different exposure patterns by assessing changes on the transcriptional level (protein coding mRNAs) in the peripheral blood over the days, weeks and months after irradiation. Altogether, 17 baboons were irradiated, reflecting five total-body and 12 partial-body exposure patterns. Along with these categories, we also examined gene expression changes corresponding to the percentage of body area exposed. Hereby, we focused more on qualitative differences and, if differences in gene expression can be detected, whether different body areas are exposed to ionizing radiation.
In previously published work we evaluated the post-transcriptome and identified 22 miRNAs (mostly at day 7 postirradiation) which appeared to indicate the exposure pattern and possibly the percentage of the exposed body area (14 ). Our current examination on the transcriptional/mRNA level complements and extends the previous findings. We have now identified hundreds of mRNA species that also showed a wave-like increase in the number of deregulated mRNAs as were found for miRNAs in previous studies on the same baboons. However, the number of deregulated miRNAs in previous studies peaked at day 7 postirradiation and for mRNA this peak shifted to day 28 postirradiation, with the number of deregulated genes increased by an order of magnitude (Fig. 1). This pattern in time supports a causal miRNA-mRNA association, which is very well established in the literature (24, 25 ).
The large number of mRNAs that were significantly associated with a certain exposure pattern encouraged the notion that a set of mRNAs could be used independent of the time after irradiation to determine a partial- or total-body exposure pattern. Also, the changes in gene expression (up to a 40-fold difference between different exposure patterns) appeared promising (Fig. 3). However, the overlap in deregulated mRNA species at multiple time points was very low (Fig. 2). This finding, along with the wave-like increase and decrease in the number of deregulated mRNA, points to the activation of several transient biological processes which in most instances return to normal during the examined time frame. Corresponding bioinformatic analysis for gene enrichment of related biological processes revealed almost no overlap within the examined time points (Supplementary Table S2; https://doi.org/10.1667/RADE20-00121.1.S2). As a consequence, the determination of exposure pattern based on mRNA gene expression changes is complicated because a different set of candidate mRNAs appears to be required at each day after irradiation. We acknowledge that gene expression changes occurring between the large time points examined are unknown. Regardless, the identified genes might be applicable over several days, e.g., 7–10 days but not at 28 days.
Low sample numbers per exposure group made significant findings by chance likely. Independent validation in another group of animals was not possible, so we performed different comparisons with the expectation of finding similar predictive capabilities of the candidate mRNAs. For instance, when identifying mRNAs to discriminate TBI from PBI, we assumed a considerable overlap in mRNAs which discriminate TBI from LHB as well, since LHB represents a subset of PBI baboons (Table 1). Indeed, approximately 30% (131/441) of the identified mRNAs discriminated both exposure pattern comparisons (Fig. 2D).
As an alternate measure to assess the validity of our results, we analyzed the association of gene expression changes corresponding to the percentage of the exposed body area since partial or total-body irradiation should affect target organs (e.g., bone marrow or the peripheral blood) proportionally with corresponding gene expression changes. This allowed simultaneous examination of all 17 samples with a concomitant increase in power (linear regression analysis instead of constructing categories and using non-parametric tests). This resulted in strong associations and low P values with most (80%) of 10 candidate mRNAs surviving conservative Bonferroni correction (Table 2, Supplementary Table S1; https://doi.org/10.1667/RADE-20-00121.1.S1).
Some limitations of our study should be noted. We performed gene expression measurements using human genomic sequences, because the baboon transcriptome was not publicly available. Given the high homology of both genomes (93%) and the previously reported work of other groups (26–28 ), we proceeded as described above. Furthermore, in previously published studies with another focus (e.g., prediction of ARS severity) using the same baboon database, we performed an independent validation of the whole genome microarray data (16 ). For that we employed another methodology (qRT-PCR) in split blood samples and showed good agreement between the microarray and qRT-PCR data [linear regression, rsq = 0.92, P < 0.0001, (16 )]. These results provide evidence for the current approach based on microarray gene expression measurements only.
Nevertheless, due to the small sample size, we judge our study more as a “proof of concept” study in which we explore the potential of easily identified mRNAs in the peripheral blood and their usefulness to determine the exposure pattern (partial vs. total-body irradiation). Knowledge of differing body areas exposed to radiation will complement information on the absorbed dose and may better predict a clinical outcome.
In summary, our data suggest mRNA changes could be used to presume the pattern of radiation exposure on the human body. Compared to our previously performed miRNA study, mRNA gene expression changes likely can be applied over a much larger time frame (2–75 days after exposure), but due to the transient gene expression changes, a different set of candidate mRNAs appears to be required at each day after irradiation. Future work will require larger sample sizes and results could be validated in additional species.
Table S1. Summary of qRT-PCR results [normalized threshold cycle (Ct) values] of candidate mRNAs per exposure group comparison (subtitles) over time after irradiation. Numbers after the subtitles in parenthesis indicate the number of candidate mRNAs per exposure group comparison. Descriptive statistics reflect the distribution and the fold change (FC) difference among mRNAs of the first exposure group (e.g., TBI) relative to the second exposure group (e.g., PBI used as the reference). P values corrected for and surviving multiple comparisons (Bonferroni) are shown in bold face. Results from logistic or linear regression are shown in the right-side columns of the table. Linear regression results for the groups TBI vs. LHB and UBE vs. LHB are omitted due to the low number of candidate mRNAs. TBI = total-body irradiation; PBI = partial-body irradiation; LHB = left hemibody exposure; UBE = upper body exposure; stdev = standard deviation; KW = Kruskal Wallis; OR = odds ratio; 95% CI = 95% confidence interval; ROC = receiver-operator characteristic; comp. sep. = complete separation; log = logistic; lin = linear; stderr = standard error of the mean.
Table S2. Summary of the bioinformatics analysis on all genes associated with a >2-fold gene expression difference (up or down) for TBI vs. PBI groups for each time point, separately. Gene set enrichment analysis was performed using PANTHER pathway software (version 15.0, http://www.pantherdb.org). The PANTHER tool groups genes with similar molecular or biological function based on their annotation (reference list was the current Homo sapiens GO database). Results were created using Fisher's exact test with a false discovery rate (FDR) multiple test correction displaying only results with FDR <0.05. TBI = total-body irradiation; PBI = partial-body irradiation; LHB = left hemibody exposure; UBE = upper body exposure; FDR = false discovery rate.
The sophisticated and carefully performed technical assistance of Sven Doucha-Senf and Eva Grumpelt are very much appreciated. This work was supported by both the French and the German Ministry of Defense.
Editor's note. The online version of this article (DOI: https://doi.org/10.1667/RADE-20-00121.1) contains supplementary information that is available to all authorized users.