A multi-omics longitudinal study of the murine retinal response to chronic low-dose irradiation and simulated microgravity

The space environment includes unique hazards like radiation and microgravity which can adversely affect biological systems. We assessed a multi-omics NASA GeneLab dataset where mice were hindlimb unloaded and/or gamma irradiated for 21 days followed by retinal analysis at 7 days, 1 month or 4 months post-exposure. We compared time-matched epigenomic and transcriptomic retinal profiles resulting in a total of 4178 differentially methylated loci or regions, and 457 differentially expressed genes. Highest correlation in methylation difference was seen across different conditions at the same time point. Nucleotide metabolism biological processes were enriched in all groups with activation at 1 month and suppression at 7 days and 4 months. Genes and processes related to Notch and Wnt signaling showed alterations 4 months post-exposure. A total of 23 genes showed significant changes in methylation and expression compared to unexposed controls, including genes involved in retinal function and inflammatory response. This multi-omics analysis interrogates the epigenomic and transcriptomic impacts of radiation and hindlimb unloading on the retina in isolation and in combination and highlights important molecular mechanisms at different post-exposure stages.

Spaceflight hazards such as microgravity and space radiation are anticipated to impact biological systems and could potentially affect crew health and performance both in mission and once astronauts return to Earth. These unique hazards have been thought to induce adverse changes in the central nervous system (CNS), with a frequently reported impact on ocular function [1][2][3][4] . Resources aboard space assets that have a sustained human presence like the International Space Station (ISS) are limited and it is difficult to recreate spaceflight conditions on the ground for controlled experiments. Therefore, characterization of the type and magnitude of these changes is challenging. Despite these challenges, analogs can be used to inform the understanding of the spaceflight environment impacts on biology and support development of mitigation strategies to reduce spaceflight health risks. Ground-based analogs can also support the interrogation of biological effects from disparate stressors in isolation from one another that cannot be untangled in the spaceflight environment. Furthermore, advances in bioinformatics have expanded the potential scientific return of prior experiments provided data is accessible for re-evaluation. Fortunately, NASA's GeneLab is an open-access omics database for space-relevant biological experiments that provides access to uniformly formatted data 5 .
The current study leverages transcriptomic and methylation sequencing data acquired as part of a larger experiment to understand the longitudinal effects and underlying mechanisms of low-dose ionizing radiation (IR) exposure and hindlimb unloading (HLU; an analog for microgravity) on neurovascular oxidative stress and brain function 6 . Global transcriptome profiling with RNA-seq quantifies gene expression at the time of sample acquisition while reduced representation bisulfite sequencing (RRBS) analyzes DNA methylation, which is one of the most studied epigenetic modifications involved in disease development [7][8][9] .
Previous studies analyzed spaceflight-induced transcriptomic and epigenomic response in murine retina 10,11 or evaluated response to different ground-based hazard analogs in murine brain 6 and spleen 12 .
The spaceflight studies reported expression and methylation changes in genes known to be involved in visual perception, phototransduction, and macular degeneration pathways in murine models. However, these studies were performed in space where effects of both space radiation and microgravity on the retina were intertwined. In addition, these studies were performed 35-37 days after a spaceflight and there was no analysis done on longitudinal changes in the retina due to radiation or microgravity.
To our knowledge, the current analysis is the first multi-omics study of longitudinal changes in the murine retina due to microgravity, low-dose radiation, or combined exposure to both. We report differentially expressed genes (DEG), differentially methylated genes (DMG), and biological processes (BP) that are exclusive to, or shared across, multiple exposure conditions and post-exposure time points. We observed similar epigenetic and transcriptomic discordance in genes reported in the NASA Twins Study 13 and in the retina of mice flown to space 11 , with some changes persistent even 4 months after exposure. This study provides an insight into the retinal response to exposures that are frequently used as analogs for individual spaceflight hazard and their interplay at different post-exposure stages, and contributes towards identifying key pathways to interrogate that may provide translational relevance across exposure types.

