Adverse outcome pathways and linkages to transcriptomic effects relevant to ionizing radiation injury

Abstract

Background

In the past three decades, a large body of data on the effects of exposure to ionizing radiation and the ensuing changes in gene expression has been generated. These data have allowed for an understanding of molecular-level events and shown a level of consistency in response despite the vast formats and experimental procedures being used across institutions. However, clarity on how this information may inform strategies for health risk assessment needs to be explored. An approach to bridge this gap is the adverse outcome pathway (AOP) framework. AOPs represent an illustrative framework characterizing a stressor associated with a sequential set of causally linked key events (KEs) at different levels of biological organization, beginning with a molecular initiating event (MIE) and culminating in an adverse outcome (AO). Here, we demonstrate the interpretation of transcriptomic datasets in the context of the AOP framework within the field of ionizing radiation by using a lung cancer AOP (AOP 272: https://www.aopwiki.org/aops/272) as a case example.

Methods

Through the mining of the literature, radiation exposure-related transcriptomic studies in line with AOP 272 related to lung cancer, DNA damage response, and repair were identified. The differentially expressed genes within relevant studies were collated and subjected to the pathway and network analysis using Reactome and GeneMANIA platforms. Identified pathways were filtered (p < .001, ≥3 genes) and categorized based on relevance to KEs in the AOP. Gene connectivities were identified and further grouped by gene expression-informed associated events (AEs). Relevant quantitative dose-response data were used to inform the directionality in the expression of the genes in the network across AEs.

Results

Reactome analyses identified 7 high-level biological processes with multiple pathways and associated genes that mapped to potential KEs in AOP 272. The gene connectivities were further represented as a network of AEs with associated expression profiles that highlighted patterns of gene expression levels.

Conclusions

This study demonstrates the application of transcriptomics data in AOP development and provides information on potential data gaps. Although the approach is new and anticipated to evolve, it shows promise for improving the understanding of underlying mechanisms of disease progression with a long-term vision to be predictive of adverse outcomes.

Keywords:

• Gene expression
• radiation
• risk assessment
• AOPs
• adverse outcome pathways

Introduction

Adverse outcome pathways (AOPs) represent a series of sequential events that collectively culminate in a negative health outcome. Structured AOPs provide a visual landscape of the current understanding of disease progression at increasing levels of biological organization observed at the macromolecular, cell, organ, tissue, organism, and/or population level (Ankley et al. 2010). This depiction of biology begins with a molecular initiating event (MIE) induced from an external stressor(s) (i.e. chemical, ionizing radiation, etc.) followed by induction of a set/series of key events (KEs) linked by key event relationships (KERs) to an adverse outcome (AO) of interest to support regulatory decision-making (OECD 2016).

A fundamental aspect surrounding AOPs is that they are developed and supported using a weight of evidence approach that considers the vast amount of data/literature surrounding the components of the AOP. Through various literature mining techniques, these linkages are identified using evidence streams defined by the modified Bradford Hill criteria (Becker et al. 2015). The benefit of an AOP is that, as a living document, it allows areas that lack sufficient evidence to be further supported when additional evidence becomes available. Within the last few years, high-throughput and high-content technologies such as genomics have received increased attention in elucidating the mechanism(s) of disease, with the added possibility of supporting predictive toxicology (Rosenstein 2017).

Data generated from rapidly advancing new investigative methodologies are valuable for supporting AOPs by providing additional data to strengthen causal connectivities across the AOP continuum. Furthermore, mechanistic incorporation of high content “transcriptional” data may help to elucidate gene signature profiles and associated pathways underlying/linking the KEs to the AO. These signature pathways can be thought of as causal subnetworks denoting each KE and the AOP framework can be annotated with this information to enrich it in terms of mechanistic information. Integration of transcriptional data into AOPs potentially allows for the delineation of novel KEs and/or mechanisms/systems where the initial weight of evidence to support their entry into an AOP framework was lacking in traditional non-genomic literature. Indeed, in the chemical toxicity field, work is in progress on considerations for best practices for use of transcriptional data in risk assessment. Currently, there is recognition that there is a lack of available studies to demonstrate linkage to apical endpoints, and agreement that transcriptional data offer the potential for informing health risk assessment when complemented with a systems biology approach (Cassman 2005). Nonetheless, case examples need to be developed to help move efforts forward in this direction (Brockmeier et al. 2017).

Recently, our group has demonstrated the potential application of the AOP framework within the field of ionizing radiation (IR), illustrating a working case example of lung cancer (AOP 272: https://www.aopwiki.org/aops/272) (Chauhan et al. 2021). Lung cancer is relevant to a variety of environmental and occupational exposures, particularly radon gas, the second leading cause of lung cancer, next to smoking (Darby et al. 2005). Previous studies have shown that an early biological response to stressors, such as radiation exposure, is the initiation of key genes in DNA repair pathways (Hoeijmakers 2001; Jackson and Bartek 2009; Jeggo et al. 2016). This AOP was developed using data from decades of research in the area of DNA damage/response and repair. It begins with deposition of energy (considered as an MIE), to DNA lesions (i.e. base damage, single- or double-strand breaks) (KE1), followed by quick induction of DNA repair mechanisms which in some cases can be inaccurate (KE2). This can result in chromosomal aberrations (KE4) or sequence alterations (i.e. mutations) (KE3), which can have detrimental consequences, especially when found in critical genomic regions (i.e. those containing proto-oncogenes, tumor suppressor genes, and/or DNA damage response genes). The increase in mutations in critical regions and/or chromosomal aberrations can lead to genomic instability and proliferative cellular events (KE5) that eventually result in lung cancer (the AO). These KEs could be viewed as apical endpoints and could be correlated to transcriptional changes as they are the observable outcomes indicating the disease state. The AOP could be a working example for the integration of gene-based data, particularly at the transcriptional level, as genes involved in multiple DNA repair pathways have shown promise to be used as predictive biomarkers relevant for assessing environmental radiation exposure and the associated health risk (Zhao et al. 2020; Costa 2022).

The feasibility of incorporating available transcriptomics data in the area of DNA damage sensing and repair into the AOP analytical construct would help delineate a visual representation of gene profiles and networks underlying KEs. In particular, it would further guide current molecular-level understanding and provide direction for the interpretation of transcriptomics in AOP development. Therefore, the objective of the present study was to incorporate transcriptomics data as they relate to AOP 272. Under this premise, a workflow was built using data mining and informatics tools (e.g. Reactome and GeneMANIA), to identify gene signature profiles and pathways (e.g. subnetworks) underlying/linking KEs, thereby providing a further mechanistic understanding of lung cancer as it relates to IR exposure. The results from this study contribute to the knowledge of the benefits and challenges of developing transcriptomics-informed AOPs.

Materials and methods

Data retrieval tools

DistillerSR platform (Evidence Partners: https://www.evidencepartners.com/) was used to facilitate screening and data extraction and to identify relevant transcriptomic-based studies that support either new or existing KEs in the proposed AOP pathway to lung cancer (AOP 272) that has been deposited in the AOP-Wiki (https://www.aopwiki.org/aops/272) (DistillerSR 2021). This AOP consists of five KEs including the formation of DNA double-strand breaks, inadequate repair, mutations, chromosomal aberrations, cell proliferation, and seven adjacent KERs. Reactome, a widely cited and accepted database was applied for pathway analysis (Fabregat et al. 2017). GeneMANIA was used alongside Cytoscape v3.4.0 to visualize the genetic and physical network to explore the correlation and interaction among radiation-responsive genes (Shannon et al. 2003; Warde-Farley et al. 2010).

Literature search strategy and criteria

A literature screening of relevant radiation studies in the context of transcriptional changes was undertaken using scoping review methodology. Published peer-reviewed literature was identified using two electronic databases: MEDLINE and Embase. Literature search strategies and terms can be found in Supplementary Document A. The searches were not limited by language or publication date. Studies underwent two rounds of screening. The level-1 screening included the assessment of the title and the abstract of the literature using PECO(E) elements (Population, Exposure, Comparator, Endpoints, and Outcomes) (Table 1). For studies that met the eligibility criteria, a full-text review of PECO(E) elements was performed at level-2; studies that did not meet the PECO(E) criteria were excluded (Supplementary Table 1). Additionally, studies of all types (in vitro, in silico, in vivo, epidemiological, cohort), all levels of biological organization (molecular, cellular, tissue, organ, individual, population), and transcriptional studies that contain pathway connectivity and networks were included. Abstracts-only publications, conference abstracts, posters, theses/dissertations, and presentations were not included. Non-English studies were considered, provided the essential data could be easily identified within the abstract description. The level-1 screening (title and abstract) of the search results was conducted independently by two reviewers to include relevant studies. Conflicts were resolved by an independent reviewer.

Table 1. PECO(E) elements.

Element	Evidence
Populations	Animal, human, ex vivo, in vivo, cell, cell-lines, epidemiological, cohort, population, organ, tissue, cellular, molecular, cellular components.
Exposures	Radiation, radon, alpha particles, high linear energy transfer
Comparators	No exposure vs exposure
Endpoints	Genomic changes, gene expression, transcriptional responses
Outcome	Lung cancer, DNA damage, DNA repair, cell proliferation, mutations, chromosomal aberrations

Note: The eligibility of the study was based on the population, exposure, comparator, endpoint, and outcome statements provided in this table. This was applied for both abstract level and full-length context level screening to collect relevant references. To be included in the scoping review, studies must demonstrate primary research investigating the link at each element, but not all elements need to be fulfilled.

Data extraction

At level 2-full-text review, the following information was extracted from the included studies for further analysis: reference information, title, study type including in vivo (animal model and sex) and/or in vitro (tissue, cell type, and morphology), tissue source, platform, gene type, conditions (radiation type, gene editing, and treatment), dose, dose concentration, dose rate, time, fraction information, and fold change or expression level. Fold change data points were manually collected from figures by the Matlab R2017b (version 9.3.0.713579) with package GRABIT (version 1.0.0.1) (Jiro 2021) and WebPlotDigitizer (version 4.5) (Rohatgi 2021). If available, raw transcriptional response data were downloaded from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). All data for this review were exported to Microsoft Excel and subsequently analyzed by R (version 4.1.1) (R-Core-Team 2021). Data collected from different assay platforms, including microRNA, mitochondrial DNA, and single-nucleotide polymorphism were further filtered to identify those studies examining mRNA only, and within these all data units were normalized by converting them to a standard radiation unit (Gy), dose rate (Gy/min), and incubation time (hours). Studies that did not provide adequate information relating to specification on stressor type, dose range, tissue, time-points, model, and experimental procedures used were not considered.

Data analysis and visualization

Pathway analysis

In the data extraction step, differentially expressed genes reported in selected references, derived through appropriate statistical methods (e.g. correction for false discovery rate) with information on fold change or gene expression between experimental conditions were collected manually. After combining the genes from each study and removing repetitive genes, a name list of 988 genes in total from 43 eligible studies based on PECO criteria, was exported and tabulated into a master file (Supplementary Table 2). HUGO (Human Genome Organization) symbols were searched for each gene to reconcile the identification of the same genes with different names across studies that used different model systems (e.g. rat, mouse, human) (Tweedie et al. 2021). From the master file, studies were further narrowed to those that used mRNA microarray platform on human-derived cells (12 studies in total, corresponding to 404 genes, Supplementary Table 3) and analyzed in Reactome Knowledgebase (https://reactome.org) (version 79, February 3, 2022) (Griss et al. 2020; Jassal et al. 2020; Gillespie et al. 2022). The HUGO symbols were applied as input for Reactome pathway analysis. With cross-reference with other resources such as NCBI, Ensembl, UniProt, KEGG, ChEBI, PubMed, and GO, Reactome provides the pathways mapped with submitted genes. Results from the analysis were exported and pathways with a p-value less than .001 and genes ≥3 were queried against AOP 272 to identify matching/similar KEs or potential new ones. The analysis yielded many potential pathways that could inform KEs, among which the top 92-mapped pathways (p-value < .001 & mapping genes ≥3) were selected. There are in total 29 high-level biological processes presented as locations in the PathwayBrowser on Reactome analysis, among which seven were relevant to the 92 mapped pathways. The associated events (AEs) were suggested/identified by the 7 high-level locations of the mapped pathways, including the immune system response, cell cycle, TP53-mediated transcriptional regulation, DNA replication, DNA repair, cellular responses to stimuli, and disease.

Gene-based network development

Among the 404 genes uploaded for pathway analysis, 228 genes were mapped to the seven potential AEs mentioned above (Supplementary Table 4). All 228 genes were analyzed in GeneMANIA, an online database for biological network integration for gene prioritization and prediction of gene function (Warde-Farley et al. 2010). By removing genes that did not have any connection to the main network, the gene network was still too large to visualize. Therefore, we removed genes that had little weighting factor of edges for better visualization (e.g. genes with the score for the interaction > 0.019 in GeneMANIA were kept). Thus, the list was reduced to 61 genes from high-dose study and 24 genes from moderate dose studies that were still associated with the 7 AEs in the main network. The strength of the relationships was depicted visually using a pictorial diagram. Arrows drawn on edges show the direction of interaction from the source gene to the target gene, and the thickness of the edges indicated the strength of the interaction. The nodes represent the composition of AEs. To increase the readability of the network, it was further customized from GeneMANIA into Cytoscape (version 3.9.0) (Shannon et al. 2003), which provided more features by merging a pie chart of AEs with corresponding genes as well as creating compound groups of genes sharing the same AE.

Quantitative data analysis

Not all of the 43 studies had data that could be quantified, with only 8 studies having GEO datasets available that could be considered for quantitative analysis in this study. Five studies were excluded for quantitative analysis as the datasets were not appropriate for this study, such as SNP data, drosophila study, or extremely high radiation dose (90 Gy). The remaining 3 studies (in vivo mice model or in vitro human-derived cell model, high and moderate radiation doses) were selected, including Jackson et al. (2016) (GSE85359), Ghandhi et al. (2011) (GSE21059), and Ding et al. (2013) (GSE44282).

Jackson et al. applied two murine models (“radiation sensitive” C57L/J mice and “radiation resistant” normal background strain C57BL/6J mice) to compare gene expression in lungs at 24 hours after a single dose of 0 Gy, 12.5 Gy, or 15 Gy of whole thorax lung x-ray irradiation (dose rate 67 cGy/min) (Jackson et al. 2016). Ghandhi et al. used IMR-90 human lung fibroblast cell cultures that were irradiated with alpha particles at 0 Gy or 0.5 Gy and incubated for 30 minutes, 1, 2, 4, 6, and 24 hours after exposure, to characterize changes in gene expression (Ghandhi et al. 2011). In addition, Ding et al. performed transcriptome profiling in human bronchial epithelial cells (HEBC3KT) exposed to cesium-137 gamma rays to explore the molecular response to dose (0, 1, 3 Gy) and time post-irradiation (1, 4, 12, 24 hours) (Ding et al. 2013).

The datasets were grouped based on radiation exposure levels. Ghandhi et al. (2011), and Ding et al. (2013) were grouped as moderate doses, leaving Jackson et al. (2016) as the sole high-dose dataset. The data from Ghandhi et al. (2011) and Ding et al. (2013) underwent cross-platform normalization (XPN) to allow for direct comparison (Shabalin et al. 2008). XPN was performed using the R package CONOR (version 1.0.2) (Rudy and Valafar 2011). The data from Jackson et al. (2016) required a log transformation but did not undergo any further normalization. The fold change of genes for both high and moderate dose data was then found using differential expression analysis with R package limma (version 3.50.0) (Ritchie et al. 2015) that utilizes a moderated t-test-based approach. Multiple testing correction was done using a false discovery rate with a threshold of 0.05.

All moderate dose studies were performed on human cells, however, the high dose study utilized Mus musculus. Therefore, the expression levels of the orthologous genes were used for this study. These orthologs were found manually using GeneCards (https://www.genecards.org/) and MGI (http://www.informatics.jax.org/).

Results

Study composition

The literature search resulted in initial 700 studies, consisting of 538 from Pubmed/Medline, 155 from NERAC service of Embase and International Pharmaceutical Abstracts (IPA), and 7 additional recommended studies. After duplicate removal, abstract and full-text review, only 43 studies met the PECO(E) statement and selection criteria for data collection, among which 12 studies were the basis of further pathway analysis (Table 2).

Table 2. Summary of selected papers for pathway analysis.

Reference	Radiation type	Model type	Cell/Organ type	Dose (Gy)	Dose rate (Gy/min)	Assay platform	Exposure time	Post-irradiation incubation time
(Liu, et al. 2019)	Alpha particle	In vitro	BEP2D, BERP35T1, BERP35T4	1.5, 0.0	NA	Microarray, qRT-PCR	NA	NA
(Ding et al. 2015)	^28Si, ^56Fe, gamma-ray	In vitro	HBEC3KT, HBEC-P53KRAS	NA	NA	Microarray	NA	1, 4, 12, 24 hours
(Chauhan and Howland 2014)	Alpha particle, X-ray	In vitro	HFL-1	0.5, 1.0, 1.5	0.016	Microarray, qRT-PCR	0, 0.5, 1, 1.5 min	24 hours
(Chauhan et al. 2012)	Alpha particle	In vitro	A549	0.90, 0.30, 0.03	0.016	Microarray, qRT-PCR	0, 0.03, 0.3, 0.9 hours	4, 24 hours
(Sokolov et al. 2006)	Gamma ray, ^125I-IUdR	In vitro	IMR-90	2.0000, 0.0016	1e + 00, 8e-04	Microarray, qRT-PCR	2 min	2, 24 hours
(Guo et al. 2005)	Gamma ray	In vitro	A549, NCI-H446	5, 0	NA	Microarray	0, 3.125 min	6, 24 hours
(Khodarev et al. 2001)	X-ray	In vivo & and in vitro	HEL, U87, and U87 xenografts from athymic nude mice	3, 1, 10	1.18	Microarray	0.85, 2.54, 8.47 min	5 hours
(Xie et al. 2011)	X-ray	In vitro	SPCA1, A549	0, 80	0, 3	qRT-PCR, microarray	26.667 min	4 weeks
(Ghandhi et al. 2011)	Alpha particle	In vitro	IMR-90	0, 0.5	NA	Microarray, qRT-PCR	NA	0.5, 1, 2, 4, 6, 24 hours
(Ghandhi et al. 2008)	Alpha particle	In vitro	IMR-90	0, 0.5	NA	Microarray, qRT-PCR	NA	4 hours
(Ghandhi et al. 2014)	Alpha particles	In vitro	mouse ES, H1299	1	NA	Microarray, qPCR	NA	4 hours
(Ding et al. 2013)	^28Si, ^56Fe; gamma ray	In vitro	HBEC3KT	0, 0.5, 1; 0, 1, 3	NA	Microarray, qRT-PCR	NA	0, 1, 4, 12 and 24 hours

The selected 43 references comprised exposures across different radiation types including x-rays, gamma rays, alpha particles produced by radon and other emitters, and energetic heavy ions (Figure 1A). Most studies were focused on x- or gamma rays (13 and 11, respectively), followed by alpha particles from radon gas and other emitters such as americium-241 (²⁴¹Am) or plutonium-238 (²³⁸Pu) (7 and 11, respectively). Studies directly relevant to energetic heavy ions exposure were minimal, including three studies on iron, two studies on silicon, one study on space radiation, and one on iodine. Regarding the model types, around 1/3 of the collected data were from in vivo studies, including mice, rats, and fruit flies; 2/3 of the studies were in vitro conducted mostly on human-derived cell models, either immortalized or cancer cell lines or primary cells harvested from patients (Figure 1B and C). Transcriptomic data were derived using various experimental platforms, but mainly qRT-PCR (64.3%) and microarray (26.8%) (Figure 1D).

Figure 1. Summary of studies. The number of studies retrieved across (A) different radiation stressors (B) Model type (C) Taxonomy (D) Assay Platform. HiCEP stands for high coverage gene expression profiling.

Pathway analysis

As the datasets were quite diverse, it was difficult to derive meaningful pathway results, therefore the studies were prioritized based on the model type and the platform used for transcriptomic analysis. Pathway analysis was conducted using only human data/studies derived from in vitro data, followed by selecting the information from microarray platforms. A total of twelve studies met this eligibility criterion. Differentially expressed genes from these studies were collected leading to a final gene list of 404 with available HUGO symbols. These studies were subjected to Reactome analysis, in which 92 biological pathways were hit by at least three of the genes with p-value < .001, and 228 mRNAs were involved in these significant pathways (Supplementary Table 5).

A stable identifier is assigned to each mapped biological pathway in the Reactome analysis result, from which a hierarchy of locations in the Pathway Browser could be found. The top of the locations is identified as the high-level biological processes, including seven locations covering the 92 biological pathways. Therefore, the seven high-level biological processes of the significant biological pathways with corresponding genes were classified as potential AEs to the proposed lung cancer AOP 272 (Figure 2). The AEs were categorized as follows: twenty-three biological pathways are involved in cell cycle processes, such as mitotic G1 phase, G1/S transition, G2/M transition, etc. Seventeen biological pathways are involved in immune system regulation, including antigen presentation, endosomal/vacuolar pathway, ER-phagosome pathway, etc. Five pathways were associated with DNA replication and DNA repair, respectively. Additionally, oxidative stress-induced senescence was identified in the high-level biological process of cellular responses to stimuli.

Figure 2. Pathway analysis. Pathway analysis of 404 genes was performed using REACTOME and correlated to AOP 272 with respect to changes involving molecular, cellular, and tissue level events. In brackets, the number of biological pathways identified using p-value < .001 and greater than 3 genes is provided. The orange represents high-level biological processes and the number of genes involved in the pathway. Lower level categories of pathways are provided as examples, which were selected by low p-values. The figure is created with Biorender.com.

Genes involved in these significant pathways were further analyzed by GeneMANIA to develop network connectivities as AEs in the context of the KEs in AOP272. The relationships between each gene and AEs are illustrated in Figure 3 for moderate-dose and Figure 4 for high-dose. It was noted that some genes were shared across multiple AEs. For example, CDK1 (cyclin-dependent kinase 1) was shared by seven high-level biological processes, TP53 (tumor protein p53) and PCNA (proliferating cell nuclear antigen) were shared by six and four AEs, respectively, while CDC20 (cell division cycle 20) and KIF20A (kinesin family member 20 A) were shared by two AEs (Figure 3A and 4A). Transcripts with strong weighted connections such as TP53 to COP1 (COP1 E3 ubiquitin ligase), CYFIP2 (cytoplasmic FMR1 interacting protein 2) to TK1 (thymidine kinase 1), and CDK1 to CCNB2 (cyclin B2) may be suggestive of a high correlation in connectivity. TP53, CCNA1 (cyclin A1), and CDK1 were included in both high- and moderate-dose gene networks, as well as genes (e.g. CDC6 (cell division cycle 6), GINS2 (GINS complex subunit 2), etc.) involved in the cell cycle and DNA replication process were found at both dose groups (Figure 3A and 4A). For example, LIF (LIF interleukin 6 family cytokine) identified only in the immune system response group was found in the network at moderate dose but not in the high-dose network. SOD1 (superoxide dismutase 1) was involved in both immune response and cellular response to stimuli. It was noted that the moderate-dose network had predominately genes associated with pathways related to cell cycle and DNA replication while the higher doses had in addition to these, relatively more pathways related to the immune system and disease processes.

Figure 3. Gene network and quantitative analysis of moderate-dose radiation study. (A) 24 genes were subjected to network analysis for the moderate dose groups. The color of nodes represents the pathway associated with each gene. The strength of the connection was shown by the thickness of the edges. Genes were grouped with similar colors together in order to visualize AEs. (B) The studies from Ding et al. (2013) (GSE44282) and Ghandhi et al. (2011) (GSE21059) were applied for quantitative analysis. The graph shows the intersection of each gene among the 7 AEs. The dots correspond to AE hit by genes: AE1 cell cycle; AE2 DNA replication; AE3 DNA repair; AE4 TP53-mediated transcriptional regulation; AE5 cellular responses to stimuli; AE6 immune system responses; AE7 disease (C) shows the fold change level of genes at moderate doses (0.5Gy alpha particle and 1Gy gamma ray) at 24 hours post-irradiation.

Figure 4. Gene network and quantitative analysis of high-dose radiation study. (A) 61 genes were subjected to network analysis for the high doses groups. The color of nodes represents the pathway associated with each gene. The strength of the connection was shown by the thickness of the edges. Genes were grouped with similar colors together in order to visualize AEs. (B) The study from Jackson et al. (2016) (GSE85359) was applied for quantitative analysis. The graph shows the intersection of each gene among the 7 AEs. The dots correspond to AE hit by genes: AE1 cell cycle; AE2 DNA replication; AE3 DNA repair; AE4 TP53-mediated transcriptional regulation; AE5 cellular responses to stimuli; AE6 immune system responses; AE7 disease. (C) Shows the fold change level of genes at high doses (12.5Gy and 15Gy) of radiation at 24 hours post-irradiation.

Quantitative analysis

To determine the causal connectivity of genes across the AEs, a gene map was generated that listed the genes and the AEs with the quantitative expression levels. Of the 43 selected studies, only three had dose-response data (mRNA microarray data available in the GEO datasets) that could be used to extract the pattern of quantitative response. The studies that were eligibly comprised one using mice lungs exposed to a high dose of radiation (12.5 Gy and 15 Gy of x-rays) and two studies that used human lung cell cultures (apparently normal fibroblast cell line IMR-90 and immortalized bronchus epithelial cell line HEBC3-KT, respectively) exposed to moderate doses (0.5 Gy alpha particles and 1 Gy gamma-ray). Several genes showed a similar pattern in fold change, such as MCM4 (minichromosome maintenance complex component 4), MCM2 (minichromosome maintenance complex component 2), and CDC6. These genes increased expression at both 0.5 Gy of alpha particles and 1 Gy of gamma-ray and were associated with both cell cycle and DNA replication (Figures 3B and C). Additionally, in the high dose group, MCM2, MCM4, MCM6 (minichromosome maintenance complex component 6), and CDC45L (cell division cycle 45) had a similar pattern of fold change (reduced at both 12.5 Gy and 15 Gy of x-ray) and were associated with both AEs hit by the above genes at moderate doses (Figures 4B and C). Similar patterns of connectivity were not observed for many of the other genes.

Discussion

The current investigation explored the application of gene signatures and pathway ontologies in the context of AOPs and radiation-induced injury at the macromolecular to organ level through the integration of transcriptomics data using an existing AOP as a case example (AOP 272). This guided the development of a gene network that supported underlying KEs across the AOP. To achieve our objective, we utilized the wealth of data surrounding genomic-based studies in the area of DNA damage response/repair and lung cancer. Studies in this area are plentiful, with the past few decades generating a plethora of data that were sometimes translated to the functional level as they related to the adverse outcome (Sekido et al. 2003; Weir et al. 2007; Brockmeier et al. 2017; Rosenstein 2017). Although radon exposure is the more well-known cause of lung cancer development from environmental radiation, we also included in the analysis, transcriptomic studies generated from exposures to x-and gamma-rays, and energetic heavy ions, known to induce lung damage (Gaskin et al. 2018; Kennedy et al. 2018).

The list of differentially expressed genes from the identified studies was subjected to pathway analysis and mapped to AEs that may connect to and enrich the AOP 272 across certain KEs (i.e. DNA double-strand breaks, DNA repair) and also identify potentially new KEs (including cellular response to stimuli and immune system response) and AEs (i.e. DNA replication, cell cycle) that may be relevant to lung cancer. It is important to highlight that we recognize that certain genes may be connected with several pathways. Therefore, with this meta-analysis, it was more of interest to provide the landscape of pathways and gene responses that may elicit the KEs in AOP 272 to the accepted understanding of lung cancer development while highlighting those genes that appeared across multiple AEs.

Among the genes that were identified to be differentially expressed across different datasets, 61 and 24 genes available in GEO datasets were used for the development of the network for high- and moderate doses, respectively, using GeneMANIA (Supplementary Table 6). Among these genes, we note that some common genes were found within several AEs. One example is the TP53 gene (encoding P53 protein, a tumor suppressor) that is associated with six AEs, including cell cycle, DNA repair, immune system responses, TP53-mediated transcriptional regulation, cellular responses to stimuli, and disease. We also note that some genes had bi-directional connectivity using GeneMANIA. For example, the proteins encoded by the identified genes involved in DNA replication and cell cycle can regulate each other, such as MCM4 regulates CDC45L, and CDC6 regulates MCM2 (Figures 3 and 4). These specific genes and their corresponding proteins may act as modulators to AEs and are therefore valuable in developing strategies for preventing disease progression.

In developing AOP 272, the focus was on KEs relevant to DNA damage and response, as this is well-studied and also relevant to stressors such as radon gas that induces clustered damage leading to double-strand breaks (Chauhan et al. 2021). However, we recognize that oxidative stress induced by IR may also serve as an important KE that was not included in AOP 272. Interestingly, our transcriptional analysis did identify expressed genes related to oxidative stress (e.g. SOD1, SOD2, HIF1A), highlighting these to be important underlying events. The oxidative stress response, which involves the activation of antioxidant enzymes and scavengers of free radicals and other oxidizing species, has been demonstrated to be an important response to radiation injury and also correlated to disease (Liou and Storz 2010; Cecchini and Cecchini 2020; Saleh et al. 2020; Shenoy 2020; Alam and Czajkowsky 2022; Kalyanaraman 2021). It was also noted that among the AEs identified by pathway analysis and applied to the gene-informed AOP diagram (Figure 5), the immune system response is mapped with seventeen biological pathways in the Reactome database. Specifically, 9 out of 10 top pathways are located in the immune system. These pathways are related to fundamental immune processes including interferon and other cytokine signaling, antigen presentation and processing, phagocytosis, and internal events initiated (i.e. endosomal/vacuolar and ER-phagosome), amongst others. Although the immune system response is technically classified as organ/system level of organization in current AOP under development, we recognize that it works ubiquitously in the human body environment.

Figure 5. Gene-informed AOP diagram. A schematic summary of identified associated events (AEs) and select high expressing genes identified using Reactome from radiation-induced transcriptomic studies. The AEs (light green) are grouped to support KEs (dark green) in AOP 272 (https://www.aopwiki.org/aops/272). Some AEs (e.g. cellular response to stimuli and immune system response) are identified as new KEs (red) not represented in AOP 272.

Many of the gene products grouped in the immune system not only interact with each other but also express redundancy with genes located in other AEs, such as TP53-mediated transcriptional regulation and cell cycle. Our results indicate that the immune system response plays a fundamental role in the AOP diagram as a silent but omnipresent player. For instance, potential candidates identified include the gene CDKN1A (cyclin-dependent kinase inhibitor 1 A) implicated in the immune system and DNA damage response and is highly expressed following exposure to high and low doses of IR, and the multipurpose/pleiotropic genes GINS2 and KIF20A (Kim et al. 2016). However, more thorough and systematic research or new experiments are necessary to prove any causal relationship between transcriptional responses and connectivity to apical endpoints to KEs in the AOP.

Several challenges were encountered in collecting genomics studies to support AOP 272. Although we successfully collected data from different platforms, such as sequencing, qRT-PCR, and microarray investigations, the analysis of quantitative data across these platforms was met with great difficulty. Aside from experimental assays themselves, there was no consistency in methodological approaches concerning IR type, fractionation, exposure time, and incubation (recovery) time, leading to acute or chronic alterations in gene expression. Although various types of irradiation may induce the same AO, the specific level of detriment may be potentially different at the quantitative level. With limited studies, the quantitative analysis was conducted on select datasets that used multiple doses, which were limited to those using alpha particle sources simulating the effects of radon (Table 2). Although ²³⁸Pu and ²⁴¹Am are widely applied to mimic natural radon gas in an experimental environment, the real situation is not fully recaptured which is a limitation. Nonetheless, quantifying cell-level data is informative and through future studies, this information could be correlated to tissue and organ level responses. In addition to radiation types, doses collected to induce lung damage ranged from low and medium (0.01–5 Gy), high (15 Gy), to extremely high (35 Gy, 90 Gy), and the dose rate of exposure was not guaranteed equivalency, which limited comparison across studies. Thus, we focused on identifying studies that used the same assay platform (microarray) to help reduce external variation; however, under this umbrella, the data resource availability represented an additional challenge in data collection. In this regard, among the 43 papers selected with valuable transcriptomic profiles in response to irradiation, only three datasets from GEO could be utilized to analyze the quantitative data. These datasets include relatively high doses of radiation (12.5 Gy and 15 Gy) (Jackson et al. 2016) and moderate doses (0.5 Gy and 1 Gy) (Ghandhi et al. 2011; Ding et al. 2013) only, and the quantitative analysis for low-dose radiation is not available. Furthermore, the dose range, number and spacing are not ideal for meaningful dose-responsive analysis, therefore this was not undertaken. However, fold change data across the genes in response to two doses from moderate and high dose irradiation (0.5 Gy vs 1 Gy, and 12.5 Gy vs 15 Gy, respectively) is compared to identify the quantitative responses. With the limited data, we did observe some causal connectivity across genes, however, more relevant data are needed to draw valid conclusions. In particular, there is a need for best practices for standardized experimental approaches for generating, storing, processing, and interpreting omics data that allow researchers to accurately connect omics perturbations to phenotypic alterations. This will increase confidence in data integration for regulatory risk assessment. Moreover, with high-quality continuous dose-response data, thresholds of dose activation can be derived using applications such as the DRomics tool and BMDExpress similar to the ecological and chemical toxicology fields (Crump 2002; Larras et al. 2018; Phillips et al. 2019; Mezencev and Auerbach 2020). This will be particularly relevant for the use of hematopoietically humanized mouse models as gene expression patterns can differ between ex vivo irradiated human blood cells compared with in vivo irradiated mice (Ghandhi et al. 2019).

The challenge of deriving causal relationships is not only in the context of transcriptional responses but also at the protein and metabolite levels. The overbearing complexity and the compensatory physiological mechanisms hinder the identification of the most relevant interactions that contribute to the pathophysiology (Yang 2020). For this reason, models need to be developed that are designed to study the interconnectivity across biological levels of the organization. An approach that could be explored is the use of Multiorgan microphysiological systems (MOMPSs). These technologies mimic human physiology and can be applied in the context of AOP studies. MOMPs allows for the continuous monitoring of alterations induced by stressors at the transcriptomic, proteomic, and metabolomics levels. The information derived from the whole process informs interactive networks and causal relationships for computational prediction by machine learning (Trapecar 2022). Although experiments with MOMPs are still at the proof-of-principle stage, the future premise to apply these technologies as multicellular tissue models would provide systematic research with ex vivo options that reduce dependence on in vivo animal models (Trapecar 2022).

Herein we present, how transcriptional responses following radiation exposure can be correlated to AEs to support KEs in an AOP, using a case example of lung cancer from stressors such as radon gas. The overall approach identified some interesting underlying genes/pathways that could be correlated to AEs, however, many challenges were encountered. There is a need for stronger evidence in the form of dose-response concordance data that can provide information on causal connectivity. In identifying relevant studies, it was evident that better coordination in the generation of genomics studies across institutions would help extract more meaningful data. The best use of omics data may be to improve the underlying mechanistic understanding of KEs in an AOP. This information could also help fill evidence gaps for incomplete KERs that lack sufficient evidence. Moreover, one of the advantages of doing this additional mechanistic layering is that the AOP will be refined in terms of insult-specificity (dose-related, dose-rate-specific effects based on the type of interactions). The work also highlights that in order for gene data to be incorporated into AOPs, there is a need to develop a coordinated framework, supported by genomic, transcriptomic, and proteomic studies. The knowledge gaps identified through this study warrant further research to advance AOP models for successful incorporation in predictive toxicology and help provide the causal connectivity to disease and thus support risk assessment.

Acknowledgments

The authors would like to acknowledge for insightful comments and edits to the manuscript. Special thanks to Dr. Edouard I. Azzam and Dr. Ruth C. Wilkins for their review, suggestions, and encouragement.

Disclosure statement

The authors declare they have no competing interests. The funders had no role in the design of the study, in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Additional information

Funding

This work was supported by Atomic Energy of Canada Limited’s Federal Nuclear Science & Technology Program, Genomics Research and Development Initiative, NSERC Discovery Grant, and Canada Research Chair Program.

Notes on contributors

Jihang Yu

Jihang Yu is a biologist at Canadian Nuclear Laboratories.

Wangshu Tu

Wangshu Tu Ph.D. is a post-doctoral fellow at Carleton University.

Andrea Payne

Andrea Payne is a 4th-year honor student at Carleton University.

Chris Rudyk

Chris Rudyk Ph.D. is a course instructor/lecturer at Carleton University and a scientific evaluator at Health Canada.

Sarita Cuadros Sanchez

Sarita Cuadros Sanchez is a research assistant at Health Canada.

Saadia Khilji

Saadia Khilji Ph.D. is a postdoctoral researcher at Health Canada.

Premkumari Kumarathasan

Premkumari Kumarathasan Ph.D. is a research scientist at Health Canada.

Sanjeena Subedi

Sanjeena Subedi Ph.D. is an assistant professor in the School of Mathematics and Statistics and a Canada Research Chair in Data Science and Analytics at Carleton University. Her research focuses on clustering and classification of high-dimensional data with application in bioinformatics.

Brittany Haley

Brittany Haley is a library assistant IV at Canadian Nuclear Laboratories.

Alicia Wong

Alicia Wong is a 4th-year honor student at McMaster University.

Catalina Anghel

Catalina Anghel Ph.D. is a computational research scientist at Canadian Nuclear Laboratories.

Yi Wang

Yi Wang Ph.D. is a research scientist and biologist at Canadian Nuclear Laboratories and an adjunct professor at the University of Ottawa.

Vinita Chauhan

Vinita Chauhan Ph.D. is a Senior Research Scientist at the Consumer and Clinical Radiation Protection Bureau of Health Canada. She is a Canadian delegate of the High-level group on low-dose research (HLG-LDR) and Extended Advisory Group on Molecular Screening and Toxicogenomics (EAGMST) of the OECD. She co-chairs the HLG-LDR Rad/Chem AOP Joint Topical Group and is the co-founder of the Canadian Organization of Health Effects from Radiation Exposure (COHERE) initiative.