Introduction

Space travel induces various physiological changes that can affect cardiovascular, metabolic, endocrine, cognitive, and immune functions1. Among the risk factors associated with space travel, space irradiation (IR) from Solar Particle Events (SPE) that contain high energy protons and Galactic Cosmic Rays (GCR) that contain high charge and energy (HZE) particles, is one of the primary health risks2. The effects of irradiation on carcinogenesis and neurodegenerative risks have been relatively well characterized3,4. However, there are limited animal and human studies on the risk of cardiovascular disease (CVD) development due to space IR, especially in the context of future exploration-type space missions.

A systematic review analyzing 35 space-associated human and animal studies investigating the effects of space-IR and microgravity on CVD showed that microgravity causes cardiac atrophy and orthostatic intolerance, while space-IR may potentially accelerate atherosclerosis5. However, the evidence is not entirely conclusive. Another study reported an increased relative risk per Gy for all CVDs and its subtypes per unit of radiation dose (Gy)6. This risk was found to be higher for lower doses and with fractionated exposures4. Furthermore, a recent meta-analysis of 93 studies reported a significant association between high-dose IR exposure and CVD risk. However, the association was less clear at low doses, and only modest differences between acute and chronic exposures were observed6. However, strong inter-study heterogeneity and interpersonal variability complicated causal data interpretation.

Animal studies focusing on the assessment of health risks from exposure to simulated space-type IR are of special importance in estimating CVD risks for future deep-space missions7,8,9,10,11,12. Our previous transcriptional profiling study in female CB6F1/Hsd mice at 16 months after different types of low and very low doses of gamma (137Cs 40–160 cGy) and particle IR (14Si-IR 4–32 cGy or 22Ti-IR 3–26 cGy) revealed 12 common differentially expressed genes (DEGs) associated with CVD, and metabolic diseases as well as cancer and aging across all 3 IR types. This study suggested significant and long-lasting perturbation in gene expression following low and very low doses of IR in the whole heart tissue, in an IR type- and dose-dependent manner.

Our recent study in male C57Bl6/J mice showed acute and chronic changes in LV systolic function in response to γ- and simGCRsim-IR at 14-, 28-, 365-, 440- and 660- days post-IR8. Interestingly, male mice exposed to simGCRsim-IR exhibited a cardiac phenotype resembling heart failure with preserved ejection fraction (HFpEF) at 660 days post-IR, while cardiac function was restored in γ-IR mice. At 660 days post-IR, increased expression of markers of cardiac fibrosis, inflammation, and hypertrophy suggests that simGCRsim-IR might continuously affect cardiac remodeling during the post-IR period8. In contrast, cardiac function in females wasn’t affected during similar post-IR monitoring (unpublished data). These findings indicate sex-depended dichotomy in cardiac function alterations in response to space-type IR.

We hypothesized that biological responses induced by γ- and simGCRsim-IR are chronic, IR type-dependent, and potentially can increase the relative risk for developing CVD during and after long-duration space missions. Further, there may be sex-specific differences in IR-associated alterations in CV function, structure, and gene expression in heart tissue. To the best of our knowledge, there are no published lifetime murine studies comparing gene expression in the hearts of male versus female mice.

In this study, we aimed to evaluate the changes in transcriptome landscape in the LV of male and female C57BL6/J mice at 440/440 and 660/550 days after single-dose full-body exposure with γ- and simGCRsim-IR. Our goal was to determine the potential sex-dependent molecular changes induced by different IR-types and aging, providing a clearer understanding of the underlying molecular mechanisms and the long-lasting effects of these IR in the heart tissue.

Results

Transcriptome landscape

We performed RNA-seq of LV tissue collected at 440/440 and 550/660 days after exposure to a single full-body dose of simGCRsim- and γ-IR; sham/non-IR mice served as controls. Compared to male mice, between 365 and 440 days post-IR, a higher proportion of female mice spontaneously died (i.e., found dead in cages overnight). Autopsies of mice that were found spontaneously dead or euthanized between 365 and 440 days post-IR showed a larger distribution of organ-specific pathologic diagnoses, but none had a high prevalence to provide a unifying etiology for the higher mortality observed within this mouse cohort (Supplementary Tables 1 A, B). Therefore, we limited the collection of female LV samples to 550 days post-IR to ensure an adequate sample size for analysis.