Results
Differential gene expression. Post-exposure timepoints had different global expression patterns that were visible in principal component analysis where samples clustered more closely based on time point/age of mice instead of the exposure condition ( Supplementary Fig. 1A). For both gene expression and methylation, each experimental group (Supplementary Table S1) was compared to a control group (no radiation or HLU) with matched post-exposure duration, and significant differential expression was identified in eight out of the nine comparisons (Table 1, Fig. 1A,B). While the samples clustered by time point/age based on global gene expression, they could be separated by exposure condition within each time point when genes were selected based on differential expression ( Supplementary Fig. 1B, Fig. 1B).
The number of DEGs peaked at 7 days for groups exposed to only hindlimb unloading (HLU) or irradiation (IR) and showed a subsequent decline. The 275 DEG identified in 7-days IR were enriched in ocular disease terms (Supplementary Table S2) whereas the 71 DEG identified in HLU alone had no significant enrichment of disease pathways. In the 1-month IR group, 7 out of the 31 DEG showed > twofold change (down-regulated: Pon1, Sult1c1, Ermn, Stra6, Tmem72, Gm24514; up-regulated: Bc1). Stra6 and Pon1 are known to be involved in microphthalmia and diabetic retinopathy, respectively and Bc1 affects vascular development 14 and structural plasticity 15 . Prss56, a gene implicated in human and murine refractive development and myopia 16 , was the only DEG at 4 months after exposure to IR alone.
The highest number of DEG within HLU + IR (combination of hindlimb unloading and irradiation) was detected at 1 month post-exposure with no DEG at 7 days. The 64 DEG at 1 month included 20 genes with enrichment (adjusted p-value 8.24e-08) in the somatodendritic compartment of the neuron where synaptic inputs are received and integrated, and 15 genes known to be targets of mmu-miR-466, a miRNA that regulates age-related macular degeneration (AMD) 17 (Supplementary Table S3).
Overlapping differentially expressed genes across exposure groups. Nine DEG were common between the HLUonly and IR-only groups at 7 days (down-regulated: Arl2bp, Chrna6, Egln1, Map2k1, Stmn2; up-regulated: Ccn2, Mbd6, Sox9, St6galnac2, Zfp36l1) (Fig. 1C). Bc1 and Pgk1 were significantly up-regulated in all exposure groups at 1 month. Zinc finger protein Zc3h11a was significantly down-regulated in the HLU only group at 7 days but up-regulated at 1 month. Zc3h11a expression has been shown to be altered in the murine retina and sclera in an experimental myopia mouse model 18 . Only a single common DEG was observed between the HLU + IR and www.nature.com/scientificreports/ HLU groups at 4 months (mt-Co2; up-regulated). mt-Co2 is in the mitochondrial inner membrane and is an essential subunit of cytochrome c oxidase, the terminal component of the mitochondrial respiratory chain 19 . An integrated analysis of mammalian spaceflight data has reported the kidney to show an induction of mt-Co2 20 .
Differential DNA methylation. The 1-month HLU + IR group showed the highest total number (957) of differentially methylated loci (DML) and regions (DMR) in known genes ( Table 1, Fig. 2A). In all exposure groups, hypermethylation was more frequent than hypomethylation. Among the hypermethylated loci, a higher proportion of methylation differences were associated with CpG shores than with islands except in all three exposures at 4 months ( Fig. 2A). Conversely, more hypomethylated loci overlapped with CpG islands than with shores in most of the groups except in all three groups at 7 days. The exposure groups at 4 months contained a higher proportion of total DML within CpG islands than shores compared to the other time points. Top Differentially methylated loci (DML) with an adjusted p-value ≤ 0.05 were selected from each of the nine groups and pairwise correlation for methylation difference was plotted for any given pair. The diagonal shows distribution of methylation difference (%) between exposure and matched control group. Pearson correlation value and significance are displayed on the right of the diagonal, and scatter plots with methylation differences for common DML from any two given groups are displayed on the left. Only significant correlations (> 3 shared DML and p-value ≤ 0.05) are displayed. (C) Upset plot showing differentially methylated genes (DMG), defined as genes with their promoter or gene body containing a DML or DMR (adjusted p-value ≤ 0.05, |methylation difference|≥ 10%), shared between at least two of the nine groups. The Y-axis represents the number of genes in a given intersection set. The X-axis lists different intersection sets and is ordered by the number of intersecting groups in a given set. Dots are colored by exposure condition (magenta for IR, green for HLU, orange for combination). R package MethylKit was used for differential methylation analysis, genomation and GenomicRanges were used for annotation of DML, and ggplot2, ggpairs and ComplexUpset were used for visualization. DML/R/G differentially methylated loci/region/gene, HLU hindlimb unloading, IR irradiation. Correlation in methylation differences across exposure groups. Pairwise analysis of all methylated loci in genes were compared across our 9 exposure groups (Fig. 2B). The highest correlation in methylation difference across different exposures was observed at the same post-exposure time point. Hypermethylation of Bclaf3 promoter was observed in four groups (~ 15% in 1m HLU, 4m IR and 4m HLU + IR, and 26% in 1m HLU + IR), and another CpG island in Tmc8 promoter was common across all three exposure groups at 7 days (10% hypermethylation). The top DMLs shared across all three exposures were in Cav2 at 7 days (30% hypermethylation), Grip2 at 1 month (20% hypermethylation), and Sorbs2 at 4 months (8% hypermethylation).
Overlapping differentially methylated genes across exposure groups. A number of genes containing DML or DMR were observed in multiple exposure groups and post-exposure time points (Fig. 2C). Bcl11b, Bclaf3, Gse1, Necab2, Plec and Tafa5 were each differentially methylated in five out of the nine groups. Six genes (Adamts5, Gse1, Rbm15b, Fscn2, Satb1, Tenm3) were shared across different time points after radiation exposure, with Adamts5 and Fscn2 promoters being hypermethylated even after 4 months of exposure to radiation alone. The promoter region of Pcdh19, a protocadherin, was hypermethylated (> 20% methylation difference) at all three time points for HLU + IR and also in HLU alone at 7 days. Bclaf3 was hypermethylated in five groups (HLU at 7 days and 1 month, HLU + IR at 1 month and 4 months, and IR and 4 months), and CpG islands in Hmgb3 promoter were hypermethylated in HLU and HLU + IR at both, 7 days and 1 month.  Table S4). Processes related to epithelial morphogenesis were activated in the exposure groups, whereas all other processes were suppressed across the three conditions. All shared processes were activated in the 1-month groups, with the exception of axon guidance and cell junction assembly, which were suppressed in the radiation only group. Common processes in the 1-month groups were related to amine transport, regulation of membrane potential, synaptic vesicle exocytosis, and neurotransmitter transport. At 4 months post-exposure, the top processes shared across the three exposure conditions were suppressed. These included lens development in camera-type eye, wound healing, angiogenesis, and sulfur compound metabolic process.

