Dose and Dose-Rate Effects in a Mouse Model of Internal Exposure to ¹³⁷Cs. Part 1: Global Transcriptomic Responses in Blood

The potential health risks to humans after radiological accidents such as the Fukushima and Chernobyl nuclear disasters (1–3 ) have led to concerns about radiation exposure levels and potential treatments (4 ). In addition, in the event of a possible detonation of a radiological dispersal device (RDD) or explosion of an improvised nuclear device (IND) in a highly populated area (5 ), there would be widespread health concerns and a need for rapid biodosimetry to measure exposure levels in potential victims. All of these scenarios may produce fallout containing radioactive isotopes, including ¹³⁷Cs, ⁹⁰Sr and ¹³¹I. The internal exposures resulting from inhalation or ingestion of these isotopes may complicate assessment of radiological injuries. Biomarkers of radiation exposure have been developed that accurately assess acute exposures to radiation and help stratify at-risk individuals who may benefit from immediate treatment to mitigate early detrimental effects on their health (6–10 ). Internal contamination from radioisotopes in the environment may be incorporated into tissues and organs, however, with the resulting exposure continuing for extended periods from weeks to years (1 ). The cumulative doses from protracted exposures can be comparable to acute doses, but because of longer time periods after the initial event, dose-rate effects may significantly impact the health and treatment needs of such individuals.

Several groups, including ours, have focused on the development of gene expression signatures as biomarkers of acute radiation exposures and dose (11–20 ). Much less is known, however, about the effect of different dose rates on such gene expression signatures (12, 14, 21 ), or how this may complicate evaluation of exposure when internal emitters are involved. As an initial study of the in vivo biological response to internally deposited ¹³⁷Cs, we previously gave mice a single ∼8 MBq injection of cesium chloride (¹³⁷CsCl) solution (5, 13, 22, 23 ) and reported a large number of differentially expressed genes in a 30-day study. In general, expression patterns were more complex than those seen after acute exposures, although some genes were consistently underexpressed throughout the study. Another pattern of expression included many known radiation-response and TP53-regulated genes, which showed increased expression at early times followed by a drop to near or below control levels at later times despite the higher accrued dose (13 ). Blood and biofluids from the same animals were also analyzed for cytogenetic (23 ) and metabolomic (24 ) changes.

We have now extended these studies to include four injection dose groups followed over the course of two weeks, to broaden our understanding of the dose and dose-rate effects on the blood cell transcriptomic response to ¹³⁷Cs in vivo. The mice were injected with 0, 5.71, 6.78, 7.67 or 9.29 MBq of ¹³⁷CsCl, which resulted in maximum doses of 0, 4.3, 5.6, 7.1 and 12 Gy at the two-week point, respectively. The study was designed to achieve an accrued dose of ∼4 Gy at days 3, 5, 7 and 14 in the time course, for comparison with the initial study. Some animals in the study received similar doses but at different dose rates (and times). Other animals received different total doses with similar dose rates (at different times).

Mouse experiments were performed under IACUC approved protocol no. FY15-087 following the guidelines of The Lovelace Biomedical Research Institute (LBRI; Albuquerque, NM), which is fully accredited by the Association for Assessment and Accreditation of Laboratory Animal Care (AAALAC). This study was in compliance with all applicable sections of the Final Rules of the Animal Welfare Act regulations (9 CFR parts 1 and 3) and the Guide for the Care and Use of Laboratory Animals, National Research Council (25 ). Male C57BL/6 mice (∼8–10 weeks old, ∼20–30 g in weight, from Charles River Laboratories, Frederick, MD) were initially quarantined (at least 14 days) at LBRI, and then randomized by body weight into five injection groups of 40 animals per group. The mice were given a single intraperitoneal (IP) injection to deliver target levels of 0, 5.74, 6.66, 7.65 or 9.28 MBq of ¹³⁷CsCl solution. To minimize cross irradiations between animals, housing cages were shielded with lead. Mice were euthanized at days 2, 3, 5, 7 and 14 after injection, with the experiment being split into three independent runs to achieve a final number of 8 mice for each combination of injection amount and sacrifice time [for further details of mouse handling and treatment, see Turner et al. and supplementary files therein (22 )]. Animals were examined twice per day throughout the study. Onset or progression of abnormal clinical signs, and identification of dead, weak or moribund animals were carefully recorded. Special attention was paid to dehydration, hunched posture, severely ruffled coat, reluctance to move, lateral recumbency and excessive weight loss (greater than 20% compared to randomization), none of which were observed in this study. No differences were found between the three independent runs for any measured physiological or morphological parameters, therefore, results are presented for all animals in an injection group as the average for the group. The vehicle group received the same pH-adjusted saline solution that was used as the vehicle for the ¹³⁷CsCl injections. The study was designed to target a cumulative 4 Gy dose at day 14 (5.74 MBq injection amount), day 7 (6.66 MBq injection amount), day 5 (7.65 MBq injection amount) and day 3 (9.28 MBq injection amount) after delivery of the radioisotope. This study follows up on previously published work from our laboratory, in which the cumulative dose at day 5 after injection was 4 Gy (13 ).