RNA-seq metrics showed an average (mean ± SD) of 47,591,797 ± 5,544,387 100 bp single-end reads per sample. The clean reads/raw reads percentage were (mean ± SD) 98%±0.5; the percentage of Q20 and Q30 bases were 97.88%±0.27 and 93.79%±0.62, respectively. The CG content of the library was on average 47.58%±0.87. Detailed RNA sequencing QC data are presented in Supplementary Table 2.

After raw read preprocessing, we performed a portrayal of transcriptome landscapes using the self-organizing map (SOM) approach13. Evaluation of deregulated spot distribution on the group-centered SOM portraits demonstrated that both types of IR exposure trigger variations in the transcriptome landscape compared to sham mice. Moreover, the highest variability of gene expression changes is associated with simGCRsim-IR at 440 days in males and 440 days in females, suggesting sex-, time point- and IR-type-associated specific SOM spot distributions (Fig. 1).

Fig. 1
figure 1

Visualization of transcriptome landscapes studied animal groups. Inspection of portraits showed considerable differences in gene expression profiles depending on all studied variables (sex: 5 male/ 5 female mice per group; time point: 440 days/550 days/660 days post-IR; irradiation: Sham/simGCRsim/γ), which was especially noticeable between males and females exposed to simGCRsim radiation at 440 days post-IR.

Full size image

To assess the contribution of these factors to the variability of LV gene expression, we performed PCA and UMAP analyses. PCA eigenvectors were used to perform ANOVA for each variable. Analysis indicated that sex and post-IR time were significantly associated with PC4 and PC6 components, respectively, while IR-type was significantly associated with PC24 suggesting a higher contribution of the sex and post-IR time to overall gene expression variability (Fig. 2). This is also confirmed by the visualization of our samples based on the studied variables in the UMAP two-dimensional plot.

Fig. 2
figure 2

The contribution of sex, radiation type, and collection time on gene expression variability. (A) Association of sex, radiation type, and collection time with principal components on the Elbow plot. One-way ANOVA was used to assess the contributing factors and PCs (F-test p < 0.05 considered significant). (B–D) UMAP Clustering based on left ventricular gene expression data. Each dot represents a sample colored according to sex (B), radiation type (C), and collection time (D).

Full size image

Differential gene expression and functional mining

Next, we analyzed sex-stratified DEGs using the DESeq2 R package. We performed 8 comparisons for 12 experimental groups. The complete list of DEGs is presented in Supplementary Data 1. DEG analysis of LV samples from males and females showed significantly different sex-associated responses in DEGs. In total, 54 DEGs (24 upregulated and 30 downregulated) were observed in the simGCRsim-IR males at 440 days (Fig. 3A). In contrast, no differential expression was observed in the γ-IR males at 440 days post-IR (Fig. 3B). At 660 days post-IR, 12 DEGs (6 upregulated, 6 downregulated) were observed in simGCRsim-IR males (Fig. 3C), and only 5 upregulated DEGs in the γ-IR group (Fig. 3D). In females at 440 days post-IR, we detected 65 DEGs (44 upregulated and 21 downregulated) in the simGCRsim-IR group (Fig. 3E) and 59 DEGs (37 upregulated and 22 downregulated) in the γ-IR group (Fig. 3F). At 550 days, 59 DEGs (50 upregulated, 9 downregulated) were detected in simGCRsim-IR females (Fig. 3G), whereas no DEGs were observed in γ-IR females at the same time point (Fig. 3H). Overall, we observed more DEGs in simGCRsim-IR compared to γ-IR mice at any time point examined (Supplementary Data 1).

Fig. 3
figure 3

Analysis of differential expression across the eight comparisons (n = 5 animals per group, 12 study groups) - (A) Sham vs. simGCRsim for males at 440 days, (B) Sham vs. γ-IR for males at 440 days, (C) Sham vs. simGCRsim for males at 660 days, (D) Sham vs. γ-IR in males at 660 days, (E) Sham vs. simGCRsim for females at 440 days, (F) Sham vs. γ-IR for females at 440 days, (G) Sham vs. simGCRsim for females at 550 days, (H) Sham vs. γ-IR for females at 550 days. Each column represents one radiation type and each row one collection time point. DESeq2 package was used to identify significantly differentially expressed genes. The numbers of upregulated and downregulated genes displayed at the top of each plot represent the counts of genes filtered using a Benjamini–Hochberg FDR adjusted p-value (padj < 0.05).