Biological processes affected by low-dose radiation and hindlimb unloading across multiple
Processes that were significantly enriched at all post-exposure time points for a given exposure condition were evaluated (Fig. 3B, Supplementary Table S4). Nucleotide metabolism was suppressed at 7 days and 4 months but activated at 1 month in mice subjected to microgravity alone, whereas epithelial tube morphogenesis and lens development were activated at 7 days but suppressed at all later time points, and detection of light stimulus was suppressed at 7 days and subsequently activated at 1 month and 4 months. Retina homeostasis was suppressed in the radiation only groups at 7 days and 1 month but activated at 4 months, eye morphogenesis was activated at 7 days and suppressed subsequently, and synaptic vesicle endocytosis was suppressed at 7 days but activated at all later time points. In the combination groups, postsynapse assembly was activated at all time points, nucleotide metabolism was activated only at 1 month, and Rho signal transduction was activated at 7 days and 1 month followed by suppression at 4 months. Nucleotide metabolism and related processes were significantly enriched in all nine groups with consistent activation at 1 month and suppression at 7 days and 4 months. Pdk1 was within the gene set for nucleotide metabolism processes across all nine groups, and its expression was decreased at 7 days and 4 months but increased at 1 month in all exposure groups ( Supplementary Fig. 3). Pdk1 is a kinase known to play a key role in regulation of glucose and fatty acid metabolism, and is involved in cellular response to hypoxia 21 .
As a complementary approach to examining shared processes, BP exclusive to each exposure condition (enriched in ≥ 2 post-exposure time points), or exclusive to each post-exposure time point were also examined ( Table 2, Supplementary Table S5). Processes related to microtubule dynamics were enriched exclusively in HLU alone. IR-specific BP included DNA conformation changes, Golgi organization, protein-DNA complex assembly, and type I interferon response. The BP activated exclusively in IR but not enriched in HLU + IR at 7 days included response to wounding, innate immune response, regulation of interleukin-6 production, heart morphogenesis, lipid metabolism, and epithelial cell proliferation. The processes suppressed exclusively in IR were related to telomere organization, DNA repair, ncRNA metabolism, detection of light stimulus, and organelle localization by membrane tethering. Lastly, BP exclusive to HLU + IR were postsynaptic density organization and regulation of AMPA receptors. Each post-exposure time point was associated with distinct BP enriched across multiple conditions only at the given time point. At 7 days, BP related to proteasomal protein catabolism were suppressed in all groups, and mammary gland morphogenesis was activated in IR and HLU + IR. Oncogenesis in the retina and breast have been linked by the regulation of cyclin D1 22 ; the gene encoding cyclin D1, Ccnd1, was significantly up-regulated 7 days after exposure to IR but was not differentially expressed at the later time points. mRNA catabolic process (suppressed) and ion transmembrane transport (activated) were among the BP exclusive to 1 month post-exposure. At 4 months, glycerolipid metabolism, pigmentation, and brown fat cell differentiation were suppressed.
Over-represented biological processes based on differential methylation. Over-representation analysis (ORA) was applied separately to hypo-and hypermethylated genes, and only the hypermethylated genes resulted in significantly over-represented processes (Supplementary Table S6 www.nature.com/scientificreports/ the number of genes shared between a preselected input gene set (e.g., genes differentially methylated based on cutoffs for significance and magnitude of change) and each annotated gene set (e.g., all genes in a GO category), and applying a hypergeometric test to determine the statistical significance of the overlap 23 . Since all the groups showed differential methylation in a sizable set of genes compared to differential expression, ORA was suitable for methylation whereas gene set enrichment analysis (GSEA) was more appropriate for gene expression as it does not require an arbitrary cutoff for selecting significant genes. All genes overlapping with at least one covered CpG site were used as the gene background for ORA. Using a Benjamini-Hochberg adjusted p-value of 0.05 as the cutoff, nine processes were enriched in 4-month HLU + IR and included regulation of Wnt signaling, neurogenesis, nervous system development, cell projection organization, cell surface receptor signaling, pattern specification process and morphogenesis of a polarized epithelium. Genes with hypermethylation and involved in the regulation of these processes included oncogenes and tumor suppressors such as Wnt3, Trp53, Trp73, Sox10, Rap1gap, Paqr3, Kit and Bex1. The promoter region of the brain expressed X-linked gene, Bex1, contained four hypermethylated sites overlapping a CpG island. The gene has been reported to show promoter hypermethylation in malignant glioma cell line specimens 24 . Analysis based on KEGG pathways also showed an enrichment for Wnt signaling (raw and adjusted p-value of 0.005 and 0.3, respectively). Only one other group, 1-month HLU, showed over-representation for a biological process at an adjusted p-value cutoff of 0.05 (actin filament-based process).

Integrated analysis of genes and biological processes impacted at transcriptomic and epigenomic levels.
A total of 23 genes showed changes in expression and methylation at an adjusted p-value cutoff of 0.05 (Table 3). Eef1a1, a protein involved in peptide chain elongation, was down-regulated and hypomethylated in groups exposed to radiation or microgravity 7 days post-exposure. Crystallin Crybb3 www.nature.com/scientificreports/ and signal-induced proliferation-associated Sipa1l3, both implicated in cataract, were up-regulated and hypermethylated 7 days after exposure to IR or HLU, respectively. Fibroblast growth factor 1 (Fgfr1) has been known to mediate photoreceptor rescue effect in response to retinal injury 25 , and the gene was up-regulated and hypomethylated 7 days after radiation. Ncor2, a gene in the Notch signaling pathway, was up-regulated with hypermethylation in the gene body in 7 days IR. Plectin Plec, a cytolinker protein involved in cytoskeletal organization, was up-regulated and hypermethylated. Filamin Flnb was also up-regulated and hypermethylated www.nature.com/scientificreports/ and is known to be involved in microvascular development 26 . No common differentially expressed and methylated genes were found in any of the exposure groups at 4 months. The biological processes and constituent genes shared within a given group with over-representation based on methylation data and enrichment in expression analysis were also obtained ( Table 4, Supplementary Table S7). The highest number of overlapping processes were observed in 4-months HLU + IR, with Wnt signaling pathway being one of the shared BP. Actin filament-based process was enriched in 1-month HLU and 4-months HLU + IR based on both, differential expression and methylation. Processes related to ion transmembrane transport were activated in the 1-month HLU + IR group. Shared genes included ephrin Ephb2, which is also dysregulated in the murine brain following hindlimb unloading 6 . Calcium binding protein Cabp1 and calcium voltage-gated channel subunits Cacna1c and Cacna1h were among the genes involved in ion transport; another calcium voltage-gated channel subunit, Cacna2d4, was also impacted in the murine retina after spaceflight 10 .

Discussion
Spaceflight is an environment to which Earth-based biological systems have not yet adapted. Changes in gravity as well as radiation exposure that differs from terrestrial radiation both in terms of intensity and type, are the two primary naturally occurring hazards the human body contends with during space travel. It is crucial to characterize the impacts of these hazards to ensure the crew and NASA as an agency are appropriately informed. While hazards will be experienced concurrently in spaceflight, important information about the potential interactions between hazards can reveal how impacts can be appropriately mitigated. Different countermeasure strategies will need to be employed if hazards modulate similar or divergent processes.
It is important to note the current study leverages data from a previously performed animal study where mice were exposed to gamma-rays and microgravity as analogs of space radiation and microgravity respectively. However, even with the inherent limitations of using analogs, the observed results still provide insights on processes related to retinal dysfunction due to individual and combined exposures, at different time points, allowing us to gain a preliminary understanding of longitudinal changes. Compared to expression results from other tissue types collected in accompanying studies, spleen data from 7 days post-exposure showed highest number of DEG in HLU alone 12 , brain data from 4 months had most DEGs in HLU + IR 6 , whereas the retina Table 3. Genes with statistically significant expression and methylation changes. Genes that are differentially expressed in a given group (adjusted p-value ≤ 0.05) and also contain at least one differentially methylated locus or region (adjusted p-value ≤ 0.05) are listed. Log 2 (fold-change) for expression difference, percentage methylation difference, and location of DML or DMR within gene context (promoter vs. gene body) are also included. Adjusted p-values for expression and methylation differences are also provided. HLU hindlimb unloading, IR irradiation, DML differentially methylated locus, DMR differentially methylated region. www.nature.com/scientificreports/ results presented here showed DEGs peaked in IR alone at 7 days, and HLU + IR at 1 month. Combined stressors can induce complex responses further regulated by factors such as radiation type, post-exposure duration, and tissue type. Given the conflicting results from various studies addressing synergies between microgravity and radiation and their impact on DNA damage repair and cellular repair 27 , further work is needed to understand their interplay. Biological processes shared across, or exclusive to, exposure conditions and post-exposure stages were identified (Fig. 1, Table 2, Supplementary Table S5), and nucleotide metabolism was enriched in all the groups with Table 4. GO biological processes and related genes common between RNA-Seq and methylation. Top ten BP (or lower if there are < 10 processes) enriched in a given group based on both, differential expression and differential methylation, are listed. The total number of BP for each group, and constituent genes that overlap for a given BP are also included. Full list of BP, constituent genes, and associated methylation difference and fold change values are included in Supplementary Table S7. BP biological process, HLU hindlimb unloading, IR irradiation. www.nature.com/scientificreports/ activation at 1 month, suppression at 7 days and even lower at 4 months for all exposures. The suppression of nucleotide metabolism has been shown to play a role in the establishment and maintenance of the stable growth arrest of oncogene-induced senescence 28 . Further work is required to understand the mechanism behind the activation at 1 month with a return to being suppressed at 4 months, and the specific role of Pdk1 in process modulation. Several genes known to be implicated in ocular diseases were found to be impacted at the transcriptomic and epigenomic levels. A subset of the 7-days radiation only DEG (Cryaa, Cryba1, Cryba2, Cryba4, Crybb1, Crybb2, Crybb3, Crygs, Bsfp1) were enriched in nuclear and congenital cataracts, however, these genes were not dysregulated at 1 month or 4 months. A key component of the MAP kinase signal transduction pathway, Map2k1, was down-regulated in the IR and HLU + IR groups at 7 days but up-regulated at 1 month; MAPK signaling has been shown to be involved in AMD pathogenesis 29 . Bc1, a non-coding RNA known to affect vascular development 14 and structural plasticity 15 was significantly up-regulated in all three exposures at 1 month. Four months after radiation, retention of hypermethylation was observed in genes involved in eye development and pathogenesis including Myc-associated zinc-finger transcription factor Maz 30 , disintegrin Adamts5 involved in AMD, and Fscn2 with a role in photoreceptor disk morphogenesis. Interestingly, Adamts5 was previously observed to be dysregulated in the murine retina after spaceflight 11 . Additionally, early differential methylation (only at 7 days and 1 month) was observed in Hmgb3, a member of the high-mobility group superfamily and a transcription factor known to be impacted by spaceflight 10 .
A total of 23 genes shared differential expression and methylation (Table 3). Of these, 14 had differential methylation in the promoter region and 9 in the gene body with negative correlation between expression and methylation changes in 3 and 6 genes, respectively. The rest showed a positive correlation. Recent studies have expanded the understanding of effects of DNA methylation on gene expression, especially in cancers 31 . Along with the traditional methylation-induced gene silencing, there have been patterns of consistently positive or negative correlations for all CpG sites associated with specific genes. More work is needed to understand the variation in the impact of DNA methylation on expression across different tissue and cell types, genes, and the location of CpG sites. High-mobility group containing Sox9 is known to be involved in retinogenesis 32 and was up-regulated in 7-days HLU and IR groups, and hypomethylated in 7-days HLU. Ppm1a has been shown to play a crucial role in the wound healing-inflammation-angiogenesis axis in mice 33 , and was down-regulated and hypermethylated 7 days after radiation but not later. Sphk1 plays a key role in the regulation of inflammatory responses 34 ; Sphkap, a modulator of Sphk1, was up-regulated and hypermethylated in 1-month HLU + IR. The glycoprotein B4galt6 regulates astrocyte activation during CNS inflammation 35 and was up-regulated and hypomethylated. Interestingly, Sphk1 and glycoproteins B4galt2 and B4galt3 were also up-regulated in 7 days HLU + IR based on spleen data from the same set of mice 12 . Transcriptional co-repressor Tle3, involved in Wnt and Notch signaling pathways, was up-regulated and hypermethylated in 1-month HLU + IR.
Biological processes related to Wnt signaling pathway were significantly enriched in multiple groups based on differential expression, with activation in IR and HLU + IR at 7 days, and suppression in 1-month HLU and all three exposures at 4 months. Wnt signaling pathway is critically involved in cell-cell communication and regulates tissue homeostasis 36,37 . Altered activities may promote tissue degeneration 38 . Wnt signaling plays a role in eye organogenesis 32 and genes in the pathway have been reported to be impacted by spaceflight 11 . The pathway has also been shown to be affected by spaceflight in 3D-cultured neural stem cells 39 and human cardiovascular progenitor cells 40 . Bclaf3 has been proposed to regulate proliferation/apoptosis by suppressing Wnt signaling in the mouse gastric epithelium 41 , and showed hypermethylation in the promoter within four of the groups including IR and HLU + IR at 4 months. Interestingly, genes involved in Wnt signaling also appeared as members in processes enriched in 4-months HLU + IR based on methylation data (Wnt3, Wnt11, Wnt7b, Lrp6, Tnik; Supplementary Table S6).
Among the BP shared between differential methylation and expression (Table 4), pathways related to Wnt signaling were also over-represented in 4-months HLU + IR. Collagen fibril organization was over-represented with three overlapping constituent genes: Col5a1 encodes an alpha chain of type V collagen and the gene family forms a major component of the basement membrane of the corneal endothelium and related genes were altered in the mouse retina after spaceflight 11 , lysyl-oxidase like Loxl1 has been implicated in pseudoexfoliation syndrome which is a major cause of glaucoma and cardiovascular complications 42 , and Comp plays a role in vascular wall remodeling 43 . Actin filament-based process was enriched in 1-month HLU and 4-months HLU + IR. Some of the gene members based on differential methylation in 1-month HLU are also known to be involved in signaling by Rho GTPases (Rhoc, Rhobtb1, Arhgef10l, Cyfip2, Fmnl3, Fgd1, Iqgap3; Supplementary Table S6). Rho GTPases are key players in adaptation to microgravity and actin cytoskeleton dynamics in mammalian cells 44,45 . Temperature homeostasis was enriched in 7-days HLU, and has been shown to be impacted by hindlimb suspension in a previous study where the authors noted that unloading prevented the mice from curling up at night to regulate body temperature 46 . This difference in temperature regulation between the control and experimental groups should be considered while designing unloading experiments and interpreting the results.
Interestingly, the NASA Twins Study reported Notch3 to show epigenetic discordance based on spaceflight blood samples 13 and the gene was significantly up-regulated (twofold change) at 7 days post-irradiation in the mouse retina. Rbx1 and Ncor2, also in the Notch signaling pathway, showed 33% hypermethylation in 4-months HLU + IR, and up-regulation and hypermethylation 7 days post-irradiation, respectively. The NASA Twins Study also reports a decrease during flight in leucine rich alpha-2-glycoprotein 1, a protein involved in retinal vascular pathology. Lrg1 was down-regulated by more than fourfold (unadjusted p-value = 0.006) in 1-month HLU and significantly hypermethylated (16.8% methylation difference) in 4-months HLU + IR.
This multi-omics murine analysis highlights changes in gene expression, DNA methylation, and regulation of biological processes exclusive to, or shared between, exposure groups at different timepoints. Furthermore, there were expressed genes and processes that only presented themselves at the 4-month time point, which may www.nature.com/scientificreports/ be indicative of a longer post-exposure duration to understand the effects of low-dose radiation and microgravity on the retina. Wnt signaling has been shown to be perturbed in the murine retina after spaceflight 11 and our results show discordance in the pathway in the 4-months HLU + IR group based on changes in expression and methylation. Notch3 and Lrg1 showed discordance during flight in The NASA Twins Study blood samples 13 and also exhibit differential patterns in the murine retina in the current study. Hence, despite limitations mentioned below, our study provides confidence that murine models can provide an understanding of molecular mechanisms and pathways contributing to ocular function as a response to analogs for spaceflight stressors. Our study comes with limitations. First, the number of samples was limited with only a total of 59 mice. Therefore, each exposure group at each time point was limited to 3-6 mice, thus limiting power. A follow-up study with more animals should be conducted to confirm these observations. Second, the animals were radiated by low-dose gamma rays which do not represent space radiation. Ground-based low-dose gamma or x-ray radiation studies are often used to study the effect of space-relevant radiation exposures (in terms of total dose) on different tissues/systems 6,47,48 as to inform future space radiation experiments and form an initial understanding of the impacts across radiation doses and types. However, further studies should be implemented with higher dose gamma radiation and GCR simulations to leverage a more realistic space radiation analog. Third, a key limitation of applying bulk RNA-Seq and RRBS to a complex tissue such as the retina with ~ 140 cell types 49 is that the approach averages signals over heterogeneous cell types at different transcriptomic states. As a result, true perturbations in, or correlation between, expression and methylation in a small cell population can be obscured by aggregated profiles. Fourth, caution is needed while comparing results between this study and The NASA Twins Study as biospecimen type and analysis time points differ between studies, and correlating time points across species with such different lifespans is challenging. Therefore, it is important to note that this analysis with its limitations only provides a preliminary signal for further investigation into understanding the longitudinal effect of spacelight hazard analogs on retina.
Overall, the results of this study support the development of a more robust understanding of the longitudinal molecular mechanisms implicated in the murine retina response to low-dose gamma irradiation and simulated microgravity. Further investigation into the interplay between these mechanisms and validation of results in multiple models may provide insights into key pathways with translational potential across biological systems and exposure types as well as targets for mitigation strategies.

Methods
Experimental conditions. The original study was conducted at Loma Linda University following recommendations in the Guide for the Care and Use of Laboratory Animals 50 and was approved on April 30, 2014, by the Institutional Animal Care and Use Committee (IACUC) at Loma Linda University (Protocol number 8130028). A detailed description of the experimental setup has been published previously 51 . Briefly, a total of 59 six-month-old, female C57BL/6 J mice were housed in cages of same dimensions and subjected to one of four experimental conditions: sham control with no radiation or unloading (N = 15), hindlimb unloading (N = 14), low-dose irradiation (N = 15), and a combination of hindlimb unloading and low-dose irradiation (N = 15) for 21 days (Supplementary Table 1). Hindlimb suspension was used to model the unloading, fluid shift, and physiological stress aspects of the microgravity component. Gamma radiation was delivered using 57 Co plates placed underneath the cages and a dose rate of 0.01 cGy/h was used to deliver a total dose of 0.04 Gy over a period of 21 days to simulate a space radiation environment. For the combination group, mice were exposed to unloading and simultaneously irradiated by 57 Co plate placed underneath the cages. Mice from each exposure condition were then allowed a post-exposure period of 7 days, 1 month, or 4 months, following which mice were euthanized and retinas were collected and frozen. RNA and DNA libraries were constructed and used for RNA-Seq and RRBS. The resulting data was deposited to GeneLab where RNA-seq output was processed per NASA GeneLab RNA-seq pipeline 52 , and RRBS results are available in the form of raw sequencing reads. Unnormalized RNA-seq counts and raw RRBS FastQ files were downloaded from GeneLab and used for the current analysis.
Differential gene expression analysis. Unnormalized RNA-seq counts were downloaded from GeneLab. Out of a total of 59 samples (Supplementary Table S1), three were excluded from expression analysis due to low RNA integrity number values or high ribosomal RNA contamination. These counts were generated with GeneLab RNA-seq consensus pipeline, which includes quality control, read trimming, mapping, and gene quantification steps. Full details of the pipeline including the parameters used for each individual tool have been previously published 52 . RNA-seq counts were subsequently analyzed with DESeq2 (v1.30.1) 53 where counts were first normalized using the median-of-ratios method followed by differential gene expression analysis comparing each experimental group to the corresponding post-exposure time point-matched control group. Differentially expressed genes were determined using cutoffs of 0.263 and 0.05 for log 2 (fold-change) and Benjamini-Hochberg adjusted p-value, respectively, to be consistent with analyses of brain and spleen tissues from the same set of mice 6,12 . R packages ggplot2 54 , ComplexHeatmap 55 and pheatmap 56 were used for visualization and generating heatmaps.
Differential methylation analysis. Raw RRBS reads were downloaded from GeneLab and processed with nf-core/methylseq pipeline (https:// github. com/ nf-core/ methy lseq; v1.6.1 accessed on June 23 rd , 2021). In short, FastQC (http:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc) was used for running quality checks on reads, adapter sequences were trimmed using Trim Galore! for Cutadapt 57 without the -rrbs option as recommended for NuGEN RRBS Ovation kit, NuGEN's diversity trimming script was run, and Bismark 58 was used for alignment and methylation call extraction. Deduplication based on overlapping coordinates was skipped as it is not recommended for RRBS 59 . Two samples were excluded from methylation analysis due to www.nature.com/scientificreports/ GC content bias. Mean Phred quality score across each base position in the read was confirmed to be ≥ 20 ( Supplementary Fig. 4), and MethylKit (v1.16.1) 60 was used for detection of differentially methylated loci with a cutoff of 0.05 and 10% for a sliding linear model adjusted p-value 61 and magnitude of methylation difference, respectively 11 . As a complementary approach, methylation differences were also calculated over 100-bp tiling windows to obtain differentially methylated regions 60 . Sites and regions were mapped to known genes (mm10 assembly) and known CpG islands using genomation 62 and GenomicRanges 63 packages. Promoter regions were defined as 2 kbp upstream and downstream of annotated transcription start sites. Differentially methylated genes were defined as genes containing at least one DML or DMR within their promoter region or gene body. For determination of overlap between differential expression and methylation, no thresholds for fold-change or methylation difference were used to allow for genes with modest, but still significant (adjusted p-value ≤ 0.05), changes to be observed.
Gene set enrichment and over-representation analyses. GSEA for gene expression and overrepresentation analysis for methylation were conducted with the ClusterProfiler (v3.18,1) package 64 against the biological process GO category using gseGO and enrichGO functions, and against KEGG pathways using enrichKEGG. Wald test statistics from DESeq2 differential gene expression analysis were used to rank genes for enrichment analysis. Appropriate background gene lists were created by selecting all genes expressed in the retina or covered by at least one CpG site (within ± 2 kb of the gene) for expression (GSEA) and methylation (ORA) data, respectively. A Benjamini-Hochberg adjusted p-value cutoff of 0.05 was used for selecting enriched terms. REVIGO 65 was used for reducing significant biological processes to a representative non-redundant set of terms. ORA was performed separately on hypo and hypermethylated genes to obtain significantly enriched GO biological processes (adjusted p-value ≤ 0.05) based on differential methylation in each group. Following adjustments were made to GSEA and ORA to allow for the discovery of overlapping candidate processes and genes between expression and methylation that can be further validated: (a) differentially methylated genes were defined as genes with their promoter or gene body containing a DML (adjusted p-value ≤ 0.05, |methylation difference|≥ 10%) or DMR (adjusted p-value ≤ 0.05, |methylation difference|≥ 2%); regulated methylation targets are routinely clustered into short genomic regions 66 and a lower threshold for methylation difference in DMR (2%) is recommended for enrichment analysis since perturbed processes can be a cumulative effect of multiple DMR with small effects 67 , and (b) an adjusted p-value cutoff of 0.25 was used for detecting overlapping biological processes based on differential methylation and expression; GSEA guidelines recommend a 25% false discovery rate (FDR) for exploratory analyses as a more stringent threshold can exclude significant results even before they can be validated 68 . Enrichr 69 and ToppGene 70 were used for disease term and functional enrichment analysis.

Data availability
Transcription profiling and methylation data were downloaded from NASA's GeneLab platform (https:// genel ab-data. ndc. nasa. gov/ genel ab/ acces sion/ GLDS-203/; version 6). All software packages are open source and have been cited in the Methods section.