A radiological/nuclear (RAD-NUC) incident, especially in an urban setting, results in diverse radiation-induced injuries due to heterogeneities in dose, the extent of partial-body shielding, human biodiversity and pre-existing health conditions. For example, acute radiation syndrome (ARS) can result in death within days to weeks of exposure to 0.7–10 Gy doses and is associated with destruction of the bone marrow, known as hematopoietic ARS (H-ARS). However, partial-body shielding that spares a portion of the bone marrow from exposure can significantly reduce the occurrence of H-ARS, but delayed effects of acute radiation exposure (DEARE) can still occur within months or years after exposure depending on the individual. In a mass casualty event, ideal triage must be able to pre-symptomatically identify individuals likely to develop radiation-induced injuries and provide an appropriate treatment plan. Today, while there are FDA approved treatments for hematopoietic ARS, there are no approved diagnosis for radiation injury and no approved treatments for the broad spectra of injuries associated with radiation. This has resulted in a major capability gap in the nations preparedness to a potentially catastrophic RAD-NUC event. Circulating microRNA (miRNA) are a promising class of biomarkers for this application because the molecules are accessible via a routine blood draw and are excreted by various tissues throughout the body. To test if miRNA can be used to predict distinct tissue-specific radiation-induced injuries, we compared the changes to the circulating miRNA profiles after total-body irradiation (TBI) and whole thorax lung irradiation (WTLI) in non-human primates at doses designed to induce ARS (day 2 postirradiation; 2–6.5 Gy) and DEARE (day 15 postirradiation; 9.8 or 10.7 Gy), respectively. In both models, miRNA sequences were identified that correlated with the onset of severe neutropenia (counts <500 µL–1; TBI) or survival (WTLI). This method identified panels of eleven miRNA for both model and assigned functional roles for the panel members using gene ontology enrichment analysis. A common signature of radiation-induced injury was observed in both models: apoptosis, DNA damage repair, p53 signaling, pro-inflammatory response, and growth factor/cytokine signaling pathways were predicted to be disrupted. In addition, injury-specific pathways were identified. In TBI, pathways associated with ubiquitination, specifically of histone H2A, were enriched, suggesting more impact to DNA damage repair mechanisms and apoptosis. In WTLI, pro-fibrotic pathways including transforming growth factor (TGF-β) and bone morphogenetic protein (BMP) signaling pathways were enriched, consistent with the onset of late lung injury. These results suggest that miRNA may indeed be able to predict the onset of distinct types of radiation-induced injuries.
INTRODUCTION
Effective triage after a mass casualty RAD-NUC incident would benefit from an assay capable of identifying individuals likely to develop radiation-induced injuries at a time when therapeutic intervention would be efficacious. The Food and Drug Administration (FDA) approved the G-CSF-based drugs Neupogen and Neulasta in 2015, and more recently, Leukine (2018) and Nplate (2021) to treat individuals exposed to myelosuppressive doses of radiation (1). However, the onset of hematopoietic ARS (H-ARS) at a given dose is severely reduced by partial bone marrow shielding (2–4), so without effective, injury-specific biodosimetry, there is a risk of unnecessarily treating individuals that may not benefit from these drugs. In other words, how an individual's body responds to ionizing radiation exposure not only depends on the dose received but also the parts of the body that were irradiated, the presence of pre-existing conditions, diet, genetic makeup, and a range of other variables often termed ‘biodiversity’. While many research groups have focused on developing assays for determining the dose of ionizing radiation received, we propose here a more holistic and system-centric approach where we predict injuries that result from exposure to ionizing radiation. We believe an assay that can predict tissue-specific injury from exposure to ionizing radiation, from the time of triage to long-term surveillance of survivors, would provide a clinician with the tools needed to tailor treatment planning through the continuum of care.
To achieve this, we, and others, have shown that microRNA (miRNA), short (∼22 nt), non-coding RNA sequences that post-translationally regulate gene expression (5), may be useful biomarkers for biodosimetry and predicting radiation-induced injuries (6–13). In a triage scenario, circulating miRNA, miRNA sequences produced by various tissues in the body and released into circulation in exosomes or protein complexes (14, 15), have the advantage of being accessible via a simple blood draw and can be assayed using conventional techniques, such as RT-qPCR (16).
One of the challenges associated with establishing miRNA panels for radiation-induced injury is assigning functional roles for the identified sequences. A miRNA can potentially regulate several, perhaps hundreds of genes, and a single gene can be regulated by several miRNA. Moreover, circulating miRNA may affect gene expression in cooperation with each other. As a result, approaches for validating potential miRNA biomarkers that examine the functional roles individually may not be effective. Instead, understanding the expected biological activity of the set of miRNAs may be more insightful. One tool to do this is gene ontology enrichment analysis (GOEA). GOEA, like similar pathway analysis tools, assigns genes to functional roles (GO terms, e.g. signal transduction) (17). GO terms are organized in a graph, where each GO term is a node and the relationship between the terms are edges. GO is hierarchical with parent terms being broader than child terms. For example, signal transduction is the parent term of cell surface receptor signaling. This structure enables identifying the most specific functional role for the set of miRNA. While some miRNAs have been functionally annotated, the GO annotations remain incomplete. Instead, miRNA can be mapped to their gene targets, using interaction databases (e.g. miRWalk), and the genes can be used to identify enriched GO terms where the frequency of genes in the data set associated with a GO term of interest is compared to the frequency of all organismal genes for the GO term (18, 19). P values are calculated for the GO terms using Fisher's exact test. Using this algorithm, GOEA has been used to identify and rank potential biological functions for sets of miRNA in a variety of contexts (20–24).
To determine if circulating miRNA can be used to differentiate between the various types of radiation-induced injuries, we compared changes to the circulating miRNA profile in non-human primates (NHP) after total-body irradiation (TBI) at H-ARS-inducing doses (2–6.5 Gy) to whole thorax lung irradiation (WTLI) at DEARE-inducing doses (9.8 or 10.7 Gy). Circulating miRNA were isolated from plasma taken before the onset of symptoms, at day 2 post-TBI and day 15 post-WTLI. NHP exposed to TBI developed severe neutropenia (neutrophil counts < 500 µL–1) by day 7 and NHP exposed to WTLI developed radiation-induced lung diseases (including pneumonitis and fibrosis), that resulted in reduced aerated lung volume by day 60 by computed tomography (CT), increased collagen deposits in the lung tissue by histochemistry, and reduced survival (25). Circulating miRNA that correlated with the onset of post-TBI neutropenia and post-WTLI survival were identified. Predicted functional roles for the identified miRNA were assigned using GOEA.
MATERIALS AND METHODS
Total-Body Irradiation Study
Male and female rhesus macaques (Macaca mulatta) NHPs were exposed to 0 (N = 4), 2 (N = 6), 5.5 (N = 6), or 6.5 Gy (N = 7). All animals underwent routine quarantine and health checks before inclusion into the study as described previously (4, 26). Animal housing and care was compliant with the Animal Welfare Act and recommendations set forth in the Guide for the Care and Use of Laboratory Animals (National Research Council 2011). All procedures were carried out in compliance with institutional policies on the Humane Care and Use of Laboratory Animals and were approved by the Institutional Animal Care and Use Committee (IACUC) of CiToxLAB North America (study number 2011-0603). Prior to inclusion into the study, all animals were acclimated to the facility for 14 days. Animals found to have abnormal baseline complete blood count (CBC), evidence of lung disease, injury or poor heath were excluded from the study. The day of irradiation was defined as day 0. At the time of irradiation, the age of the animals ranged from approximately 2 to 5 years. On the day prior to irradiation (day –1), the body weights ranged from 3.4 to 5.9 kg and from 3.8 to 6.5 kg, for males and females, respectively. All animals were declared healthy before inclusion and were evaluated for mortality and clinical signs (twice daily). Detailed examinations, including body weight and temperature, were performed, prior to animal assignment, during the week prior to irradiation, and on days 1, 3, and 7. On day 0, the animals were subjected to a single uniform total-body dose of gamma radiation from a 60Co source (Theratron 1000), with a targeted dose rate of 100 cGy · min–1. Sequential anteroposterior and posteroanterior exposures with whole-body irradiation were used. Animals were not sedated during irradiation and were breathing ambient air. Dosimeters (nanoDot™; Landuer Glenwood, IL) were placed on the animals for confirmation of exposure, however there was no real-time dosimetry performed via Farmer's chamber. In this study, all animals reached scheduled termination on day 7. Blood samples were collected, without sedation or anesthesia, by venipuncture into a Vacutainer® containing EDTA as an anti-coagulant for hematology three occasions prior to irradiation (within 2 weeks prior to irradiation) and daily for the first seven days postirradiation. Plasma was isolated and frozen (–80°C) prior to miRNA extraction and analysis.
Whole Thorax Lung Irradiation Study
For this study, l4 males (3.93 ± 0.18 kg body weight, 3–4 years old) and 14 females (4.05 ± 0.26 kg body weight, 3–4 years old) rhesus macaques (Macaca mulatta; Chinese origin) were used. All animals underwent routine quarantine and health checks before inclusion into the study as described previously (4, 26). Animal housing and care was compliant with the Animal Welfare Act and recommendations set forth in the Guide for the Care and Use of Laboratory Animals (National Research Council 2011). All procedures were carried out in compliance with institutional policies on the Humane Care and Use of Laboratory Animals and were approved by the Institutional Animal Care and Use Committee (IACUC) of Alta Sciences (study number SNBL 1089.01). Prior to inclusion into the study, all animals were acclimated to the facility for 14 days. Animals found to have abnormal baseline complete blood count (CBC), evidence of lung disease, injury or poor heath were excluded from the study. During acclimation period, baseline ungated thoracic CT scan was acquired per CT-scan acquisition procedure.
Animals were randomly assigned to receive 9.8 Gy (LD20/180; 6 males, 6 females) or 10.7 Gy (LD75/180; 8 males, 8 females) WTLI at a dose rate of 0.60 ± 0.10 Gy min–1. Prior to WTLI, animals were administered antiemetic (Ondasentron HCL, 1.9 mg kg–1) and anesthetized with Ketamine/Xylazine. Sedated animals were restrained on the linear accelerator (LINAC) couch. Radiation was delivered using a 6 MV photon source LINAC (Varian CLINAC 21EX, Varian Associates, Palo Alto, CA) in an anterior-posterior (AP) and posterior-anterior (PA) technique with approximately 50% dose contribution from both the AP and PA beams (4, 27).
After irradiation, cage-side clinical observations were performed twice daily. Observations included non-sedated respiratory rate (NSRR), activity, stool consistency, vomiting and posture. In addition to clinical observations, each animal was weighed twice during the acclimation period and weekly postirradiation. Rectal pulse oximetry was performed twice weekly post irradiation on non-sedated animals. Blood samples for hematology, and serum chemistry were collected by venipuncture in EDTA containing tubes at day –2, and every 3 days postirradiation through day 30, days 40, 50 and 60, and then every 30 ± 3 days, thereafter. Blood samples for miRNA analysis were collected before irradiation, and at days 6, 12, 15, 30, 60 and 90 postirradiation.
Animals that developed tachypnea and respiratory distress indicative of pneumonitis, defined as NSRR >80 breaths per minutes, received medical supportive therapy (dexamethasone) following a scheduled taper described below (4, 28). In the schedule, the first course of treatment was administered intramuscularly (IM) using the following regimen: 1.0 mg kg–1 twice daily (BID) on the first day of treatment; 0.5 mg kg–1 BID for 3 days; 0.5 mg kg–1 once daily (SID) for 3 days; 0.5 mg kg–1 every other day for three doses. For animals that developed respiratory distress within 7 days of stopping the first course of treatment, the second course of treatment was as follows: 1.0 mg kg–1 SID on the first day of treatment; 0.5 mg kg–1 SID for 4 days; 0.5 mg kg–1 every other day for 10 doses. For animals that developed respiratory distress outside 7 days of stopping the first course of treatment, the second course of treatment was as followed: 1.0 mg kg–1 BID on the first day of treatment; 0.5 mg kg–1 BID for 3 days; 0.5 mg kg–1 SID for 3 days; 0.5 mg kg–1 every other day for 10 doses. For animals that completed two courses of dexamethasone treatment and required additional medical supports, a veterinary consultation was obtained and treatment with dexamethasone was conducted by veterinary recommendation. Animals with repeated episodes of respiratory distress and non-responsive to treatment met criteria for euthanasia. Onset of radiation-induced pneumonitis was measured by the timing to first dexamethasone administration, displayed as mean ± standard error of the mean.
Euthanasia criteria included any of the single criterion or two or more of the combination criteria. Single criterion included observation of (i) respiratory distress (e.g., open mouth breathing, cyanotic appearance, labored breathing, respiratory rate ≥ 80 breaths/min) and non-responsive to dexamethasone treatment, (ii) inactivity (e.g., recumbent in the cage for at least 15 min, or non-responsive to touch), (iii) uncontrolled hemorrhage from any orifice, (iv) unrelieved pain or distress following administration of two consecutive increased doses of buprenorphine (0.02 mg/kg IM BID), or (v) severe dehydration. Combination criterion included two or more observations of (i) hyperthermia (rectal temperature ≥41°C), (ii) hypothermia (rectal temperature ≤35°C) weight loss (≥25% of baseline) for two consecutive days, (iii) severe injury or condition (e.g., minor bone fracture, progressive tissue necrosis, non-healing wound), or (iv) complete anorexia for 48 h. Surviving animals were euthanized at day 270 postirradiation. Animals that met the criterion for euthanasia were sedated and blood was collected for hematology, serum chemistry and biomarker discovery. If possible, CT scan was acquired prior to blood collection. At euthanasia, animals were weighed and euthanized by intravenous (IV) injection of a pentobarbital overdose, followed by necropsy.
Computer tomography (CT) scans and image analysis
To evaluate the severity of the radiation-induced lung injury, serial computer tomography (CT) scans were acquired every 30 ± 3 days postirradiation until endpoint (day 270) or until IACUC criteria for euthanasia were met. Non-contrast CT scans were acquired with the animal sedated and positioned supine, headfirst into the scanner, arms overhead and restrained on the table. A spiral CT scan (SomoTom Scope, Siemens Medical Solutions, USA, Inc) with slice thickness of 1.5 mm was acquired with a lung window reconstruction using 80 kV potential, 111 mAs, 1.0 rotation time, 0.4 pitch level and B70s kernel. CT scans were read and analyzed by an experienced veterinarian radiologist on the ClearCanvas Workstation (Version 2.0.12729.37986 SP1; Synaptive Medical, Toronto, Canada). To calculate lung volumes and quantify extent of lung injury, the radiologist performed a three-dimensional volumetric reconstruction of the thoracis CT radiographs with 3DSlicer (Version 4.7.0-2017-08-12 r26243; http://www.slicer.org) using a semiautomated radiodensity (Hounsfield units, HUs) thresholding (4, 27): aerated lung volume (–880 to –220); active pneumonitis and fibrosis (–300 to +88). Portions of lung boundaries that were falsely included in the initial lung boundaries were manually excluded and portions of the lung boundaries that were falsely excluded were manually included. Pleural effusions and atelectasis were separately identified manually. Data are presented as aerated lung volumes, excluding normal soft tissue, vasculature and areas of pneumonitis/fibrosis from the volume calculation.
Histopathology
To assess the degree of tissue damage, lung and heart tissues were collected, trimmed, formalin-fixed, paraffin-embedded, sectioned, and stained with hematoxylin and eosin (H&E) or Masson's Trichome (MT). Sections of the right and left cranial and caudal lung lobes, the right medial lung lobe, and longitudinal sections of the right and left ventricle and atria and the septum of the heart, through the middle of these structures, were taken. Slides were then microscopically examined and scored by an experienced board-certified veterinarian pathologist for inflammatory cell infiltrates and collagen deposition in irradiated tissue. Histopathological alterations were quantified following methods described by Parker et al. (2). Briefly, to convey the extent of histological alteration, each observation was given a severity score from 0 to 4, corresponding to undiscernible, minimal, mild, moderate and marked, respectively (2).
Sequencing
Plasma was checked for excess hemolysis using UV-vis spectroscopy and samples with an absorbance value greater than 1.2 A.U. at 415 nm, corresponding to 0.3% hemolysis (29), were excluded from the study. Circulating plasma miRNA was isolated from 100 µL of plasma using the miRNEasy Serum/Plasma Advanced Kit (Qiagen). Libraries for sequencing were prepared using the QIAseq miRNA Library Kit (Qiagen), with 5.8 µL of miRNA extracts as input, a 1:10 dilution of the 3′-adaptor, a 1:5 dilution of the 5′-adaptor, a 1:10 dilution of the RT primer, and 22 amplification cycles. The concentrations of the prepared libraries were determined via Bioanalyzer analysis (2100 Electrophoresis Bioanalyzer, Agilent). Libraries with an adaptor dimer peak (∼160 nt) at least five times greater than the library peak (∼180 nt) were not sequenced. miRNA counts for 2 nM samples were determined via next-generation sequencing (NextSeq 550, Illumina) using 76 read cycles. Demultiplexing, trimming (read lengths between 18–40 bp, 5′-end base quality ≥30, read score ≥20, and 3′-end adaptor sequence to trim of AACTGTAGGCACCATCAAT), and miRNA alignment (using “Homo sapiens/hg19” as the species) was performed using Base-Space (Illumina), using the Small RNA v1.0.1, FASTQ Toolkit v2.2.0, and FASTQ Generation v1.0.0. Sequencing runs with less than 100,000 miRNA reads were rejected.
Analysis
Raw sequencing counts were normalized by total library size to obtain the reads for a given sequence per million total reads (RPM), then by quantile normalization of the log2 RPM. miRNA that correlated with neutropenia (TBI, measured at day 2 postirradiation) or survival (WTLI, measured at day 15 postirradiation) were identified by splitting the animals into two cohorts based on the relative abundance of a given miRNA. Animals that expressed above the median value were placed into one cohort and animals that expressed the miRNA at or below the median value were placed into the second cohort. Kaplan-Meier survival curves of the two cohorts were compared using the log rank test (Mantel-Cox) (30), and uncertainties (95% confidence intervals) estimated according to Greenwood (31). For miRNA where the log rank P value was less than 0.05 between the cohorts were used for GOEA. Blood count data were compared over time using an ANOVA model, and P values were obtained by comparing postirradiation values to pre-irradiation values using Dunnett's method. Statistical differences between prevalence score distributions were quantified using the Mann-Whitney U test. GOEA was performed by identifying miRNA-gene interactions using miRTarBase (32, 33) and miRWalk (34, 35) (SI file 1). For each GO term in the “biological_process” namespace, the genes associated with the GO term were identified using Homo sapiens GO annotations ( http://current.geneontology.org/products/pages/downloads.html). Enrichment was calculated using Fisher's exact test, as previously described (36).
RESULTS
TBI Results in Severe Neutropenia at High Doses
Total-body irradiation resulted in a decrease in neutrophil counts for every dose ≥2 Gy (Fig. 1). Severe neutropenia by 7 days post-TBI was observed in animals that received doses ≥5.5 Gy, with a frequency proportional to the dose. No animals that received 2 Gy TBI, half of the animals that received 5.5 Gy, and most (5 out of 7) animals that received 6.5 Gy developed severe neutropenia.
FIG 1
TBI at doses ≥2 Gy results in loss of neutrophils and severe neutropenia by 7 days postirradiation. Sham-irradiated animals did not show any significant changes in neutrophil counts, while 2 Gy irradiated animals showed a significant drop in neutrophils (ANOVA P value = 5.0 × 10–5), but no instances of severe neutropenia (counts <500 µL–1) were observed. 3 out of 6 animals that were exposed to 5.5 Gy, and 5 out of 7 animals exposed to 6.5 Gy developed severe neutropenia. This dose-dependent response is consistent with our previous findings (6). NHPs exposed to 9.8 or 10.7 Gy WTLI did not develop severe neutropenia by day 7 (25).
WTLI Results in Loss of Aerated Lung Volume
WTLI did not result in a significant decrease in neutrophil counts, which remained constant within the first few weeks postirradiation followed by an increase by month 2 ( Fig. S1 (547_rare-196-06-06_s01.pptx); https://doi.org/10.1667/RADE-21-00043.1.S1). This increase is consistent with neutrophil infiltration of the lung due to radiation-induced injury and correlates with the observation of loss of aerated lung volume by CT (Fig. 2) accompanied by increased respiratory distress as measured by changes to the respiratory rate that may need dexamethasone intervention ( Fig. S2 (547_rare-196-06-06_s01.pptx); https://doi.org/10.1667/RADE-21-00043.1.S1). CT scans showed significant increases in density within the lung, indicating a loss in aerated lung volume, at two months post-WTLI for animals exposed to 10.7 Gy (0.7-fold change, P value = 0.0004), and at three months post-WTLI for animals exposed to 9.8 Gy (0.8-fold-change, P value = 0.009). Aerated lung volume continued to decrease until 4–5 months post-WTLI, where the loss stabilized and remained low.
FIG. 2
WTLI resulted in loss in aerated lung volume and eventual pulmonary and cardiac fibrosis. Panel A: CT scans taken prior to irradiation and two months post-WTLI show a pronounced, dose-dependent loss in aerated lung volume. The red contour lines indicate boundary areas of lung injury (pneumonitis/fibrosis) using the HU range of –300 to +80. The green contour lines indicate areas of normal lung. Panel B: Aerated lung volume after WTLI was reduced for both dose cohorts, with animals exposed to 10.7 Gy WTLI presenting a statistically significant loss in aerated lung volume at day 60 (0.7-fold change, P = 0.0004), compared to day 90 for animals that received 9.8 Gy (0.8-fold change, P = 0.009). Panels C and D: H&E and MT staining of lung and heart tissue extracted after the death of the animal. MT staining showed significant fibrosis in the lungs of both animals and in the heart tissue of animals exposed to 10.7 Gy. Panel E: Quantification of the prevalence of the indicated markers of heart and lung tissue from both dose cohorts, where 0 indicated the marker was not present and 5 indicates that the marker was ubiquitous. In the heart, dose-dependent (Mann-Whitney U test P < 0.05) infiltration of mononuclear cells, myocardium degeneration and interstitial fibrosis was observed. In the lung, dose-dependent (Mann-Whitney U test P < 0.05) infiltration of neutrophils into the alveolar space and interstitial invasion of mononuclear cells was observed.
Identification of miRNA Corresponding to the Onset of Neutropenia (TBI) or Survival (WTLI)
Radiation-induced changes in the circulating miRNA profiles of NHPs exposed to either TBI or WTLI were identified using a survival-based analysis with the onset of severe neutropenia (TBI) or death (WTLI) as the endpoint. For each model, 11 miRNAs were identified where there was a statistically significant difference in neutropenia onset or survival (log rank P value <0.05) based on the abundance of the miRNA. In TBI, let-7a-5p, miR-16-2-5p, miR-93-5p, miR-122-5p, miR-186-5p, miR-193a-5p, miR-223-3p, miR-320a-3p, miR-375-3p, miR-423-5p, and miR-627-5p were identified from plasma isolated at day 2 postirradiation (Fig. 3). Similarly, miR-29c-3p, miR-181a-5p, miR-185-5p, miR-191-3p, miR-200c-3p, miR-206, miR-301a-5p, miR-320d, miR-340-3p, miR-433-3p, and miR-496 were identified in animal exposed to WTLI (Fig. 4).
FIG. 3
The relative abundance of certain circulating miRNA isolated two days post-TBI was associated with the onset of severe neutropenia. For a given miRNA, NHPs were divided into two cohorts based on expression above (blue) or below (orange) the median abundance. For the indicated miRNA, the difference in survival between these two cohorts was statistically significant (log rank test P < 0.05), suggesting a role for these miRNAs as potential biomarkers for predicting severe neutropenia.
FIG. 4
The relative abundance of certain circulating miRNA isolated at day 15 post-WTLI was associated with reduced survival. For a given miRNA, NHPs were divided into two cohorts based on expression above (blue) or below (orange) the median abundance. For the indicated miRNA, the difference in survival between these two cohorts was statistically significant (log rank test P < 0.05), suggesting a role for these miRNAs as potential biomarkers for early identification of late organ injury following exposure to ionizing radiation.
Gene Ontology Enrichment Analysis
The biological effects of the identified miRNAs were estimated using GOEA. The set of TBI miRNA resulted in 829 enriched (Fisher's exact P value <0.05) GO pathways, and the set of WTLI miRNA resulted in 990 enriched GO pathways. 676 GO pathways were common to both models (Fig. 5). The top enriched GO terms were clustered using keywords (e.g., apoptosis or interleukin) or sub-terms (e.g., DNA damage, p53 class mediator). The frequencies these keywords and sub-terms appeared in the enriched GO terms versus all GO terms was calculated (Fig. 6). Apoptotic signaling pathways were the most enriched terms for both TBI and WTLI. Similarly, GO terms containing the sub-terms p53 class mediator, DNA damage, NF-κB, and toll-like receptor displayed similar frequencies of enrichment after both TBI and WTLI. This analysis also identified keywords/sub-terms that were only enriched in TBI or WTLI. Ubiquitination, and specifically ubiquitination of histone H2A, was only observed after TBI. In WTLI, GO terms associated with bone morphogenetic protein (BMP) and transforming growth factor-β (TGF-β) signaling were enriched. As GO terms are organized in a directed acyclic graph (DAG), where more specific GO terms are related to more general terms using a parent-child relationship, called an “is_a” relationship. To better understand the role of the enriched sub-terms, we visualized the sub-DAG for all GO terms containing the enriched keyword/sub-term and all parent terms. The nodes of the sub-DAGs, representing specific GO terms, were colored based on the Fisher's exact test P value.
FIG. 5
Venn diagram of enriched GO terms associated in miRNA corresponding with the onset of severe neutropenia (TBI) and DEARE survival (WTLI). The majority (676) of enriched GO terms (Fisher's exact test P < 0.05) were found both after TBI or WTLI. TBI resulted in 153 unique GO terms, and WTLI resulted in 314.
FIG. 6
Frequency of keywords or sub-terms for enriched GO terms after TBI or WTLI. To group GO terms by function, keywords or sub-terms were identified to cluster GO terms into functional groups. Frequencies for each keyword/sub-term were calculated by dividing the number of occurrences of the sub-term in the enriched data set by the total number of occurrences in all GO terms. The most enriched GO categories were associated with known responses to radiation injury (apoptosis, p53 class mediator, DNA damage, NF-κB, etc.). Some categories of GO terms were enriched in the TBI (ubiquitination, histone H2A) and WTLI (BMP signaling pathway, transforming growth factor) cohorts.
Ubiquitination
GO terms associated with ubiquitination were enriched in TBI but not WTLI (Fig. 7). Three parent GO terms were enriched by the TBI miRNA: actin ubiquitination (GO:0007014, P value = 0.0016), histone ubiquitination (GO:0016574; P value = 0.049), and protein polyubiquitination (GO:0000209, P value = 0.013). In addition, enriched child GO terms of histone ubiquitination included histone H2B ubiquitination (GO:0033523, P value = 0.049), histone H2B conserved C-terminal lysine ubiquitination (GO:0071894, P value = 0.049), histone H2A ubiquitination (GO:0033522, P value = 0.049), histone H2A-K15 ubiquitination (GO:0036352, P value = 0.049), histone H2A-K13 ubiquitination (GO:0036351, P value = 0.049), histone monoubiquitination (GO:0010390, P value = 0.039), histone H2A monoubiquitination (GO:0035518, P value = 0.028), histone H2A-K119 monoubiquitination (GO:0036353, P value = 0.022), and histone H2A K63-linked monoubiquitination (GO:0070535, P value = 0.0017), which was the only ubiquitination-related GO term also enriched in the WTLI miRNA (P value = 0.044). Enriched child GO terms of protein polyubiquitination included protein K29-linked ubiquitination (GO:0035519, P value = 0.013), protein K27-linked ubiquitination (GO:0044314, P value = 0.013), protein K48-linked ubiquitination (GO:0070936, P value = 0.0086), protein K11-linked ubiquitination (GO:00790979, P value = 0.015), protein K6-linked ubiquitination (GO:0085020, P value = 0.013), protein linear polyubiquitination (GO:0097039, P value = 0.014), protein K33-linked ubiquitination (GO:1990390, P value = 0.013), and protein K63-linked ubiquitination (GO:0070534, P value = 0.0086).
BMP Signaling Pathway
GO terms associated with BMP signaling pathway were enriched in WTLI but not TBI (Fig. 8). Both types of radiation injury produced a miRNA response that implicated significant changes in signal transduction (GO:0007165, TBI P value = 0.019, WTLI P value = 0.0033). However, only the WTLI miRNA implicated significant changes to BMP signaling pathway (GO:0030509, P value = 0.048), and the child GO terms specification of steric bud anterior/ posterior symmetry by BMP signaling pathway (GO:0072101, P value = 0.03), BMP signaling pathway involved in ureter morphogenesis (GO:0061151, P value = 0.048), BMP signaling pathway involved in renal system segmentation (GO:0061151, P value = 0.048), BMP signaling pathway involved in nephrite duct formation (GO:0071893, P value = 0.048), BMP signaling pathway involved in heart development (GO:0061312, P value = 0.047), and negative regulation of branch elongation involved in ureteric bud branching by BMP signaling pathway (GO:0072097, P value = 0.023). Both the TBI and WTLI enriched GO terms included positive regulation of Wnt signaling pathway by BMP signaling pathway (GO:0060904, TBI P value = 0.03, WTLI P value = 0.0081).
Transforming Growth Factor
Like BMP signaling pathway, GO terms associated with the sub-term “transforming growth factor” were more enriched in WTLI than TBI. In both TBI and WTLI, signal transduction, and positive regulation of cytokine secretion (GO:0050715, TBI P value = 0.0076, WTLI P value = 0.00012). However, its child GO term, positive regulation of transforming growth factor-beta secretion (GO:2001203), is only predicted to be enriched by the WTLI miRNA (P value = 0.0001). In addition, the WTLI miRNA are associated with negative regulation of cytokine secretion (GO:0050710, P value = 0.00094) and its child negative regulation of transforming growth factor-beta secretion (GO:2001202, P value = 0.00091). In the TBI, but not WTLI, GO terms, transforming growth factor-beta secretion (GO:0038044, P value = 0.05) was enriched. Negative (GO:0010629, P value = 0.047) and positive (GO:0010628, P value =0.019) regulation of gene expression are enriched in the WTLI GO terms, and these terms are linked to positive (GO:0001819, P value = 0.003) and negative (GO:0001818, P value = 0.014) regulation of cytokine production, positive (GO:0071636, P value = 0.0025) and negative (GO:0071635, P value = 0.013) regulation of transforming growth factor beta production, positive (GO:0032914, P value = 0.0033) and negative (GO:0032911, P value = 0.011) regulation of transforming growth factor beta1 production, positive (GO:0032916, P value = 0.0023) and negative (GO:0071635, P value = 0.013) regulation of transforming growth factor beta3 production, negative regulation of transforming growth factor beta2 production (GO:0032912, P value = 0.01), and positive regulation of transforming growth factor beta1 activation (GO:1901394, P value = 0.0037). Transforming growth factor beta receptor signaling path way (GO:0007179, P value = 0.01) and its child, transforming growth factor beta receptor signaling pathway involved in heart development (GO:1905313, P value = 0.01), were also enriched.
DISCUSSION
Both TBI and WTLI resulted in the expected phenotypical outcomes associated with H-ARS and DEARE, respectively. TBI resulted in higher incidence of severe neutropenia (neutrophil counts <500 µL–1) by day 7 in a dose-dependent manner (Fig. 1). WTLI resulted in loss in aerated lung volume and respiratory distress, due to pneumonitis and/or fibrosis, leading to death (Fig. 2). These outcomes have been previously described in independent studies in NHPs (6, 25), mice (7), and other species. Thus, the outcomes described here are consistent with the well-described response to these types of radiation exposures.
Circulating miRNA were identified in both models by survival analysis using the time from irradiation to severe neutropenia for NHP exposed to TBI and survival time for NHP exposed to WTLI (Figs. 3 and 4). For each miRNA, the animals were divided into two cohorts based on the abundance of the miRNA compared to the median abundance. miRNA where the difference in survival between the cohorts was significant (log rank P value <0.05) were used for subsequent GOEA. Using this method, 11 miRNAs were identified for both TBI and WTLI cohorts. This consistent list of miRNA facilitated the comparison of the results of GOEA, as different numbers of miRNA may change the number of miRNA-associated genes.
Several identified miRNAs have been previously shown to be correlated with radiation-induced outcomes. For example, miR-375 and miR-16-2-3p have been identified to be related to survival after TBI in NHPs (9). miR-375 has been shown to induce ROS and apoptosis (37–39), and miR-16-2-3p has been shown to inhibit cell proliferation and migration, and induce apoptosis (40). Similarly, miR-181 and miR-200 family members have been shown to be differentially expressed in mice following WTLI (7). Moreover, miR-181a has been implicated in myocardial fibrosis in rats (41) and miR-200 family members have been implicated in pulmonary fibrosis (42), with TGF-β signaling associated with both outcomes. Thus, several miRNA identified here have been previously shown to be associated with specific radiation injuries and mechanistically linked to the respective outcomes.
Using these sets of miRNAs, GOEA was performed for each type of radiation injury. For both radiation models, the most significant pathways were associated with gene transcription (Table 1). This is not surprising, as radiation-induced DNA damage, as well as cytokine and growth factor signaling pathways, are expected to result in changes in transcription. In both radiation models, the most significant GO term was the SREBP signaling pathway (GO:0032933; TBI P value = 3.00 × 10–28, WTLI P value = 3.49 × 10–30). This pathway is associated with changes in lipid homeostasis (43) which has been shown to occur after irradiation (44). In a WTLI study in mice, we have shown that miRNAs are linked to lipid metabolism using KEGG pathway analysis (7, 25). Many of the other most significantly enriched GO terms are associated with regulation of transcription from RNA polymerase II in response to stress and were common to both types of radiation injury. Thus, the most significantly enriched GO terms are common to radiation response in general, and are not specific to the eventual radiation injury.
TABLE 1
Top GO Terms Enriched by miRNA that Correlate to the Onset of Neutropenia (TBI) and Survival (WTLI)
Indeed, most of the enriched GO terms were found in both the TBI (82%) and WTLI (68%) datasets (Fig. 5). To understand which pathways differentiated TBI-induced neutropenia from WTLI-induced pneumonitis and fibrosis we clustered terms by a keyword or phrase (e.g., NF-κB or p53 class mediator), which conveyed relevant biological information. Occurrences of these keywords were counted in the TBI and WTLI datasets and divided by their occurrence in all GO terms to get the frequency each term appeared in the dataset (Fig. 6). For example, over 50% of the GO terms containing the phrase “apoptotic signaling pathway” were predicted to be dysregulated by the TBI or WTLI miRNA. This analysis indicates that apoptosis (apoptotic signaling pathway and apoptotic pathway), p53 signaling (p53 class mediator), DNA damage response, NF-κB signaling, and toll-like receptor signaling are dysregulated after TBI and WTLI and may be affected by changes in the circulating miRNA profile. Given that it is well established that radiation injury affects these pathways, these data provide an indication that miRNA may be a powerful class of biomarkers for understanding systemic dysregulation.
A consequence of our approach, however, is that it may overestimate the number of genes impacted by changes in the circulating miRNA profile. Since the functional outcomes of miRNA remains poorly understood, our GOEA for miRNA involved mapping miRNA to genes using interaction databases. It is not clear that these interactions occur in vivo within the relevant tissues or if these circulating miRNAs can exert functional biological effects, such as significantly affecting the abundance of mRNA transcripts or translated proteins for genes of interest. Another potential concern is that this analysis did not take the direction of changes (up- versus down-regulation) into account. However, the analysis does provide bases for testable hypotheses for later validation, including mechanistic studies or characterization of other biomolecules (e.g., miRNA), and the data suggest that outcomes after TBI or WTLI may correlate with dysregulation of specific pathways.
The TBI miRNA panel specifically enriched GO terms associated with ubiquitination and ubiquitination of histone H2A (Fig. 7). Ubiquitination is a post-translational modification of proteins that signals protein transport machinery to send the tagged protein to the proteasome for degradation. Ubiquitination is also involved in endocytosis, autophagy, transcription, immunity, inflammation, and DNA damage repair. Histone H2A ubiquitination, in particular, is associated with DNA damage repair pathways (45). And, the role of ubiquitination in radiation response has been recognized for years (46). Histone H2A K63-linked ubiquitination (GO:0070535) was observed in both TBI and WTLI, but 19 GO terms specific to ubiquitination were enriched in the TBI cohort. Consistent with our findings, six of our eleven (55%) miRNAs identified in the TBI cohort (let-7a-5p, miR-16-2-5p, -93-5p, -122-5p, -223-3p and -375-3p) were identified in an independent study focused on identifying miRNA with a role in histone modification, including H2A. Only two of our eleven (18%) miRNAs from the WTLI cohort (miR-181a-5p and -200c-3p) were present (47). In addition to DNA damage repair, K48 ubiquitination (GO:0070936) regulates p53 homeostasis (48, 49).
WTLI resulted in more unique GO terms than TBI [314 vs. 153 (Fig. 5)]. This may be due to the different time points used in each study (D15 vs. D2); however, the D15 time point is not practical for TBI, since the onset of severe neutropenia occurs within 7 days after high-dose irradiations and earlier time points after WTLI did not correlate as well with survival. Many of the GO terms uniquely enriched in WTLI were associated with terms also enriched in TBI [e.g., TNF and interleukin pathways (Fig. 6)]. This suggests that these pathways are not selective for WTLI survival. TGF-β and BMP signaling were highly specific to WTLI and have a well-established mechanistic relationship with radiation-induced pulmonary fibrosis (50–54), as observed by histology at the study endpoint (25). miRNA linked to these pathways were also identified in mice after WTLI (7). Positive regulation of Wnt signaling pathway by BMP signaling pathway (GO:0060804) was the only BMP-specific GO term enriched by TBI in the BMP DAG (Fig 8). Similarly, transforming growth factor-beta secretion (GO:0038044) was the only TGF-specific GO term in the TGF DAG (Fig 9). Signal transduction (GO:0007165), a parent term in both the BMP and TGF DAGs, was enhanced in both TBI and WTLI, suggesting irradiation widely disrupts signaling pathways [p53, NF-κB, Wnt, MAPK, TLRs, SMAD, interferons, interleukins, etc. (Fig 6)]. Positive regulation of cytokine secretion (GO:0050715), a parent term for negative regulation of transforming growth factor-beta secretion (GO:2001202), was also identified by the TBI miRNA.
FIG. 9
GO DAG of enriched terms related to transforming growth factor. GO terms related to transforming growth factor activation, production, and secretion were enriched by the WTLI miRNA. GO terms are colored by P value with blue colors representing TBI and orange colors representing WTLI, as shown in the key.
This does not suggest that TGF-β or BMP signaling is not differentially regulated after TBI or that ubiquitination is not altered after WTLI. The miRNA used for GOEA were selected because of an abundance correlation with survival and have not been validated as biomarkers. Rather, these data suggest that these pathways (ubiquitination for TBI, and BMP/TGF-β signaling pathways for WTLI), in addition to the general radiation signature (apoptosis, p53 signaling, pro-inflammatory immune response, etc.) are related to the different radiation-model-specific injuries.
Based on our comparative study between two different ionizing radiation models in NHPs, we have shown that complex pathophysiology can be predicted and tracked with circulating miRNA. This work supports our thesis that depending on which part of the body is irradiated, the outcomes can be very different and cannot be predicted by dose alone. A tissue-specific predictive assay of injuries resulting from exposure to ionizing radiation after a RAD-NUC event could improve outcomes for individuals through the continuum of care.
CONCLUSION
Taken together, these data underscore the enormous potential for miRNA-based biomarker panels. Using a directed panel of only 11 markers, GOEA identified the key molecular contributors to radiation injury. Importantly, differences in miRNA signatures were used to identify unique pathways that may contribute to the onset of severe neutropenia after TBI and death due to late lung injury after WTLI. Further, if future studies are able to identify and validate a robust molecular signature associated with neutropenia, it may be possible to select patients likely to benefit from FDA-approved neutropenia mitigators for radiation exposure (55). In the case of WTLI, the mechanistic connection to the disease and the identified pathways is clear. BMP/TGF-β signaling has been shown to be an important contributor to radiation-induced fibrosis. While the role of ubiquitination in the onset of severe neutropenia, or H-ARS, is less understood, ubiquitination has been shown to be involved in DNA damage response, autophagy, transcription, immunity, and inflammation, which is consistent with our current understanding of H-ARS pathogenesis. It should be noted that this novel approach to radiation biodosimetry is based on predicting tissue-specific injury. We have shown the potential of our approach in two different radiation injury models by tying circulating miRNA biomarkers to fundamental genetic pathways impacted by ionizing radiation injury over a time scale ranging from days to months. We believe this fundamentally unique approach to understanding complex outcomes from systemic insults/injuries will open a new paradigm in future approaches to diagnostics. We anticipate that, with additional development, miRNA panels that predict tissue-specific injuries from exposure to ionizing radiation will play a critical role through the continuum of care that provides clinicians with actionable information to treat their patients.
ACKNOWLEDGMENTS
This project has been funded by the National Institute of Allergy and Infectious Disease (Contract Number HHSN272201700012C and Grant Number 2R44AI108019-04) and Biomedical Advanced Research and Development Authority (Contract Number HHSO100201000005C). ChromoLogic LLC staff adheres strictly to Ethical Scientific Conduct guidelines documented by the company in their Employee Agreement and all workers adhere to guidelines established by the National Institutes of Health. The authors would like to thank James Axtelle for administrative assistance; William McBride, Dorthe Schaue, Joel Greenberger, and Thomas Register for helpful discussions; and the staff at CiTox Laboratories (TBI) and Altascience (WTLI) for care of the animals.