Full size image

Interestingly, overlap of only 6 DEGs was (Tspan4, Or5m3b, Arntl, Spon2, Tef, Mid1ip1) was found between males and females irrespective of the post-IR time point and IR-type (Fig. 4A). However, the direction and magnitude of change was also different across study groups. For example, Arntl was upregulated in females subjected to gamma- and simGCRsim-IR at 440 days (logFC = 1.56, padj = 0.003 and logFC = 1.22, padj = 0.042, respectively), but was downregulated in simGCRsim-IR males at 660 days (logFC=−1.71, padj = 0.003). The same for Tspan2 gene which was upregulated in male at 440 days post-simGCRsim-IR (logFC = 0.59, padj = 0.003) while downregulated in females at the same time point but for γ-IR (logFC=−0.45, padj = 0.04). For Or5m3b, upregulation at 440 days in simGCRsim-IR females (logFC = 0.98, padj = 0.019) was found while the same gene was downregulated in simGCRsim-IR males (logFC = 1.05, padj = 0.019). Further, Spon2 was downregulated in simGCRsim-post-IR at 440 days males (logFC=−1.79, padj = 0.006) and upregulated in females at the same time point (logFC = 1.51, padj = 0.027). Additionally, DEGs varied in both females and males at the same post-IR time point, depending on the IR-type. The upset plot (Fig. 4B) shows additional overlapping DEGs between comparison groups. Bars on the left side of the overlap matrix represent the total number of DEGs for each color-coded comparison group. Color-coded bars at the top of the overlap matrix show the number of unique genes in each comparison group. (Fig. 4B).

Fig. 4
figure 4

Dissimilarity of gene expression in males and females in response to radiation. (A) Venn diagram showing overlap of all differentially expressed genes in males and females (total 12 groups, n = 5 animal per group). The numbers within the diagram represent the unique and shared differentially expressed genes due to radiation exposure from DEG analysis performed separately for male and female data sets. Only six genes between males and females overlapped. (B) Upset plot showing overlaps between all comparison groups. Bars on the left side of the overlap matrix represent the total number of DEGs for each color-coded comparison group. Color-coded bars on the top of the overlap matrix show the number of unique genes for each comparison group (not overlapping with any other group). Gray bars in the top left side of the matrix show the number of overlapping genes between comparison groups. When two or more comparison groups have overlapping genes, the respective green dots are connected to each other in the overlap matrix.

Full size image

Functional enrichment using overrepresentation analysis (ORA) with Enrichr and other tools did not reveal any significant associations with biological processes, signaling pathways, etc., so we performed a search with more mouse-focused tools such as MGI MouseMine and gProfiler. Among upregulated genes in simGCRsim-IR males at 440 days, Tlr4 and Abca1 were associated with phagocytosis (padj = 9.551 × 10 − 3 and padj = 0.00009, respectively). The downregulated DEGs were associated with transition between slow and fast fibers (Gtf2i and Gtf2ird1, padj = 0.032) as well as targets for transcription factors CKROX, MAS and BTEB2 (padj = 0.0017, padj = 0.006, padj = 0.036, respectively). Furthermore, in simGCRsim-IR group at 660 days we observed upregulation of non-tyrosine kinase fibroblast growth factor receptor activity (Fgfrl1,padj = 4.987 × 10 − 2) as well as targets for transcription factor Egr-1 (Cys1, Rhobtb1, Mid11p1, Tef, Fgfrl1, Mknk2, padj = 0.049).

In γ-IR males at 660 days the gene sets associated with immune response (Gimap3, Tcf7, Kbtbd11, Itgb3, Itgb7) were upregulated (padj = 0.00007, padj = 0.00019, padj = 0.0029, padj = 0.03, padj = 0.04, respectively) while no DEGs for the earlier time point (440 days) were detected.

In females post-γ-IR at 440 days upregulated gene sets (mt-Nd4l, mt-Co3,mt-Atp6, mt-Co2) were associated with oxidative phosphorylation (padj = 0.01, padj = 0.02, padj = 0.03, padj = 0.03, respectively). We also detected 4 DEGs associated with circadian rhythm (upregulation of Arntl and downregulation of Dbp, Bhlhe41, and Tef with padj = 0.003, padj = 0.000001, padj = 0.0002, and padj = 0.04, respectively).