After irradiation, animals were analyzed on a whole-body counting system daily for the first seven days, then again on days 10 and 14. Animals were placed in a small container with air holes for approximately 5–10 min and placed on the in vivo counter to measure radioactivity for 3 min per animal. Vehicle-injected animals were handled similarly, but without counting, to serve as controls for the ¹³⁷CsCl-injected animals.

A detailed description of the calculations of dose rate is provided elsewhere (22 ), in which blood samples were shared with this study. Briefly, the parameters of the whole-body retention equation were calculated for each of the animals in the study [Supplementary Fig. S1 (https://doi.org/10.1667/RADE-20-00041.1.S1) shows calculated mean daily dose rates]. For animals allocated to shorter time points (up to 3 days), just one exponential term was observed (Eq. 1). The individual parameters of the retention equation were used to calculate the total number of disintegrations of ¹³⁷Cs in the body. To calculate the dose rate for each day after injection, the total number of disintegrations was calculated for each time period after ¹³⁷Cs injection as follows:

The dose coefficient (Gy.Bq^–1 of administered activity) was calculated for each day after ¹³⁷Cs administration, multiplying the total number of disintegrations (Bq.s^–1) by the specific effective energy (SEE) value for ¹³⁷Cs (2.41 × 10^–19 Gy.dis^–1).

The dose coefficient is calculated using Eq. (2):

where: Ã(S) are the total number of disintegrations that will occur within source organs/tissues over time, SEE (T←S) is the specific effective energy in MeV kg^–1, i.e., the energy imparted per kg of a target tissue (T) as a consequence of a disintegration of ¹³⁷Cs in source organs/tissues (S) = 2.41 × 10^–19Gy.dis^–1. The absorbed doses for each day after ¹³⁷Cs administration were calculated by multiplying the administered activity by the dose coefficient for that given day (Eq. 3). The dose rate for a given day was calculated by subtracting the committed dose calculated for that given day from that of the previous day's.

where D = absorbed dose (Gy) and A = administered activity (Bq).

Animals were euthanized by IP injection of Euthasol (diluted with sterile saline solution) following the guidelines in the LBRI SOP ACS-0334 (Euthanasia of Small Animals). Blood samples were collected immediately after euthanasia by cardiac puncture in a sterile hood. Blood (400–500 µl) was added to PAXgene® RNA solution (Becton, Dickinson and Co., Franklin Lakes, NJ) at a 1:5 ratio and inverted ≥10 times to mix. Samples were stored at –80°C for the length of the study and then shipped together at the end of the study period to the Center for Radiological Research at Columbia University Irving Medical Center (New York, NY). RNA samples were prepared using the PAXgene protocol, globin-depleted using the mouse Ambion® GLOBINclear kit (Ambion®; Thermo Fisher Scientific™ Inc., Waltham, MA), and quantified using a NanoDrop™ One spectrophotometer (Thermo Fisher Scientific). RNA quality was determined using the Agilent 2100 Bioanalyzer (Santa Clara, CA), and RIN values ranged from 8.0 to 9.5 across all hybridized samples.

Globin-cleared RNA samples were hybridized to Agilent Whole Genome Mouse 4X44 v2 (G4846A) microarrays using the manufacturer's protocol. Feature extraction of the images was performed using Agilent's FE version 10.7 software and data were imported into BRB-ArrayTools software version 4.6.0 (26 ). The gene expression dataset was deposited in the NCBI GEO database (accession no. GSE118616). We identified significantly differentially expressed genes (using the cutoffs of P < 0.001 and FDR < 5%) in all subgroups in the study by comparing them with the time-matched controls (Supplementary Table S1; https://doi.org/10.1667/RADE-20-00041.S2). Gene ontology analysis was performed using the PANTHER database with the mouse annotations database (27 ) and REACTOME pathways and biological processes categories. Multiple testing was performed based on Bonferroni correction (28 ) and adjusted P values are reported. We also analyzed sets of up- and downregulated genes separately, and looked for common biological processes among these genes. For smaller gene lists (∼100 genes), up- and downregulated genes were not separated for analysis of enrichment of biological processes.

We applied maSigPro (29 ) to the 15,008 genes that passed the BRB-ArrayTools filters (greater than 25% values ±1.5-fold from the median and less than 25% missing values for a gene) and identified 10,579 genes with significant temporal profiles. Next, we used forward regression with a false discovery rate cutoff (a = 0.01 and R² = 0.04) to group these genes into five clusters that showed distinct patterns of gene expression in the four injection groups across the time course.

This study was designed to accumulate a dose of 4 Gy from internal ¹³⁷Cs at days 3, 5, 7 and 14 in the different injection groups, with the goal of measuring biological responses at the cytogenetic and molecular levels of mRNA and metabolites. Across the study there were points that achieved similar doses at different dose rates and conversely, similar dose rates to achieve different doses. We used the daily dose rates on the day of sacrifice as the independent variable in our analyses. mRNA levels measured using microarrays from isolated total RNA are a snapshot in time and represent a balance between existing stable mRNA levels and degradation. Recently published studies have shown that mRNA level changes, in response to stresses such as dithiothreitol treatment, may occur in a pulse-like pattern going back to near-background levels at approximately 30 h after treatment, compared to the corresponding protein level changes that are more stable over time (30, 31 ). Therefore, we focused our statistical comparisons between the dependent variable, gene expression changes, and the independent variable, dose rate on the day prior to sacrifice, as opposed to the average dose rate from the beginning of the study. The difference between these two methods of calculating dose rates was an average of 0.18 Gy/day (range 0.05 to 0.43 Gy/day) for the animals in this study.

In addition to the maSigPro analysis, we also looked for genes that showed specific patterns of response: 1. genes that were dose-rate independent, i.e., genes that showed similar fold changes at the same dose (3, 4 or 5 Gy) regardless of the initial injection activity; and 2. genes that were dose-rate dependent, i.e. positively or negatively correlated with dose rate at the same 4 Gy dose. To identify genes with these response patterns, fold changes (compared with time-matched controls) and P values were calculated for the 15,008 filtered genes. For each gene with significant fold change at one or more dose levels, fold changes were fitted to a horizontal line (to identify dose-rate-independent genes) or to a dose-rate curve (dose-rate-dependent genes). Fit scores to each of the two patterns were calculated as the ratio of the mean absolute fitted residuals and the fitted value. Genes with a fit score less than 0.05 were considered a good fit and were selected as showing either a dose-rate-independent or dose-rate-dependent pattern of expression.

We used Taqman® assays (Life Technologies/ThermoFisher Scientific) to quantify expression of selected dose-rate-independent and dose-rate-dependent genes (Fcrla: Mm00520163_m1; H2-Eb1: Mm00439221_m1; Phlda3: Mm00449846_m1; Romo1: Mm01246687_g1). We prepared complimentary DNA (cDNA) from total mRNA using the High-Capacity® cDNA Kit (Life Technologies, Foster City, CA). Quantitative real-time RT-PCR (qRT-PCR) was performed using the ABI 7900 Real Time PCR System following the manufacturer's protocol, and relative fold induction was calculated by the –ΔΔCT method (32 ) using SDS version 2.4 (Thermo Fisher Scientific). Data were normalized to the geometric mean of two housekeeping genes, Actb (Mm00607939_s1) and β2-M (Mm00437762_m1), which were found to be stably expressed in mouse mRNA using geNorm (33 ).

Gene-expression immune-cell deconvolution was performed on the microarray data using the CIBERSORT algorithm and the ImmuCC immune cell gene expression signature (34, 35 ). When genes were represented multiple times in the input data, the average gene expression was used as input to CIBERSORT. Then, the gene expression matrix containing all samples was quantile normalized and unlogged (log2). Genes not shared between the microarray assay and the ImmuCC signature were removed, resulting in a set of signal values for 511 genes that was used as input for the program that compares the gene expression to a standard comprised of 25 mouse blood cell types. Output of the reconstructed cellular composition of the whole blood used in this study was represented as relative abundance of the 25 cell types.

All animals were weighed before the start of the study and daily during the time course. Animals were healthy throughout the study and showed no adverse physical, morphological or physiological changes in the two-week time course. Retention kinetics of ¹³⁷Cs in the animals were measured for all injection groups and all groups showed similar decay kinetics (Fig. 1A), with a rapid drop to 80% of the original activity within the first two days after injection and down to 20% at day 14, which was the end of the study. These results were consistent with retention rates observed in our previously published 30-day study with one injection group (13, 23, 24 ) and are summarized in further detail in Turner et al. (22 ), which describes the cytogenetic changes in shared blood samples from the current 14-day study. The dose rates ranged from 1.4 Gy/day to 0.2 Gy/day (Table 1); and cumulative doses in this study ranged from 1.2 to 12 Gy (Fig. 1B and Table 2).

Analysis of microarray gene expression data was performed using BRB-ArrayTools software to determine the genes that are significantly differentially expressed. These results are summarized in Table 2 and detailed in Supplementary Table S1 (https://doi.org/10.1667/RADE-20-00041.S2). Genes with significantly changed expression after treatment were identified using a cutoff false positive error rate (FDR < 0.05), calculated using the method of Benjamini-Hochberg (26, 36 ). Overall, thousands of genes were found to be regulated in the four ¹³⁷CsCl injection groups (Table 2). At the earlier times (days 2 and 3 after injection) there were fewer differentially expressed genes, with numbers increasing at days 7 to 14 after injection. This appeared to be a common trend in all injection groups. This general trend towards increasing numbers of differentially regulated genes with increasing time of exposure could indicate primary stress response signaling at early times (highest dose rates) leading into secondary signaling complicated by cell turnover and continuing exposure at later times. There were proportionally more downregulated genes at later time points, suggesting that certain cell subtypes may be more sensitive to radiation and are removed from circulation, contributing to a broad loss of mRNA species. Comparison of differentially regulated genes in each injection group revealed few genes that were significantly regulated at all time points, 6 genes in the 5.71 MBq group, 85 genes in the 6.78 MBq group, 15 genes in the 7.67 MBq group and 20 genes in the 9.29 MBq group. The common genes in each injection group are listed in the Supplementary Table S1 (https://doi.org/10.1667/RADE-20-00041.S1), Venn diagram by group.

Identification of biological processes and REACTOME pathways (37 ) enriched among the differentially regulated genes was performed using the PANTHER database and gene expression tools (Fig. 2 and Supplementary Table S2 (https://doi.org/10.1667/RADE-20-00041.S3): PANTHER Gene Ontology biological processes and REACTOME table). Top significant REACTOME pathways fell into broad categories of immune system signaling and gene expression that were significant at most doses and dose rates (Fig. 2A). Other pathways that were involved later in the time course across all groups were related to olfactory signaling, mRNA splicing, GPCR signaling, chromatin modification, ribosomal transcription, post-translational modification, and epigenetic regulation of gene expression (Fig. 2B). Transcriptional regulation by the TP53 pathway also showed enrichment at later times, especially in the two highest dose groups. Cell cycle and DNA repair pathways also showed this pattern of response. REACTOME pathways that were associated with early time points in the study involved hemostasis and platelet aggregation and activation pathways (Fig. 2C). We also performed REACTOME pathway analysis on up- and downregulated genes separately to identify processes that may be consistently responsive to radiation, either up- or downregulated across all doses, and we found that the early pathways were also enriched among the upregulated genes at early times (data not shown). GLUT4 translocation to the plasma membrane, vesicle-mediated transport and Rho GTPase signaling pathways also followed the trend of regulation at early times and not beyond five days after injection in upregulated genes.

The Ingenuity Pathway Analysis (IPA) predictive analysis tool was then used to identify potential upstream regulators of the observed gene expression (38 ). In Fig. 3, we show the transcriptional regulators that were predicted to be involved in transcriptional changes at more than three of the time points in the study. MYC, IRF7 and IRF3 (predicted state inhibited, z-score <–2) and KLF3, TRIM24, NKX2-3 (predicted state activated, z-score >2) were the top transcriptional factors. Unsupervised hierarchical clustering agglomerative analysis of these regulators revealed that the predicted state clustered mostly by time within different injection groups, such as MYC inhibition occurring at day 7 across injection groups, and TP53 activation occurring early at days 2 and 3 in the time series in all groups.

Next, we used maSigPro (29 ) to identify and cluster significant patterns of gene expression change as a function of time and injected activity amount. Figure 4 shows the first three cluster patterns, comprising 85% of the clustered genes, and also the heat map of the individual genes in these clusters as a function of time and injection amounts. IPA analysis and prediction of upstream transcriptional regulators for these clusters (IPA upstream regulators; Supplementary Table S3, https://doi.org/10.1667/RADE-20-00041.S4) identified ID2, XBP1, PAX5 (z-scores >2.6) as top regulators of the cluster 1 pattern (345 genes downregulated in all groups after irradiation). Known radiation response transcriptional regulators TP53, HIF1A, TP73 and FOXO3 (z-scores >3.2) were among the proteins implicated as key regulators of cluster 2 (274 genes upregulated early in the time series and returning to background levels at day 14 in all groups). Cluster 3 (249 genes with mRNA levels that steadily increase over the time course in all groups) may be regulated by NRF2, E2F1, E2F4 and RB1 (z-scores >2.5). E2F family members are known cell cycle regulatory genes and may be involved in promoting cell division after radiation. Cluster 4 (43 genes that were upregulated early after irradiation and then switched to downregulated at later times between 7 and 14 days) may be regulated by zinc finger protein, FOG family member 1 (ZFPM1; involved in erythroid and megakaryocytic cell differentiation), BRCA2 and DMTF1. Cluster 5 (106 genes that are upregulated early in the time course and stay upregulated throughout in all groups) may be regulated by TAF7L and TOX (thymocyte selection associated high mobility group box, involved in T-cell development) (z-scores >2). TP53, NFKBIA, MYC, MYCN and HDAC1 were all predicted to be involved in regulation of both cluster 2 and cluster 3 genes, which represent two very different gene response patterns.

This study was designed to accumulate radiation doses of 4 Gy at 3, 5, 7 and 14 days in the different injection groups. The daily dose rates (shown in Table 1) for the 4 Gy doses ranged from 1.3 Gy/day (highest injection group, 9.29 MBq) to 0.2 Gy/day (lowest injection group, 5.71 MBq). We compared genes across all four injection groups at the times corresponding to the 4 Gy cumulative dose and identified 161 genes that were common to all the 4 Gy exposures and with similar expression patterns in all injection groups. Expression of most of these genes showed a gradual downregulation over time in all groups (4 Gy common genes; Supplementary Table S4, https://doi.org/10.1667/RADE-20-00041.S5). Gene Ontology analysis showed that most of these genes were enriched in processes associated with antigen processing, B-cell activation and positive regulation of the immune system (Table 3).

Next, we looked for a dose-rate-independent pattern of gene expression change and identified genes that showed a similar response at 3 Gy, 4 Gy and 5 Gy doses across all groups. This analysis yielded genes that were mostly downregulated at all doses measured. The heat map in Fig. 5A shows the expression levels of top genes, selected by highest magnitude of change in signal, at different dose rates (and times) at 4 Gy. Fc receptor-like A (Fcrla) and histocompatibility 2, class II antigen E beta (H2-EB-1) illustrate this pattern, showing consistent suppression throughout the 14-day time course, as validated by qRT-PCR (Fig. 5B). Gene Ontology analysis using the PANTHER database was applied to the top 700 genes. Biological processes identified as enriched among these dose-rate-independent genes included categories of antigen processing and presentation of MHC class II (Bonferroni P value 1 × 10^–5) and also general immune system processes (Bonferroni P value 3 × 10^–6) (Dose-rate independent genes; see Supplementary Table S4; https://doi.org/10.1667/RADE-20-00041.S5).

Using a method similar to that used for identification of genes that were dose-rate independent, we identified 178 up- or downregulated genes that showed dose-rate dependence [Fig. 6A; Dose-rate-dependent genes; Supplementary Table S4; (https://doi.org/10.1667/RADE-20-00041.1.S5)]. Two of these genes, reactive oxygen species modulator 1 (Romo1) and pleckstrin homology-like domain, family A, member 3 (Phlda3) were measured by qRT-PCR to validate the response pattern (Fig. 6B). PANTHER Gene Ontology analysis of these 178 genes revealed that mediators of immune response (Bonferroni P value 10^–6) and immunoglobulin production (Bonferroni P value 10^–7) were significantly enriched among these genes.

¹³⁷CsCl is highly soluble and is rapidly distributed through the tissues of the mouse, and then eliminated over time following characteristic biokinetics (Fig. 1A). This results in prolonged exposures to low-dose-rate radiation with decreasing exposure rates. We have previously shown that IP injection of ¹³⁷CsCl solution into mice produces unique gene expression patterns over time as well as a more sustained response compared to an equivalent acute exposure (13, 23 ). Internalization of radioactive isotopes from fallout thus has the potential to both influence the extent of radiological injury and complicate the interpretation of dosimetry biomarkers.

We have now monitored the temporal response of genes in mouse blood across a range of initial injected isotope activities, to begin delineating dose and dose-rate effects on gene expression in blood. The dose rate to the mice was highest immediately after the injection, ranging from an average of 0.6 ± 0.11 Gy/day to 1.4 ± 0.08 Gy/day during the first two days after injection, and declining to 0.2 Gy/day to 0.4 Gy/day by the end of the experiment, depending on the amount of activity initially injected (Table 1). In designing this experiment, we selected injection activities to provide dose-rate profiles that would result in a 4 Gy dose to the animals at different times. Sacrifice times were also selected to enable similar comparisons of 3 and 5 Gy doses delivered at different dose rates. We focused on the first two weeks after the start of exposure, with the maximum dose in this study being 12 Gy. Figure 1B shows the four different dose curves and sacrifice times for the experimental groups.

For an initial overview of the gene expression response, we looked at general gene expression response patterns in this study. Thousands of genes were differentially expressed across the study (Table 2), many of which are also significantly changed after acute external radiation exposures (13, 14, 20, 39 ). The number of differentially expressed genes generally increased with time. However, most of the injection groups showed a decrease in the number of differentially expressed genes at the five-day time point. This was consistent with our previously published study where a decrease in the number of differentially expressed genes was also observed at day 5. In that study, which continued up to 30 days after the injection, more than 6,000 genes were differentially expressed at 20 and 30 days after injection, suggesting that the gene expression response may have continued to increase if the current study had been continued past the two-week point. There were a large number of downregulated genes at later times in the time course for each injection group. The high preponderance of underexpressed genes may suggest that not only were there broad signaling changes, but also loss of specific cell subtypes. This would be consistent with studies of non-lethal acute doses in mice, which have been shown to affect the immune system and alter the representation of specific cell subtypes in the circulating blood (40–42 ).

We applied bioinformatics approaches to illuminate the underlying mechanisms and drivers of the response to the protracted exposure and damage from internal ¹³⁷Cs. First, we performed REACTOME pathway analysis using the tools of the PANTHER database (27, 37 ). Statistical overrepresentation analysis and comparison of the gene lists across the study revealed many pathways under broad parent categories, such as gene expression mechanisms and immune system responses, which were involved throughout the two-week study. Some pathways, however, appeared to be implicated in a more time-dependent manner regardless of the amount of injected radioactivity. Surprisingly, the TP53 REACTOME pathway was identified by this analysis as significantly activated only at later times, in contrast to results from the less stringent IPA-based network approach, which suggested TP53 activation was highest at earlier times (Fig. 2). This discrepancy may be due to the difference in algorithms used for creating the databases and correcting for multiple comparisons. REACTOME is based on pathway knowledge and membership in well-studied pathways (37 ), while the IPA network-based approach casts a wider net, being based on multiple interaction types between molecules (38 ). The REACTOME-based approach here provided an overview of changes in signaling that may affect gene expression in this protracted exposure study. In the case of TP53, it is known to be activated in chronic low-dose-rate irradiation in humans, mice and murine cells (12, 43, 44 ). In the IPA analysis of our study TP53 was implicated as an activated transcriptional regulator as early as the two-day time point (Fig. 3) in all injection groups. The TP53 target network consists of diverse genes showing multiple response patterns, the balance of which may be affected by differences in dose and dose rate. In the same IPA analysis, MYC was shown to be inhibited especially at the one-week time point in the various injection groups, indicating that it may be involved in cell turnover or suppression and not restricted to response at a specific dose or dose rate. MYC transcription is downregulated by TP53, and itself down-regulates TP53 activation, playing a necessary role for deciding cell fate after DNA damage stress (45, 46 ). It may also be involved in regulating the timing of the TP53 network activation during chronic irradiation in this study.

Conventional tools for identifying time-series-based gene clusters are not designed for complex studies with more than one variable. Therefore, we used the maSigPro algorithm to identify significant clusters of gene expression response profiles across the four injection groups and five time points in this study, applying a false discovery rate cutoff <0.01. The results shown in Fig. 4 (for gene lists and their IPA analyses see also Supplementary Table S3, https://doi.org/10.1667/RADE-20-00041.S4) revealed three prominent patterns of response. Cluster 2 consists of genes that switch from overexpressed to underexpressed during the two-week time course. This cluster was most interesting, as this was a major expression pattern observed for many TP53-regulated genes in our earlier internal ¹³⁷Cs study (13 ). Consistent with that prior study, TP53 was found to be the top predicted activated transcription factor in IPA upstream analysis of the cluster 2 genes (Supplementary Table S3). The analysis also suggests that these genes may contribute to stress signaling at earlier times and then cell damage processes at later times. Defining such reversals in the signaling response may be critical in elucidating recovery points in individuals with exposure to internal radiation. Clusters 1 and 3, with consistently down- and upregulated genes, respectively, may also be useful to identify genes and cellular mechanisms that can be developed as biomarkers of radiation exposure and response.

To identify genes that might better distinguish between an acute exposure and an ongoing exposure from an internal emitter, we looked for genes that were either dose-rate dependent or dose-rate independent. In this study, a 4 Gy dose was achieved on a different day in each of the injection groups. We compared the lists of genes that responded at 4 Gy radiation (at different dose rates) and identified 161 common genes (Supplementary Table S4; https://doi.org/10.1667/RADE-20-00041.S5), most of which were consistently downregulated immune response genes. Expanding on these observations, we then investigated dose-rate-independent genes across the entire study, searching for genes that showed similar responses at 3, 4 and 5 Gy cumulative doses (regardless of dose rate) and analyzed the top genes for enrichment of biological processes that were related to specific immune functions. The pattern of expression of the majority of these genes (Fig. 5) suggests there may be a dose threshold for mRNA changes beyond which further increasing the dose does not further affect mRNA levels.

Identification of dose-rate-dependent genes at 4 Gy radiation either positively or negatively correlated with dose rate, revealed 178 genes enriched in immune function processes. The results suggest that certain immune processes dominate the functional response at all times in the two-week study and at all doses/dose rates. Two broad patterns of response to exposure are evident from this analysis; one that is changed consistently across a broad range of radiation exposures, and the other that is highly responsive to dose rate. The first response may be driven by transcriptional regulation and mRNA turnover, but may also indicate depletion of specific cell types, such as radiosensitive B cells. The second pattern suggests that there are transcripts that potentially can be used to measure and evaluate dose-rate effects in the blood.

Blood cell subtypes are known to have different radiosensitivities (47 ), which could play a role in driving the gene expression changes observed after radiation exposure. A recent study by Calvi et al. (48 ) characterized blood cell subtypes and hematopoietic repopulation after irradiation by both external total-body irradiation (TBI) and internal soluble ¹³⁷Cs. Female C67BL/6J mice received 6 Gy TBI, with or without a subsequent IP injection of 100 µCi (3.7 MBq) of ¹³⁷Cs to mimic acute exposure from an IND followed by exposure to fallout. Approximately 40% of the animals that received the combined dose (6 Gy TBI with 2.6 Gy at ∼28 days) died from hematopoietic acute radiation syndrome. The survivors were followed for ∼8 months and hematopoietic repopulation assays suggested that the combined TBI with internal ¹³⁷Cs exposure was similar to the TBI, with less effect from the internal ¹³⁷Cs alone. Interestingly, in this study of the hematopoietic niche of the mouse, the excretion kinetics of internally administered ¹³⁷Cs were decreased by TBI. This does not preclude earlier effects in peripheral blood, and in both this and our studies there is clear evidence for persistent and long-term damage to the blood and immune system.

In vivo¹³⁷Cs internal emitter studies have the disadvantage of generating radioactive biofluids that are not amenable to follow-up end points. In the case of blood, directly measuring blood cell sub-populations at early times after ¹³⁷Cs injection is problematic. To shed some light on the potential changes in blood cell subtypes, we used CIBERSORT to deconvolute the microarray data and reconstruct relative fractions of different blood cell subtypes through the study. The results of this analysis suggested that B cells were significantly depleted after injection of ¹³⁷Cs in all groups (Fig. 7) and remained depleted throughout the time course (q values shown in Supplementary Fig. S2, statistical analysis for B cells; https://doi.org/10.1667/RADE-20-00041.1.S1). Other cell types show fluctuations, but relative abundance was not significantly changed by ¹³⁷Cs injection, consistent with the gene expression and ontology results and prior reports that suggest B cells are the most sensitive to internal or low-dose-rate radiation (11, 47, 49 ).

To enable a broader range of experiments relevant to exposure to ¹³⁷Cs from radioactive fallout, our group has developed a low-dose-rate external beam irradiator that can match the dose-rate profile of internally deposited ¹³⁷Cs (50 ). This approach eliminates the generation of radioactive biological specimens, simplifying downstream assays and lowering the cost of experiments. Using this approach, we plan to more fully characterize and validate genes identified here as potentially useful biomarkers of internal radiation exposure, dose and dose rate. These experiments will also allow direct immune-phenotyping of blood cell sub-populations to verify the in silico results throughout the time course of exposure.

In this study we investigated the molecular response to internal radiation using mice injected with ¹³⁷Cs. The goal was to evaluate dose and dose-rate effects on gene expression and the potential confounding effect of prolonged irradiation on the biological response of the animal. The results revealed both a sustained effect of low-dose-rate radiation on immune response and an evolution of signaling over time, with early processes related to blood cell damage of platelets and B cells; later processes related to cell turnover, persistent DNA damage response and RNA processing. Gene expression results, which indicate that B cells are indeed the most vulnerable cell types to internal radiation, support observations from nuclear disasters regarding the long-lived effects of radiation and impact on health (51 ). Taken together with the results of cytogenetic analyses, which can detect and measure internal radioactivity (22, 23 ) very rapidly after the start of exposure, persistent transcriptomic changes will be useful in evaluating prolonged exposures and may contribute to dose reconstruction (52 ). Realistic scenarios relevant to biodosimetry will likely involve exposures consisting of prompt irradiation, delayed irradiation due to fallout and groundshine, and further complication by internally ingested or inhaled radionuclides. Our work in measuring and understanding the biological response in vivo is a considerable step towards the goal of practical radiation biodosimetry.

Fig. S1. Daily dose rates for all samples.

Fig. S2. Relative abundance of NLR, T cells, B cells and NK cells. Statistical analysis is shown for B cells.

Table S1. Significantly differentially expressed genes of all subgroups compared with time-matched controls, Venn diagram by day and group, and IPA upstream regulators.

Table S2. Gene Ontology analyses of different subgroups. PANTHER Gene Ontology biological processes and REACTOME table.

Table S3. maSigPro results and gene lists for clusters. IPA upstream regulators: maSigPro clusters 1–5, upstream regulator predictions in IPA, z-score cut off 1.7.

Table S4. Gene expression patterns, 4 Gy common genes, Dose-rate-independent and -dependent genes.

We thank Johanna Steinbrecher and the Lovelace Biomedical Research Institute technical team for assistance with mouse and sample handling for this study. We thank Dr. Igor Shuryak for consultations on the biostatistical analyses of gene expression data. Analyses were performed using BRB-ArrayTools, developed by Dr. Richard Simon and the BRB-ArrayTools Development Team. This work was supported by the Center for High-Throughput Minimally-Invasive Radiation Biodosimetry, National Institute of Allergy and Infectious Diseases grant no. U19 AI067773.