Inhalation and ingestion of 137Cs and other long-lived radionuclides can occur after large-scale accidental or malicious radioactive contamination incidents, resulting in a complex temporal pattern of radiation dose/dose rate, influenced by radionuclide pharmacokinetics and chemical properties. High-throughput radiation biodosimetry techniques for such internal exposure are needed to assess potential risks of short-term toxicity and delayed effects (e.g., carcinogenesis) for exposed individuals. Previously, we used γ-H2AX to reconstruct injected 137Cs activity in experimentally-exposed mice, and converted activity values into radiation doses based on time since injection and 137Cs-elimination kinetics. In the current study, we sought to assess the feasibility and possible advantages of combining γ-H2AX with transcriptomics to improve 137Cs activity reconstructions. We selected five genes (Atf5, Hist2h2aa2, Olfr358, Psrc1, Hist2h2ac) with strong statistically-significant Spearman's correlations with injected activity and stable expression over time after 137Cs injection. The geometric mean of log-transformed signals of these five genes, combined with γ-H2AX fluorescence, were used as predictors in a nonlinear model for reconstructing injected 137Cs activity. The coefficient of determination (R2) comparing actual and reconstructed activities was 0.91 and root mean squared error (RMSE) was 0.95 MBq. These metrics remained stable when the model was fitted to a randomly-selected half of the data and tested on the other half, repeated 100 times. Model performance was significantly better when compared to our previous analysis using γ-H2AX alone, and when compared to an analysis where genes are used without γ-H2AX, suggesting that integrating γ-H2AX with gene expression provides an important advantage. Our findings show a proof of principle that integration of radiation-responsive biomarkers from different fields is promising for radiation biodosimetry of internal emitters.
There is an important need for radiation biodosimetry in response to large-scale radiological events such as improvised nuclear device (IND) detonations or nuclear power plant accidents (e.g., Chernobyl or Fukushima) (1–10). The primary tasks for biodosimetry involve secondary triage of large numbers of potentially irradiated individuals to identify those with potential radiological injuries, after initial clinical triage, and/or to help estimate radiation doses for other individuals/groups who received less critical/life-threatening doses (11, 12). To accomplish this, there is a need to quantitatively reconstruct radiation dose and to inform exposed individuals as quickly as possible about the magnitude of their exposure and potential consequences and treatment options (13).
Importantly, situations where radiation biodosimetry can prove useful are not limited to acute external exposures, but include more complex exposure scenarios like inhalation and ingestion of long-lived radionuclides such as 137Cs into the body (14). These situations of internal exposure can result in short-term toxicity (radiation sickness) and/or delayed effects such as carcinogenesis (15–17). Compared to acute external exposures, internal emitter exposures offer several challenges for biodosimetry-based dose reconstruction, because dose and dose rate from internal exposure change over time in a nonlinear manner due to radionuclide pharmacokinetics (e.g., distribution from blood to various organs, elimination from the body). Dose and dose rate are not independent variables because both depend on the initial administered activity and on time (18–20).
In collaboration with Lovelace Biomedical Research Institute (LBRI; Albuquerque, NM), we developed a mouse model system where various amounts of 137CsCl were injected into the animals, and excretion kinetics, dose rates, doses, γ-H2AX fluorescence (20, 21), and transcriptomic (22, 23), metabolomic (24) and lipidomic (25) biomarkers were measured over time after injection. More recently, we used a γ-H2AX end point in peripheral blood mononuclear cells to reconstruct the 137Cs activity after injection with several different 137Cs amounts (20). Activity values were converted into radiation doses, using the time since injection and 137Cs pharmacokinetics information. The model-based reconstructions were most accurate between 2–5 days of exposure due to persistent expression of γ-H2AX at these times after 137Cs injection. This work highlighted the feasibility of using γ-H2AX as a biomarker to assess injected 137Cs activity. We proposed that combining γ-H2AX with other biomarkers may have the potential to more accurately reconstruct 137Cs activity and the resulting absorbed dose and biological damage in the blood over a longer period of time.
The additional biomarker we chose was transcriptomics response. Results of the global gene expression profiling analysis are reported in Ghandhi et al. (23). The goal of this work was to identify genes that could serve as robust and reliable biomarkers of protracted internal radiation exposure. Such genes should have a strong radiation response that is consistent across several radiation studies and is stable over days–weeks after the start of exposure. We also assessed the feasibility and potential advantages of combining such gene expression data with γ-H2AX for reconstructing 137Cs activity.
MATERIALS AND METHODS
Male C57BL/6 mice were injected with varying amounts of 137CsCl, and sacrificed at times from 2 to 14 days later for collection of blood and subsequent measurement of biodosimetry end points, as described elsewhere (20, 23). Male mice were used to build on the results of these previous studies, and the data were used as an initial test of gene expression and cytogenetic data integration. All animal studies were conducted at Lovelace Biomedical Research Institute in facilities accredited by the Association for Assessment and Accreditation of Laboratory Animal Care (AAALAC) using protocols approved by the Lovelace Institutional Animal Care and Use Committee (IACUC; approved protocol no. FY15–087).
Gamma-H2AX Fluorescence Measurements
Immunodetection of γ-H2AX was performed according to our previous work reported elsewhere (20, 21, 26). Briefly, isolated blood mononuclear cells (MNC) were blocked with 3% bovine serum albumin (BSA; Sigma-Aldrich® LLC, St. Louis, MO) for 30 min at room temperature and incubated with a rabbit polyclonal γ-H2AX (phospho S139) antibody (dilution 1:600; Abcam®, Cambridge, MA), visualized with a goat anti-rabbit Alexa Fluor® 488 antibody (dilution 1:1,000; Invitrogen™, Carlsbad, CA) and counterstained with DAPI in Vectashield® mounting media (Vector® Laboratories, Burlingame, CA). Fluorescent images of cells were captured using an Olympus epifluorescence microscope (Olympus BH2-RFCA; 60× oil immersion objective) and quantification of γ-H2AX yields was determined by measuring the total γ-H2AX nuclear fluorescence per MNC. Data were analyzed using image analysis software described elsewhere (26).
Gene Expression Measurements
A brief description is included here of the gene expression results that were the basis for the gene selection for the activity reconstruction algorithm developed in this study. RNA samples from six animals per time point per injection group were used to generate transcriptomic data using Whole Genome Mouse Arrays (Agilent Technologies Inc., Santa Clara, CA), using the recommended manufacturer's protocol. The gene expression data set is deposited in the NCBI Gene Expression Omnibus database (accession no. GSE118616) and is described in detail in Ghandhi et al. (23).
Selection of Genes for 137Cs Activity Reconstruction
The research strategy was to focus on those genes that have the largest (absolute value) Spearman's correlation coefficients with injected 137Cs activity and near-zero correlation coefficients with time after injection. The search was performed using R 3.5.2 software (27). We selected the top five genes (Atf5, Hist2h2aa2, Olfr358, Psrc1 and Hist2h2ac), ranked in decreasing order by Spearman's correlation coefficients (all statistically significant at α = 0.05 level with Bonferroni correction) with injected activity, for further analysis. Sensitivity calculations using all statistically significant candidate genes (there were eight, shown in the Supplementary_data_genes.scv file in the Supplementary Materials (491_rare-195-01-03_s01.xlsx); https://doi.org/10.1667/RADE-20-00042.1.S1) instead of the top five, and comparing the performances of these model variants by Akaike information criterion (AIC), coefficient of determination (R2) and root mean squared error (RMSE), showed very similar results, so we proceeded with the top five.
We also assessed the dose-response patterns of these five genes in other published studies: 1. Our earlier internal 137Cs study using a single injection activity (22) (Dataset: GEO accession no. GSE52690); and 2. Several published studies of acute external irradiation (28–31) (Datasets: GEO accession nos. GSE124612, GSE114142, GSE99176). Consistency of dose-response patterns across all/most of these studies would provide improved confidence that the selected genes represent reliable candidates for radiation biodosimetry in a variety of exposure scenarios.
Combining γ-H2AX and Gene Expression Data Sets
Blood samples from eight mice per exposure group were pooled for γ-H2AX analysis, and six mice per group for gene expression analysis, for technical reasons: Due to reduced blood mononuclear cell counts after radiation exposure in a relatively small volume of blood (∼100–200 µl) from each mouse, it was necessary for the γ-H2AX assay to pool the blood from several mice to obtain enough cells for robust γ-H2AX analysis. Gene expression data were originally measured for individual mice, with the results being pooled to allow integration with γ-H2AX data.
Consequently, mean injected activity values for the mouse groups differed slightly for the γ-H2AX and gene expression data sets. To combine the data sets for our biomarker integration analysis, the injected activity values were averaged from the γ-H2AX and gene expression studies. The numerical changes from this procedure were very small, since activity values for the exposure groups were modified by only 0.003 to 0.062 MBq, whereas the range of studied non-zero activities was 5.76 to 9.29 MBq.
For the injected activity, γ-H2AX and gene expression measurements, the uncertainties were fairly small on a relative scale, usually a few percentages of the respective means. For injected activity, the range of variations between mouse groups was –3.2% to +4.4% of the mean values. For γ-H2AX fluorescence, standard errors were on average 2.3% of the mean (range: 1.2 to 3.0%). For gene expression, standard errors were on average 11.1% of the mean (range: 3.5 to 23.2%).
Because the signal values (log2-transformed and normalized) of these top five selected genes were strongly correlated with each other across the samples, each gene was not treated as an independent predictor. The geometric mean of signal values for all five genes (designated as the variable Net_sig) was used as a single robust predictor of injection activity. We assessed the reasonableness of this assumption by removing each of the five genes one by one and performing the analysis using the geometric mean of the remaining four genes. Only removal of Psrc1 significantly reduced R2 for the resulting model, whereas removal of any of the other genes did not affect R2 significantly.
Gamma-H2AX mean fluorescence (designated as the variable gH2AX, in arbitrary units) was used as the second predictor. Quadratic Net_sig2 and gH2AX2 terms were also considered to assess potential non-linear relationships. Time after radioactivity injection was the final tested predictor variable. All time points (2–14 days after 137Cs injection) were used in the analysis.
Mathematical Model Development and Fitting
Multimodel inference (MMI) (32, 33) using the Akaike information criterion with sample size correction (AICc) was performed using the glmulti R package on multiple linear regressions with all possible combinations of these predictor variables (main effects and interactions). Only those predictors that had the highest MMI importance scores were retained for further analysis.
We used the retained predictor variables to construct a three-parameter nonlinear model. The model structure was not intended to represent any specific radiation response mechanisms, which are not yet well known for these genes. It was simply intended to capture the main data patterns and to generate only positive numbers for reconstructed injected activity (unlike linear or polynomial regressions, which can generate negative predictions). The model's performance was assessed by R2 and RMSE.
The robustness of these metrics was assessed by calculating them not only on the full data set, but also on 100 random 50:50 splits of the data into training and testing halves. Best-fit model parameters obtained on the training data were used to generate model predictions and corresponding R2 and RMSE values on the testing data.
We quantified the sensitivity of model predictions to ±10% changes in each predictor variable (gH2AX or Net_sig) by calculating how these changes affect R2 and RMSE metrics at each time point (2, 3, 5, 7 or 14 days after 137Cs injection), or over all time points combined. Model parameter values were kept constant at best-fit values. Potential time trends for each metric were also assessed using Pearson correlation coefficients with time after injection.
The top five genes that had large and statistically significant Spearman's correlation coefficients with injected activity and remained stable over time after injection were Atf5, Hist2h2aa2, Olfr358, Psrc1 and Hist2h2ac (Table 1). All the selected genes had similarly-shaped dependences on injected activity that were strongly correlated with each other (Figs. 1–2). Combining them into a single variable by using the geometric mean of the log-transformed gene signals (Net_sig) provided a good biomarker for injected activity (Figs. 2–3, Supplementary Fig. S1 (491_rare-195-01-03_s02.pdf); https://doi.org/10.1667/RADE-20-00042.1.S2).
Five Top-Scoring Genes Selected for the Current Biodosimetry Analysis
As an initial validation, the same polynomial regression model relating Net_sig to injected activity (Fig. 3) was applied to an independent data set from our initial 137Cs injection study (22). On this data set, the model-based reconstructed activity was 6.51 MBq, whereas the actual value was 7.67 MBq. This error of 1.16 MBq (∼15% on a relative scale) is similar to the model's RMSE on the original data set (0.99 MBq). Therefore, the five selected genes performed similarly for activity reconstruction across two independent data sets. The selected genes also tended to increase in expression with increasing radiation dose in studies of acute external irradiation (28, 31, 34). These findings further support these genes as reliable candidates for internal emitter biodosimetry.
Integration of Gene Expression with γ-H2AX
Gamma-H2AX fluorescence showed curving dependences for injected activity and time (Fig. 4), which likely reflect the kinetics of the biomarker itself as well as changes in the MNC population due to the death of heavily-damaged cells and proliferation of new ones (20). To determine the best way to combine γ-H2AX with gene expression, we evaluated all possible combinations of time since injection, gH2AX, Net_sig, and their quadratic terms and interactions to reconstruct injected activity by multiple linear regression. An importance score based on Akaike information criterion weights was calculated for each predictor variable and interaction term (32, 33).
Time and its interaction terms with other variables had the lowest importance scores, <0.012. The analysis of multiple regression combinations was then repeated without time as a variable. This showed that the quadratic Net_sig2 term and its interactions with other variables now had the lowest importance scores, <0.28. The highest scores of 0.73 and 0.65, respectively, were achieved by two interactions: Net_sig×gH2AX2 and Net_sig×gH2AX. The main effects of these variables had lower scores, <0.56.
Based on these findings, we developed a nonlinear model using the interactions Net_sig×gH2AX2 and Net_sig×gH2AX as predictors of injected activity. The model contained three parameters (k1, k2 and k3) and had the following structure, where gH2AX represents the γ-H2AX fluorescence signal and Net_sig represents the combination of five selected genes:
As described in Materials and Methods, the model is empirical rather than mechanistically motivated. The goal was to create a parsimonious and easily-applied formalism for activity reconstruction. The model was fitted to the data using robust nonlinear regression (nlrob function in R), which reduces the influence of outlier data.
Model-Based Reconstructions of Injected 137Cs Activity
Best-fit model parameters and performance metrics (R2 and RMSE) are presented in Table 2, and model behaviors are plotted in Fig. 5. The model was generally quite accurate in reconstructing the injected activity values (Fig. 6). The R2 and RMSE metrics remained relatively stable when the model was fitted to one half of the data, which had been randomly selected, and tested on the other half, repeated 100 times (Table 2). Such stability suggests that the model is not strongly over-fitting and is robust to random variations of the data.
Robustness Testing Results for the Nonlinear Model that Used γ-H2AX Fluorescence and Geometric Mean of Selected Gene Expression Levels (Net_sig) to Reconstruct Injected Activity [Eq. (1)]
Since the γ-H2AX radiation response was stronger at short times (≤5 days) after injection than at longer times (20), we calculated RMSE for the current integrated model (Eq. 1) for times ≤5 days on testing data sets. The mean RMSE values were somewhat lower for times ≤5 days than for all times: 1.058 (SD = 0.360) vs. 1.109 (0.266), respectively. However, these differences were small, suggesting that the integrated model had similarly good accuracy across the tested time range (0–14 days after injection). Such similarity of RMSE across time points resulted from the contribution of genes (variable Net_sig) because the genes selected for this analysis had stable expression patterns over time.
More detailed examination of model performance metrics R2 and RMSE over time after 137Cs injection showed that Pearson correlation coefficients with time were 0.34 (P = 0.58) and –0.35 (P = 0.56) for R2 and RMSE, respectively. Therefore, both metrics varied somewhat with time, but their variations did not reach statistical significance. These results suggest that the agreement between actual injected activity and model predictions was consistent over all the analyzed time points (up to 14 days after injection).
Sensitivity analysis of model predictions using ±10% changes in each predictor variable (gH2AX or Net_sig) showed that the sensitivity was similar for each predictor. R2 values were reduced by 0.22 to 0.33, and RMSE values were increased by 0.99 to 1.37 MBq by 10% perturbations of the predictor variables. These R2 and RMSE changes were not significantly correlated with time for any of the sensitivity analysis combinations. The Pearson correlation coefficient with time approached statistical significance (–0.83, P=0.081) only for the RMSE metric, when gH2AX was increased by 10%. This finding suggests that increasing gH2AX causes RMSE to decrease with time, probably because gH2AX is a more reliable predictor of injected activity at early times after injection (20), but its performance decreases at later times.
The mean absolute errors in activity reconstruction based on the current model [Eq. (1)] over all time points were 2.2-fold smaller than in our previously published analysis using γ-H2AX alone (20). This result suggests that integrating γ-H2AX with gene expression provides a significant advantage relative to using γ-H2AX alone (P = 0.024 by paired t test comparing absolute error magnitudes). Alternatively, to test whether or not inclusion of gH2AX as a predictor improves model performance compared to using genes only, we produced and fitted a model variant which contained only Net_sig by eliminating the gH2AX variable from Eq. (1) (replacing it with 1). This genes-only model variant had R2 = 0.904 and RMSE = 0.996 MBq, which are worse than these metrics for the original model with both gH2AX and Net_sig (Table 2). Comparison of these model variants using Akaike information criterion (AIC) showed that the genes-only model performed worse than the original model by 8.2 AIC units, which translates into an evidence ratio of exp[–8.2/2] ∼ 0.016. Therefore, using genes only instead of combining genes with gH2AX resulted in significant (evidence ratio <0.05) loss of model performance.
Radiation biodosimetry is needed for an effective response to large-scale accidental or malicious events involving ionizing radiation (1–10). The importance of high-throughput biodosimetry techniques and assays after a large-scale radiological/nuclear event is well recognized (3, 35–39). This provides a rationale for using assays such as proteomic, transcriptomic, metabolomic and cytogenetic end points, which can be automated to provide high-throughput biodosimetry approaches (4, 40–46). While standard cytogenetic biodosimetry techniques such as the dicentric chromosome (DCA) and micronuclei (CBMN) assays are known for their accuracy, γ-H2AX assay and gene expression assays using qRT-PCR potentially have the advantage of shorter time-to-result and time-to-dose reconstruction because they do not require cell culturing (47).
Combining multiple biodosimetric markers is generally thought to improve the performance of dose reconstruction algorithms (48). Early work in this area showed that combinations of blood cell counts and serum levels of one or more proteins improved discrimination of individual mice or non-human primates with severe radiation exposure from those without exposure (49–51). A similar approach using a small panel of proteins has also been shown to improve dose assessment and prediction of the severity of injury in multiple animal models (52, 53). Software platforms, such as MULTIBIODOSE (54) and the Biodosimetry Assessment Tool (55), have been developed to easily combine the results of various individual assays to produce a single consensus triage dose for use in accidents. Such a multiparametric approach, combining complete blood counts with plasma levels of Flt3, has actually been tested on 63 individuals involved in a radiation accident in Senegal, and was found to give results similar to those from scoring chromosome rings and dicentrics (56).
Here we performed a proof-of-principle study showing that integration of radiation-responsive biomarkers using two different assay approaches, gene expression and H2AX protein phosphorylation levels, is promising for radiation biodosimetry of internal emitters. The results suggest that γ-H2AX as well as specific gene signals are strongly correlated with the amount of 137Cs in the body and remain sufficiently stable over several weeks after 137Cs administration. Reconstructed injected activity values can be converted into radiation doses, using the time since injection and known 137Cs pharmacokinetics information, because dose rate is a function of the 137Cs activity remaining in the body at a given time and dose is the time integral of the dose rate (20).
The computational strategy selected here to perform the 137Cs activity reconstructions focused on identifying a small number of candidate genes with robust radiation responses and stability over time. The inclusion of multiple genes is thought to strengthen biomarker performance against potential confounding factors or individual variation (57).
Most of the genes selected for inclusion in the activity reconstruction biomarker have documented radiation responses in other studies. Two of the genes, Atf5 (activating transcription factor 5) and Psrc1 (proline/serine-rich coiled-coil 1) are linked to the Tp53 signaling network, a central part of the radiation response. Atf5 is a transcription factor that itself appears to increase radiation survival by abrogating the function of Tp53, specifically in blocking radiation-induced p53-dependent apoptosis, as shown in mouse sarcoma cells (58). Psrc1 is a direct transcriptional target of Tp53, and is involved in microtubule polymerization (59). Hist2h2aa2 and Hist2h2ac encode core nucleosome proteins, involved in wrapping DNA into chromatin, limiting DNA damage and increasing chromatin stability. Many histone genes have been shown to be downregulated in response to DNA damage both via the p53-dependent G1 checkpoint response (60) and through a p53-independent mechanism (61). Transcriptomic studies have also shown increased expression of histone 2h2a family genes after radiation exposure in mice (28), non-human primates (62) and humans (63).
Model Performance Metrics at Different Time Points after 137Cs Injection, and for All Time Points Combined
The Olfr358 gene codes for a protein that is predicted to have olfactory receptor activity via G-protein coupled receptor signaling. Genes coding for proteins in this signaling pathway, including numerous members of the Olfr family, are often found to be strongly responsive in mouse irradiation studies (22, 28). Thus, this likely represents a general signaling response related to stress rather than olfaction. The highly correlated pattern of response for the selected genes suggests they may have a common regulator or signaling pathway, perhaps triggered by sustained DNA damage and Tp53-driven transcriptional activation.
To reduce the number of predictors and adjustable parameters, and to alleviate potential problems with multicollinearity of predictors, the five selected genes were not treated as independent predictors, but combined into one predictor variable (Net_sig) using the geometric mean of their signal values. The choice of model structure for integrating the gene-based predictor Net_sig with γ-H2AX was assisted by information theoretic MMI analysis.
This approach resulted in a simple mathematical model with three adjustable parameters [Eq. (1)]. Because of its simplicity, the model was quite stable to variations in the data such as random splits of the data set into training and testing halves. The limited numbers of model parameters and predictors are reasonable to avoid overfitting, because the data set was small and inter-individual variations were “smoothed” due to the technical necessity of pooling the data from multiple samples, so that each analyzed data point represents an average for a group of mice rather than the data from each individual mouse.
While injection of 137Cs into mice provides a useful model for the study of radioactive cesium as an internal emitter, it does present experimental and logistical challenges. Because the mice and all associated biomaterials are radioactively contaminated, they require special handling and costly clean-up and decontamination procedures. The mice must be housed individually to avoid cross-irradiation between animals, and contaminated biospecimens must be analyzed using dedicated equipment (64, 65). The recent development of our innovative variable low-dose-rate gamma irradiator (VAriable Dose-rate External irradiatoR; VADER) that can provide external 137Cs irradiations at dose rates below 1 Gy/day (65) will allow us to further validate these biomarkers by simulating realistic environmental exposure scenarios. Future work should also be extended to compare the biomarker response of these protracted exposures in female mice.
In summary, this study provided a proof of principle for the feasibility of the following points: 1. Identification of several genes with robust and consistent radiation responses, suitable for biodosimetry of internal emitters; 2. Integration of these gene data with γ-H2AX fluorescence in a simple model framework, with validation of the predictions in an independent study. Based on the results, we conclude that integrating γ-H2AX with gene expression is a promising strategy for radiation biodosimetry of internal emitters and provides a significant advantage relative to using γ-H2AX alone or genes alone, extending assay utility out to at least 14 days after the onset of exposure.
Supplementary materials (491_rare-195-01-03_s01.xlsx). Supplementary_data_genes.csv. This file contains the data for all eight candidate genes and their Spearman's correlation coefficients and P values with injected activity.
Fig. S1 (491_rare-195-01-03_s02.pdf). The relationships between injected activity and gene signal intensity (log2-transformed mean intensities for six mice per sample) for each of the five selected genes with the strongest positive Spearman's correlations with injected activity (Table 1), and for the geometric mean of these gene levels.
We thank Ms. Aimee Kowell, and the technical staff at the Lovelace Biomedical and Environmental Research Institute for conducting the 137Cs dosing experiment that generated the data used in this study. We also gratefully acknowledge Mr. Shad R. Morton for assistance with the gene expression data and preparation of figures. This work was supported by the National Institute of Allergy and Infectious Diseases (NIAID grant no. U19-AI067773 to the Center for High-Throughput Minimally Invasive Radiation Biodosimetry).