In post-simGCRsim-IR females at 440 days, gene sets associated with fatty-acid metabolism (Acot7, Acaa2, Nudt7, padj = 0.038) and circadian rhythm (Tef, Bhlhe41, Hlf, Dbp, Ciart, padj = 0.018) were downregulated whereas the upregulation of genes related to cell junction and assembly (snca, lgi2, fn1, nectin1, sipai1, rapgef1, picalm) was observed. Interestingly, these upregulated genes represent targets for transcription factors Ap2 (padj = 0.03) and Movo-B (Lgi2, Kif1a, Gas7, sh3pxd2b, padj = 0.04) and are also linked to synapse organization and assembly.

Finally, in simGCRsim-IR females at 550 days, the upregulated gene Sirt4 was associated with lipoamidase activity (padj = 0.049) while downregulated genes (Tnfrsf11a, lamp1, pkn1, csf1r, trf, ctsc, ptprj, arf6, etc.) were linked to immune response padj = 0.0000018).

Please note due to word limitation the results of isoform switch analysis are presented in Supplementary Materials, Results section.

Discussion

The results of this study revealed significant shifts in gene expression signatures associated with factors, such as sex and post-IR time. Those showed the highest effect on gene expression variability compared to IR type. The mouse lifetime echocardiography showed sex-specific alterations in functional cardiac changes in the same groups of male8 and female mice (unpublished data). Specifically, in male mice 660 days after 50 cGy simGCRsim-IR, the left ventricular (LV) ejection fraction was preserved, which was associated with a smaller LV size, lower mass, and possible compensation via concentric remodeling (high RWT and low/normal LV mass index) suggesting a development of phenotype resembling heart failure with preserved ejection fraction (HFpEF)8. Interestingly, under the same experimental conditions and the follow-up period, no changes in cardiac physiology were observed in γ-IR male mice8. In contrast, long-term (550 days) cardiac function was not affected in female mice during post-IR (both γ- and simGCRsim-IR) monitoring (unpublished data). These distinct differences in LV physiology in male versus female mice hearts over the murine lifetime have prompted us to examine the life-long sex-specific changes in the transcriptome landscape in LVs of male and female mice.

The RNA-seq analyses showed that simGCRsim-IR induced more pronounced perturbations in the heart transcriptome in both sexes compared to γ-IR. Consistent with our previous published findings in female CB6F1/Hsd mice12, the effects of γ-IR in terms of differentially expressed genes were also evident in female mice at 440 days post-IR in this study, however, at 550 days post-IR no DEGs were detected. Moreover, γ-IR exposure showed only a minor effect on DEGs in male mice with no DEGs found at 440 days and a few detected at 660 days post-IR. In contrast, the number of DEGs in simGCRsim-IR at 440 and 660 days remained considerably high in both sexes at 440/440 and 550/660 days (male/female). These findings suggest that space-type IR exposure may induce long-term genes expression changes and associated (patho)physiological responses in the LV compared to γ-IR.

Interestingly, in males subjected to simGCRsim-IR, we observed differential expression of genes that could be attributed to the previously described physiological alterations. Indeed, the observed association of genes reduced transition between muscle fiber types may contribute to diastolic dysfunction seen in HFpEF with prolonged recovery due to altered muscle fiber properties14. Moreover, in later time points we also observed the increased expression of gene fgfrl1 (fibroblast growth factor receptor 1) associated with fibrogenesis15. Thus, the fgfrl1 gene is known to contribute to fibrotic transformation of tissues.

In simGCRsim-IR females we observed upregulation of target for AP2 and Movo-B transcription factors. This factor is implicated in preimplantation development in mice16. Furthermore, Bcl6 was shown to have a protective function against myocyte cell death17. These changes were paralleled with downregulation of fatty acyl-CoA hydrolase activity that is implicated in the regulation of cardiac heart metabolism18. At the late time point these changes were followed by upregulation of Sirt4 gene which promotes hypertrophic growth, the formation of fibrosis and cardiac dysfunction by increasing ROS levels upon pathological stimulation19.

