The present study was carried out for understanding the impact of ionizing radiation on gene expression (mRNA diversity) and metabolic profile in mice blood after 1 Gy and 7.5 Gy radiation exposure at 24 hours. The workflow used in the present study is represented in Fig. 1.
Transcriptomic profiling
RNA sequencing was performed on RNA samples that passed the quality control (QC) indices. QC status and mean quality value across each base position can be observed from Supplementary Table 1 and Supplementary Fig. 1 respectively. The raw reads were mapped to 52636 mouse Ensembl genes IDs. However, ~23500 uniquely mapped genes were identified post-alignment to the mus musculus genome (9GRCm38.90). The principal component analysis (PCA) score plot showed a clear separation of High Dose (HD; 7.5 Gy) group from control, however, overlapping of Low dose (LD; 1 Gy) group with control was quite visible (Supplementary Fig. 2). After normalization, a total of 19668 genes were taken for differential gene expression analysis that resulted in the dysregulation of 2837 (1595 upregulated and 1242 downregulated) and 143 (67 upregulated and 76 downregulated) genes in HD and LD irradiated groups, respectively. Supplementary Table 2A and 2B lists the differentially expressed genes with log2 fold change ≥2 and adj. p-value ≤ 0.05 in HD and LD groups post irradiation compared to the control group.
Gene Ontology and KEGG pathway enrichment of transcriptomic data
Gene Ontology (GO)-based enrichment analysis mainly showed perturbation in pathways associated with immune response, cell adhesion, and receptor activity and were also found to be more prominent in HD group (41 in biological process (BP); 11 in cellular component (CC); 10 in molecular function (MF) pathways) than LD group (18 in BP; 4 in CC pathways) when compared to control (Supplementary Table 3A and Supplementary Table 3B).
Among BP in GO, “immunoglobulin production” showed the highest significance in HD group. Among the various differentially expressed genes (DEGs) involved in this pathway, H2-Ab1 (histocompatibility 2, class II antigen A, beta 1), Lax1 (lymphocyte transmembrane adaptor 1), Cd22 (CD22 antigen), Cd40 (CD40 antigen,) Cd28 (CD28 antigen) Vpreb3 (pre-B lymphocyte gene 3), Cd27 (CD27 antigen) Siglecg (sialic acid binding Ig-like lectin G), Clcf1 (cardiotrophin-like cytokine factor 1) and Swap70 (SWA-70 protein) are the top 10 DEGs based on adj. p-value. Detailed categories and sub-categories of dysregulated pathways are presented in Supplementary Table 3 and can be visualized in Supplementary Fig. 3. Visualization of the GO enrichment results were done by GOplot based on z-score (Supplementary Fig. 3). The z-score indicates whether BP, MF or CC are more likely to increase i.e., positive value or decrease i.e., negative value. In other words, it provides expression profile of genes within a GO term.
$$z-{score}=({up}-{down})/\surd {count}$$
Up and down means number of upregulated and downregulated genes in the GO term respectively.
For further analysis, main focus was imparted to BP from the “GO knowledgebase” which includes specifically signaling pathways, biological programs, and cellular functions which could be identified as therapeutic targets. The “Mammalian Metabolic Enzyme Database” (https://esbl.nhlbi.nih.gov/) was utilized to find DEGs (obtained from BP) which encode metabolic enzymes. In HD group, 16 DEGs were found to be encoding those metabolic enzymes that were involved in lipid, nucleotide, amino acid, and carbohydrate metabolism (Table 1). These 16 metabolic enzyme genes are arylacetamide deacetylase (Aadac), 4-aminobutyrate aminotransferase, mitochondrial (Abat), adenylate cyclase type 7 (Adcy7), retinal dehydrogenase 2 (Aldh1a2), plasma membrane calcium-transporting ATPase 3 (Atp2b3), ADP-ribosyl cyclase/cyclic ADP-ribose hydrolase 1 (Cd38), ectonucleotide pyrophosphatase/phosphodiesterase family member 4 isoform X6 (Enpp4), heme oxygenase 1 (Hmox1), inosine-5’-monophosphate dehydrogenase 2 (Impdh2), nitric oxide synthase, brain (Nos1), nitric oxide synthase, endothelial (Nos3), cGMP-specific 3’,5’-cyclic phosphodiesterase (Pde5a), group 10 secretory phospholipase A2 precursor (Pla2g10), calcium-dependent phospholipase A2 precursor (Pla2g5), prostaglandin G/H synthase 1 precursor (Ptgs1), and thymidine phosphorylase (Tymp) (Table 1). However, only one gene encoding to one enzyme cGMP-specific 3’,5’-cyclic phosphodiesterase (Pde5a) could be identified in LD group.
Additionally, KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway database was also utilized to analyze the biological pathways in which DEGs were involved. Based on DEGs data, 18 KEGG pathways in LD group and 12 pathways in HD group were found to be significantly enriched (adj. p-value < 0.05) (Supplementary Table 4). These enriched pathways includes hematopoietic cell lineage, primary immunodeficiency, cell adhesion molecules, and intestinal immune network for IgA production which were found common in both the dose groups. Hematopoietic cell lineage pathway showed highest significance in HD group. DEGs associated with hematopoietic cell lineage pathway can be visualized in Supplementary Fig. 4. Among the 33 DEGs involved, only 9 DEGs (Cd59b, Gp1bb, Itga3, Kit, Itga2b, Il1r2, Gp5, Cd9, Gp1ba) showed upregulation whereas 24 DEGs (Ms4a1, Cd19, H2-Ab1, H2-Eb2, H2-Ob, H2-Aa, H2-Eb1, Fcer2a, H2-DMb2, H2-DMb1, Cd22, H2-Oa, Cd8b1, H2-DMa, Cd38, Cd8a, Cr2, Cd55, Cd5, Cd2, Cd1d1, Cd3d, Cd3g, Cd4) showed downregulation in HD group. Whereas in LD group, only 18 genes were found significant in hematopoietic cell lineage pathway.
Metabolomics/Lipidomics profiling and pathway analysis
UPLC-MS-based untargeted metabolomics was conducted to identify changes in metabolites and lipids in plasma samples post-radiation exposure. During metabolic analysis of plasma, 3426 peaks in positive and 2548 peaks in the negative ionization mode were detected, whereas lipidomics analysis identified 2315 and 1056 peaks in positive and negative ionization mode respectively. Overall, the univariate analysis identified a total of 1491 significant m/z features (lipids and metabolites together; 696 downregulated and 795 upregulated) in LD group and 1221 features (424 features downregulated and 797 upregulated) in HD group (p-value < 0.05). Supplementary Table 5 lists identified metabolites and lipids with mass error in delta (Da). Significant identified metabolites and lipids can be visualized in Fig. 2 and Supplementary Table 6. PCA plot exhibits separation of HD group from controls and some partial overlap of LD with controls and HD can also be seen in both metabolic and lipid profiles (Supplementary Fig. 2). Serine, phenylalanine, histidine, aspargine and proline was upregulated irrespective of exposed dose (Fig. 2). Lipid categories common in both radiation group includes glycerophospholipids [glycerophosphocholines (PC), glycerophosphoethanolamines (PE), LysoPE, LysoPC], fatty acyls (fatty acids/esters (FA), acyl carnitines (CAR)], glycerolipids [tri(acyl|alkyl)glycerols (TG), di (acyl|alkyl)glycerols (DG)], cholesterol esters (CE).
Scatter plot showing significantly altered metabolites (A), lipids (B) and DEGs (C). Data is presented as mean ± standard error of means (number of biological replicates = 4). (For metabolites and lipids: *p-value < 0.05, **p-value < 0.01 and ***p-value < 0.001 using student’s t test (two-tailed); For transcripts * adj. p-value < 0.05, ** adj. p-value < 0.01 and *** adj. p-value < 0.001 using Wald statistics (two-tailed).
Further, pathway analysis demonstrated that a total of 18 and 15 pathways were altered in HD and LD groups respectively. Among them, 4 pathways namely, “arginine and proline metabolism”, “histidine metabolism”, “phenylalanine, tyrosine and tryptophan biosynthesis”, and “glycine, serine and threonine metabolism” showed a significant impact (p-value < 0.05; pathway impact >0.2) in HD group. Whereas, in LD group, only phenylalanine, tyrosine and tryptophan biosynthesis pathway had a significant impact >0.2 (Fig. 3A, Supplementary Table 7).
A KEGG pathway analysis of significantly changed metabolites using MetaboAnalyst. (p values shows pathway enrichment analysis using Global Test and pathway impact values shows the pathway topology analysis). B Joint-Pathway Analysis of significantly altered transcriptomic and metabolomic data after integration. Pathways having impact > 0.2 and –log 10(p) >1.3 (cutoff is highlighted in red dotted lines) were taken into consideration. Each circle represents a single metabolic pathway with area of circle proportional to pathway impact whereas color represents the pathway significance from highest (red) to lowest (yellow). (Enrichment analysis was performed based on hypergeometric test whereas topology measure was performed by degree centrality). C Network visualization of STITCH interactions by Cytoscape in HD group. Networking shows interactions between significant metabolites (p-value < 0.05; square) and DEGs (adj. p-value < 0.05; circle) with the thickness of edge proportional to the interaction score in the STITCH database. Node color shows minimum and maximum fold change. Red, blue and dark blue edge color represents the interaction between metabolite-gene, metabolite-metabolite and gene-gene respectively.
Integration of metabolic/lipidomic and transcriptomic data
The involvement of GO enrichment in metabolic pathways was less evident along with the limited number of identified DEGs that function as metabolic enzymes. Therefore, in order to greatly improve the ability to comprehend transcriptomic data, we employed a multi-omics integrated approach utilizing Joint-Pathway Analysis (MetaboAnalyst), BioPAN software, and STITCH (Search Tool for Interactions of Chemicals) interaction to understand changes in metabolic pathways after IR exposure.
Differential networking by Joint-Pathway Analysis and STITCH interaction analysis
Joint-Pathway Analysis generated 16 and 71 altered pathways in LD and HD groups respectively (Tables 2 and 3). Pathways with high impact >0.2 in HD group included “aminoacyl-tRNA biosynthesis”, “arginine biosynthesis”, “glycerophospholipid metabolism”, “purine metabolism”, “thiamine metabolism”, “linoleic acid metabolism”, “arginine and proline metabolism”, “synthesis and degradation of ketone bodies”, “butanoate metabolism”, and “alpha-linolenic acid metabolism” (Table 3). In LD group, only histidine metabolism demonstrated a significant impact >0.2 (Table 2, Fig. 3B). It is worth noting that more number of metabolic pathways were significantly enriched under Joint-Pathway Analysis using both genes and metabolites as compared to the separate pathway analysis of transcriptomics and metabolomics. Furthermore, the transcript-metabolite interaction network for significant transcripts (DEGs) and the metabolites engaged in Joint-Pathway Analysis was also created using STITCH interaction analysis to attain more comprehensive information. This approach has thus helped in visualizing the relationship between functionally linked metabolites and genes.
The transcript-metabolite interaction in HD group consists of 114 nodes connected by 81 edges, with a clustering coefficient of 0.889. Whereas no significant interaction was observed in the case of LD group. Here, interactions mean associations between gene-gene, metabolite-gene, and metabolite-metabolite and are clearly shown using Cytoscape (Fig. 3C). Associations do not always imply that there is a physical binding present, but they might also contribute to a shared biological function. Network edge represents confidence that indicates the strength of data support. In other words, associations can be based on reactions from pathway databases, literature associations, similar structures or similar activities13. It is apparent that arginine and lysine have shown the maximum number of associations (n = 6) with other metabolites. Among various dysregulated genes present in the Cytoscape network, 3-oxoacid CoA transferase 2 A (Oxct2a) (upregulated), 3-hydroxy-3-methylglutaryl-Coenzyme A synthase 2 (Hmgcs2) (upregulated), and isopentenyl-diphosphate delta isomerase (Idi1) (downregulated) displayed high fold change (>10) and were present in the same cluster. Along with them, nitric oxide synthase 2 (Nos2) was also upregulated 17 times and exhibited an association with arginine and histidine. Altogether, from the above elucidation, it is evident that the integration analysis was successful in identifying differential pathway networks and their linked metabolites and DEGs associated with radiation-induced injury.
Lipid pathway analysis
In the present study, lipid pathways based on lipidomic data were analyzed using BioPAN. The network of lipid subclasses and fatty acid (FA) pathways of LD and HD groups is presented in Fig. 4A, B. Pathway analysis revealed the biosynthesis of PC (glycerolipids and glycerophospholipids) as the most active pathway in both the irradiated groups whereas biosynthesis of PE (glycerolipids and glycerophospholipids) was observed as the most suppressed pathway only in HD group (Fig. 4). These findings clearly indicate common as well as distinct pathways following radiation exposure at different doses. On the other hand, the FA pathway was found to be hampered only in HD group with Elovl5, Elovl6 and Fads2 as the observed predicted genes. Fads2 (fatty acid desaturase 2), Elovl5 (ELOVL family member 5, elongation of long-chain fatty acids), and Elovl6 (ELOVL family member 6, elongation of long-chain fatty acids) were also amongst the DEGs identified in HD group (Fig. 4C). It noticeably displays the complementary information obtained from different omics.
Panels A and B represents active or suppressed lipid pathways in HD and LD groups respectively. Panel C represents the active fatty acid pathway in HD group compared to the control. Green nodes represent active lipids whereas green shaded arrows represent active pathways. (p-value < 0.05 equivalent to Z-score >1.645; represents significantly altered reaction between the control and irradiated group; Green arrow: positive Z-score; purple arrow: negative Z-score).