In female mice examined 440 days after γ-IR exposure, gene sets related to oxidative phosphorylation (mt-Nd4l, mt-Co3, mt-Atp6, mt-Co2) were found to be upregulated, aligning with previous findings that showed an increase in oxidative stress due to radiation treatment20. Additionally, there was an upregulation of mitochondrial dysfunction. We also identified four DEGs related to circadian rhythm, specifically, the upregulation of Arntl and the downregulation of Dbp, Bhlhe41, and Tef, suggesting that the decreased expression of circadian clock genes could be a risk factor for cardiac tissue damage induced by ionizing radiation21. Previous studies have also reported fibrosis resulting from γ-IR exposure22,23. Overall, our findings suggest significant cardiovascular alterations and mitochondrial dysfunction at 440 days post-IR exposure, raising concerns about long-term fibrosis riskt in female mice. At 440 days post-γ-IR, in male mice no DEGs were obtained while at 660 days after γ-IR in LVs of male mice, we detected upregulation of immune response (Gimap3, Tcf7, Kbtbd11, Itgb3, Itgb7). These 5 upregulated DEGs in this animal group can affect cardiac function since integrin (Itgb3) alter the tissue architecture24, while upregulation of the Tcf7 gene is associated with cardiac hypertrophy25.

SimGCRsim is a mixed-field particle radiation protocol that consists of 6 beams of five ions that are designed to mimic, to a certain degree, the space IR environment26 in order to evaluate the excess relative risks (ERR) of a space-type IR on various parameters, including cardiac function and tissue gene expression7,8. While most of the previous studies focused on short-term effects, our results show the persistent footprint of simGCRsim exposure on the LV transcriptome up to 550/660 days (female/male) post-IR.

At 440 days post-simGCRsim-IR, downregulation of fatty-acid metabolism and circadian rhythm-related gene sets was detected in female mice. Interestingly, previously decreased expression of fatty-acid oxidation enzyme in failing heart LV tissue has been detected27. In the same animal group, the upregulation of cell junction and assembly underlying genes were upregulated. At the later time point, immune response was upregulated while lipoamidase activity was downregulated.

At 440 days in simGCRsim-IR males, we observed upregulation of gene sets (Tlr7, Clec4n, Abca1) associated with immune and inflammatory processes including phagocytosis. Among those, Tspan4, a circadian-regulated gene, has a rhythmic expression pattern in multiple organs including the liver, kidney, lung, and heart, among others28. At the same time, we observed downregulation of biological processes associated with RNA processing and calcium transport. Among the downregulated DEGs in this group, two genes, namely, Cacna1h and Cacna2d2 are known as regulators of calcium channels and neurodevelopment29. Interestingly, the Abca1 gene has been reported to inhibit inflammatory signaling pathways and has a key rate-limiting role in reverse cholesterol transport. Furthermore, the upregulation of Dusp19 has been implicated in cardiac hypertrophy and dysfunction30.

The downregulated DEGs were associated with transition between slow and fast fibers (Gtf2i and Gtf2ird1, padj = 0.032) as well as targets for transcription factors CKROX, MAS, and BTEB2 (padj = 0.0017, padj = 0.006, padj = 0.036, respectively). Furthermore, in simGCRsim-IR group at 660 days we observed upregulation of non-tyrosine kinase fibroblast growth factor receptor activity (Fgfrl1,padj = 4.987 × 10 − 2) as well as targets for transcription factor Egr-1 (Cys1, Rhobtb1, Mid11p1, Tef, Fgfrl1, Mknk2, padj = 0.049).

In γ-IR males at 660 days the gene sets associated with immune response (Gimap3, Tcf7, Kbtbd11, Itgb3, Itgb7) were upregulated (padj = 0.00007, padj = 0.00019, padj = 0.0029, padj = 0.03, padj = 0.04, respectively). In females post-γ-IR at 440 days upregulated gene sets (mt-Nd4l, mt-Co3,mt-Atp6, mt-Co2) were associated with oxidative phosphorylation (padj = 0.01, padj = 0.02, padj = 0.03, padj = 0.03, respectively). We also detected 4 DEGs associated with circadian rhythm (upregulation of Arntl and downregulation of Dbp, Bhlhe41, and Tef with padj = 0.003, padj = 0.000001, padj = 0.0002, and padj = 0.04, respectively).

Overall, across all time points combined, there were 174 DEGs in females and 70 in males of which only 6 genes (Arntl, Tef, Spon2, Mid1ip1, Tspan4, and Or5m3b) were overlapping. Almost all overlapping genes were directly or indirectly related to the circadian clock31. Disruptions in the circadian clock have been reported to significantly affect post-IR systolic dysfunction and increased fibrosis in mice21. Among these genes, Arntl (Aryl Hydrocarbon Receptor Nuclear Translocator) and Tef (thyrotroph embryonic factor), are essential components of circadian clock regulation. Artnl is an activator of the circadian clock that regulates physiological functions, including metabolism, sleep, body temperature, and blood pressure32. In turn, Tef, together with other circadian PAR-domain basic leucine zipper transcription factors (comprising of Vbp/Tef, Dbp, and Hlf) is expressed in different peripheral tissues and is known to modulate basal and inducible xenobiotic detoxification33. Another overlapping gene, Spon2, is known for its role in the development of CVD and response to fluctuations in circadian rhythms34. Midline1 interacting protein 1 (MID1IP1), a glucose-responsive gene that colocalizes with c-Myc, is involved in liver cancer growth mediated by ribosomal proteins35. A notable oscillation in both RNA expression and DNA methylation has been reported for Tspan4, a circadian-regulated gene28. Moreover, it has been shown that among other genes, Tspan4 has a rhythmic expression pattern in multiple organs including the liver, kidney, lung, and heart, among others28.

The results of this study point to several observations. On the level of DEGs, we observed sex-specificity in the gene expression response of the LV to terrestrial and space-type IR. This phenomenon has been observed previously, both in human and animal models after terrestrial occupational and accidental exposures36; however, to the best of our knowledge, there are no reports of longitudinal effects for space-type particle IR exposure on heart tissue. Across all time points, only 6 overlapping genes were found. Importantly, in males at 440 days post-simGCRsim we observed downregulation of the slow to fast fibers and at 660 days upregulation of non-tyrosine kinase fibroblast growth factor receptor activity. Moreover, gene sets associated with circadian rhythm were altered at different time points (post-γ-IR males at 660 days, post-simGCRsim-IR males at 440 days, post-simGCRsim-IR females at 440 days).

In summary, our study demonstrates that the effects of a single-dose γ-IR and simGCRsim-IR persist for as long as 550–660 days post-IR in both male and female mice and can be detected on gene expression level in the heart LV. The simGCRsim footprint appears to be more pronounced than γ-IR. Another important observation is that the effects of IR on the expression of individual genes are heavily sex-biased, though, on the systems level, similar processes are affected. Thus, the long-term persistence of IR effects on LV tissues could pose significant health risks, underscoring the need for early development of mitigating strategies. Based on these results, it can be concluded that the effect of IR is long-lasting both in males and females in all age groups, with more prominent differences observed in males at 440 days post-IR.

Methods

Mouse groups and irradiation protocol

All animal procedures were performed in accordance with the guidelines of the Guide for the Care and Use of Laboratory Animals for the National Institutes of Health and approved by the Animal Care and Use Committees at Brookhaven National Laboratory (BNL) (Upton, NY, USA) (BNL IACUC Protocol #502) and the Icahn School of Medicine at Mount Sinai (New York, NY, USA) (ISMMS IACUC Protocol #2019-0017). Sixty C57Bl/6J mice (30 males, 30 females) from Jackson Laboratory (Bar Harbor, ME, USA) were used in this study. We exposed 3-month-old male and female age-matched C57Bl/6J wild-type (WT) mice to 100 cGy, 0.662 MeV of 137Cs-γ-IR ((20 animals, male/female: 10/10), 50 cGy 500 MeV/nucleon simGCRsim-IR (20 animals, male/female: 10/10), whereas control mice (20 animals, male/female: 10/10) were not irradiated. The radiation doses were recommended by the Space Radiation Element (SRE) of NASA’s Human Research Program (HRP) based on the SRE programmatic requirements. Animals were housed in groups of 4–5 per microisolator cage, with water and food ad libitum. Both control non-irradiated (Sham) and irradiated mice were fed commercial mouse chow, referred to as Normal Diet (ND). The animal room was kept at a 12:12 h light-dark cycle at 20–22 °C with 30–70% relative humidity. The detailed experimental design is shown in Fig. 5A and the radiation ions, doses, energy, and animal procedures have been previously described. Animals were monitored at least once a day for any physical or behavioral changes that might indicate distress, discomfort, pain, or injury. Only mice showing general signs of poor health and distress (e.g., rapid heart/respiratory rate, oral/nasal discharge, wound dehiscence, marked swelling, ulcerating neoplasms >2 cm or ulcerating, inability to eat/drink, or weight loss >15%) were immediately euthanized using 100% CO2 at 20% air replacement per minute rate, followed by cervical dislocation. Animals found spontaneously dead in their cages were grossly examined, and their carcasses were frozen. In this study, mice were euthanized via exsanguination at 440/440, and 660/550 (male/female) days after irradiation. Fresh snap-frozen LV tissues from male and female mice were collected at 440 and 660/550 days from control, 100 cGy γ and simGCRsim 50 cGy (5/group, respectively), and stored at -80 °C until further use. All animal procedures were performed in accordance with the ARRIVE guidelines.

RNA extraction and sequencing

Total RNA was extracted from LV tissues using the RNeasy Mini Kit (Qiagen, USA) according to the manufacturer’s protocol and stored at -80°C until further use. A poly-A selected mRNA library was prepared using the NEBNextUltra™II RNA Library Prep Kit. The sequencing was performed on the Illumina NovaSeq-6000 platform to collect an average of 48 million 2 × 150 bp paired-end reads per sample.

Raw read processing, QC, and alignment

FastQC (available online at https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ - version 0.11.9) was used for the quality assessment of raw sequencing reads (Q30% was greater than 94%). Next, we aligned the reads to the mouse reference genome (mm39) with STAR aligner (version 2.7.8a)37 with gene count quantification mode and obtained raw gene counts. On average, 90.48% of the reads were mapped to the reference genome.

Transcriptome landscape portrayal

Transcriptome landscape visualization was performed using the oposSOM R package13. The package self-organizing maps (SOM) machine learning dimension reduction, clustering, and visualization technique to project high-dimensional gene expression data onto a two-dimensional 40 × 40 grid. During SOM training 13,352 gene expression profiles in 60 samples were transferred into 1600 metagenes (gene clusters), each representing a cluster of co-expressed/correlated genes. The metagene expression values of each sample are then visualized (expression portrayal) on the SOM grid using maroon to blue colors corresponding maximum to minimum expression values in each of the portraits12. The sample portraits in each group were further combined to represent a group portrait, highlighting gene clusters deregulated in all samples of a particular group.

Differential gene expression analysis and functional annotation

Raw gene count data was filtered to exclude genes with counts lower than 10 in more than 10 samples. We analyzed differential gene expression from raw gene counts using the DESeq2 R package (version 1.42.1) (https://bioconductor.org/packages/release/bioc/html/DESeq2.html; DOI: https://doi.org/10.18129/B9.bioc.DESeq2). After normalization, scaling, and estimation of global variance, we evaluated differential gene expression in 8 groups as follows − (1) γ vs. sham/control in males at the 440 days post-IR, (2) simGCRsim vs. sham/control in males at 440 days post-IR, (3) γ vs. sham/control in females at 440 days post-IR, (4) simGCRsim vs. sham/control in females at 440 days post-IR, (5) γ vs. sham/control in males 660 days post-IR, (6) simGCRsim vs. sham/control in males at 660 days post-IR, (7) γ vs. sham/control in females at 550 days post-IR, (8) simGCRsim vs. sham/control in females at 550 days post-IR. Genes with pFDR adjusted < 0.05 were considered as differentially expressed. The workflow of this study is shown in Fig. 5B. Functional annotation was performed using Mouse Genome Informatics browser, the international database resource for the laboratory mouse, providing integrated genetic, genomic, and biological data to facilitate the study of human health and disease (https://www.informatics.jax.org/). We also detected isoform switches between groups with a differential isoform usage test and further functional analysis described in more detail in the Supplementary Materials).

Fig. 5
figure 5

Schematic representation of study design. (A) Animal groups and treatment protocol, (B) Experimental procedures and data analysis workflow. Sixty C57Bl/6J mice were used in this study (male/female: 30/30), 5 mice per study group.

Full size image

UMAP clustering and PCA analysis of gene expression data

We visualized the sample clustering patterns based on gene expression data using the Uniform Manifold Approximation and Projection (UMAP) method from the Seurat R package38. To understand the contribution of various parameters to the variability in gene expression in our samples, we performed principal component analysis (PCA). The eigenvectors derived from the PCA analysis were then utilized to conduct ANOVA statistical tests between groups for each clinical variable, thereby elucidating which variables explain the most variability in gene expression. ANOVA was performed for each eigenvector across the groups for every variable, and subsequently, one eigenvector was selected per variable based on the lowest significance.