Using prediction models to identify miRNA-based markers of low dose rate chronic stress

• Biopredictors of chronic low grade environmental exposures such as radiation are urgently needed. • Expression analyses of miRNAs (~600) in two mouse strains led to a signature

Robust biomarkers of exposure to chronic low dose stressors such as ionizing radiation, particularly following chronic low doses and dose-rates, are urgently needed. MicroRNAs (miRNA) have emerged as promising markers of exposure to high dose and dose-rate. Here, we evaluated the feasibility of classifying γ-radiation exposure at different dose rates based on miRNA expression levels. Our objective was to identify miRNA-signatures discriminating between exposure to γ-radiation or not, including exposure to chronic low dose rates. We exposed male CBA/CaOlaHsd and C57BL/6NHsd wild-type mice to 0, 2.5, 10 and 100 mGy/h γ-irradiation (3 Gy total-dose). From an initial screening of 576 miRNAs, a set of 21 signature-miRNAs was identified based on differential expression (N± 2-fold or p b 0.05). This 21-signature miRNA panel was investigated in 39 samples from 4/5 livers/group/mouse strain. A set of significantly differentially expressed miRNAs was identified in all γirradiated samples. Most miRNAs were upregulated in all γ-irradiated groups compared to control, and functional analysis of these miRNAs revealed involvement in several cancer-related signaling pathways. To identify miRNAs that distinguished exposed mice from controls, nine prediction methods; i.e., six variants of generalized regression models, random-forest, boosted-tree and nearest-shrunken-centroid (PAM) were used. The generalized regression methods seem to outperform the other prediction methods for classification of irradiated and control samples. Using the 21-miRNA panel in the prediction models, we identified sets of candidate miRNA-markers that predict exposure to γ-radiation. Among the top10 miRNA predictors, contributing most in each of the three γ-irradiated groups, three miRNA predictors (miR-140-3p, miR-133a-5p and miR-145a-5p) were common. Three miRNAs, miR-188-3p/26a-5p/26b-5p, were specific for lower dose-rate γ-radiation. Similarly, exposure to the high

Introduction
There is an increasing demand for the development of novel biomarkers with power to identify low-stress scenarios in relation to environmental, nutritional and life-style-related stressors. Radiation induces cancer as well as a range of non-cancer effects, depending on the dose and dose rate. Nuclear accidents such as Fukushima Daiichi in March 2011, occupational and diagnostic exposure, and the generally increased threat level world-wide for larger-scale human radiation exposure has brought about a renewed demand for rapid, feasible, accurate, specific, and robust biomarkers to identify and quantify exposure of individuals, to facilitate decisions of evacuation or treatment (Borghini et al., 2017;Meineke and Dorr, 2012). Correct triage and dosimetry will prevent unnecessary and potentially life-threatening remediation decisions, such as evacuation-related deaths, as was observed following the Fukushima Daiichi-accident (Nomura et al., 2016).
There is thus a need for specific, feasible, stabile and robust biomarkers of exposure to determine i) whether an individual has been exposed to ionizing radiation including low dose rate radiation, and ii) if exposed, to determine the dose and dose rate the individual have been subjected to. Today there is a lack of suitable biomarkers to determine exposure, especially for lower doses and dose rates. Recently a comprehensive review on the status of biomarker development for ionizing radiation in epidemiological studies was published, emphasizing the shortcomings regarding low dose and low dose rate ionizing radiation (Hall et al., 2017). In this review, an array of potential biomarkers were proposed where non-coding RNAs (particularly, microRNAs) emerged as promising. Others have also pointed to miRNAs as promising biomarker (Malachowska et al., 2020). A powerful strategy to identify novel biomarkers is the combination of measuring an array of stabile biomolecules like miRNAs in conjunction with advanced prediction statistics that exhibit the power to identify predictive biomolecule signatures, which is the objective of this study.
MicroRNAs (miRNAs) are abundant classes of endogenous, small noncoding RNAs (20-25 nucleotides in length), which negatively regulate specific target genes by mRNA degradation or translational repression (Ambros, 2001;Bartel, 2004;Carrington and Ambros, 2003;He and Hannon, 2004). MiRNAs have fundamental roles in multiple cellular processes and are also implicated in the development of multiple diseases (Alvarez-Garcia and Miska, 2005;Cimmino et al., 2005;Meltzer, 2005;Reinhart et al., 2000;Sayed and Abdellatif, 2011). They are involved in the regulation of all cellular functions from differentiation and proliferation to apoptosis, and aberrant miRNA functions can lead to the activation/inhibition of multifactorial physiological processes. The strong stability of miRNAs in almost every biological specimen suggests their promising potential as powerful biomarkers for a wide range of diseases, or exposures. Since miRNAs also represents as promising therapeutic targets and candidate biomarkers in pathophysiology, it is an active area of research.
In recent years, miRNAs have been studied in different scenarios in relation to exposure to ionizing radiation, and they are proposed as markers of exposure to high dose and high dose-rate γ-radiation (Beer et al., 2014;Kraemer et al., 2011), however, the impact of chronic low doses or low dose rates γ-radiation on miRNA expression profiles has scarcely been studied. MiRNA expression profiles have been used to classify different types of cancers, and they may also have the potential to distinguish between low dose-rates γ-radiated and control samples, which is one of the objectives of this study.
We previously reported that the whole blood micronucleus flowcytometry assay is diagnostic of low dose-rate radiation in mice (Graupner et al., 2017;Graupner et al., 2016). In this study our hypothesis is that expression signature sets of miRNAs may serve as rapid, high-throughput predictive biomarkers of exposure, also for chronic low dose rate irradiation. We evaluate the miRNA responses of~600 miRNAs in two strains of mice after continuous exposure to three dose rates of γ-irradiation, including one low dose rate given over a prolonged period. We assayed miRNA expression in liver because this tissue is dominated by one cell type, hepatocytes, that are metabolically active and capable of proliferation. Moreover, radiation induces liver cancer and the liver is considered an important contributor of the miRNAs circulating in the blood stream, which is the preferred sampling medium for human biodosimetry. Thus, liver is suitable for establishing an approach for generating predictor sets of miRNA molecules as well as functional studies, although less relevant for screening purposes.
The main purpose of this study was to identify sets of differentially regulated miRNAs that can predict i) exposure to γ-radiation, ii) exposure to low dose rate γ-radiation, and iii) exposure to high dose rate X-ray exposure. Moreover, we aimed at identifying miRNAs that are differentially regulated due to exposure to ionizing radiation to get insight into their role in the biological response to low dose rate radiation.

Material and methods
2.1. Experimental design: Mice and exposure to ionizing radiation Male (specific pathogen free) mice (8-9 weeks old at start of γirradiation) CBA/CaOlaHsd and C57BL/6NHsd mice (Envigo, The Netherlands) were continuously exposed to γ-irradiation in the Figaro γ-facility (NMBU, Ås, Norway), described in Graupner et al. (2016Graupner et al. ( , 2017 to the following dose rates; 2.5, 10 and 100 mGy/h γirradiation for 1200, 300 and 30 h, respectively, receiving a precalculated total dose of~3 Gy (STable 2). Dosimetry was measured using nanoDots as described earlier (Graupner et al., 2017). The numeric value of the air kerma to whole body absorbed dose conversion coefficient for the chronic exposures was 0.932 ± 0.008 (Graupner et al., 2017), resulting in a whole body absorbed dose of 2.60 ± 0.19 Gy for the 2.5 mGy/h-group, 2.67 ± 0.16 Gy for the 10 mGy/hgroup and 2.65 ± 0.13 for the 100 mGy/h-group, all denoted as~3Gy throughout the article. The animals were kept at 21 ± 2°C at 45 ± 15% relative humidity in individually ventilated cages (Innovive, San Diego, USA) at 50 air changes/hour. Five mice were housed per disposable PET plastic cage with aspen bedding (Nestpack, Datesand Ltd., Manchester, UK). Water (tap water in PET water bottles) and feed (SDS RM1, Special Diet Services, Essex, UK), were given ad libitum to all mice. The continuous irradiation was interrupted daily for 30 min -2 h for animal care purposes, and the beam-on time was correspondingly adjusted to achieve the pre-calculated total dose of 3 Gy. All cages were moved one position to the right in the racks daily to assure uniform exposure of cages throughout the whole irradiation period. Unexposed control mice were outside of the irradiation field behind lead shielding in separate cage racks. Separate groups of mice (10 or 14 weeks of age) were acutely exposed to X-rays (whole body absorbed dose of 2.6 ± 0.1 Gy (~3 Gy) total dose; 117.6 s at 1.51 Gy/min absorbed dose to water, 225 kV, 13 mA, 0.5 mm Cu-filter, with the numeric value dose conversion coefficient of 0.87 ± 0.04 for the absorbed dose to water to whole body absorbed dose, in an X-RAD 225 irradiator from Precision X-ray, North Branford (Graupner et al., 2017)).

miRNA expression analysis
Quantitative real-time PCR (qPCR) was used to analyze the miRNA expression response of mice liver tissue samples. Total RNA was isolated from fresh frozen liver tissues using Quick-RNA™ miniPrep kit, cat. # R1055 (Zymo Research, Nordic BioSite, Norway) according to the manufacturer's instructions. In brief, liver tissue samples (~46 mg) were homogenized in 600 μl lysis buffer plus 5 mm stainless steel-bead by TissueLyser II using following program 3 × 2 min at 20/s (Qiagen, Hilden, Germany). Homogenized samples were centrifuged for 15 s at 10,000 ×g and supernatants were transferred into new 1.5 ml Eppendorf tubes. The total RNA was isolated from the supernatants. Liver RNAs were isolated from five animals from each treatment group (n = 5; i.e., 5 animals/treatment group/mouse strain; 25 liver samples per mouse strain; i.e., in total 50 liver samples). The quantity and quality of isolated RNA was determined as previously described (Aarem et al., 2016;Duale et al., 2014) using a NanoDrop Spectrophotometer (Thermo Fisher Scientific, Massachusetts, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, California, USA). RNA purity was estimated by examining the OD 260/280 and the OD 260/230 ratios. RNA integrity numbers (RIN) from 1 to 10 (low to high RNA quality) were calculated using the 2100 Expert software (Agilent Technologies, Santa Clara, California, USA). RNA samples were stored at −80°C until further use.
The cDNA synthesis was performed with 1000 ng total RNA from samples as template as previously described (Aarem et al., 2016), using the miScript II RT kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol. The cDNA concentration and quality were assessed by a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Massachusetts, USA). All cDNA samples were stored at −20°C prior to miRNA expression analysis.
MicroRNA specific qPCR analysis were carried out as previously described (Aarem et al., 2016;Duale et al., 2014) in 384-well plates using miScript SYBR Green PCR Kit according to the manufacturer's protocol (Qiagen, Hilden, Germany) on a CFX384 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, California, USA). In brief, for an initial miRNA expression profile screening, cDNA (1:80 dilution) from each treatment group, i.e., 3 mice/treatment group; in total liver samples from 12 mice, were analyzed. All liver samples and several different miRNA assays were simultaneously measured in one 384-well plate. This qPCR layout allowed simultaneous measurement of all samples in one run, reducing errors due to run-to-run variations. The expression profile of 576 miRNA assays was performed using this qPCR layout. The cycling program included an initial enzyme activation step at 95°C for 15 min, and then 40 cycles of denaturation, annealing and extension steps at 94°C for 15 s, 55°C for 30 s and 70°C for 30 s, respectively. The melting curve (Tm) analysis was included in each run. Nontemplate controls (NTC) were included in each run. Subsequent to the initial miRNA expression profile screening, we measured the expression levels of 21 selected target miRNAs and four stably expressed reference miRNAsbased on the initial screening resultsin all 49 liver samples by qPCR; i.e. 4-5 mice/treatment group/mouse strain.
For the initial miRNA screening data, the raw data Cq-values from 576 miRNAs were pre-processed and miRNAs with inadequate measurements, with multiple Tm peaks and Cq-values above 30 and Cqvalues below 15 were removed from downstream analyses. In addition, filtering criteria for missing values was set to 70%, and all patterns with b70% existing values were removed. The outcome of the quality assurance filtering criteria was an expression matrix consisting of 197 miRNAs × 12 samples, and these miRNAs were used in the downstream analysis. The Cq-values were normalized with mean expression value for individual samples; i.e., ΔCq (sample) = Cq (target miRNA) -Cq (geometric mean Cq-values of all expressed miRNAs). The normalized relative expression between exposed and control samples were then calculated as following: ΔΔCq (sample) = ΔCq (sample) -ΔCq (mean control) and then the ΔΔCq-value was transformed to linear scale; i.e., fold change (FC) = 2 −ΔΔCq(sample) . All miRNA expression data are reported either as FC or as log2-FC between the control and the exposed samples.
From the initial screening 21 target miRNAs were selected based on statistically differentially significance, fold change N ± 2.0 in one of the treatment groups and/or treatment-related trend, and four reference miRNAs based on NormFinder stability test (Andersen et al., 2004). Of the four reference miRNAs, the two most stable miRNAs (miR-17-5p & miR-431-5p) according to NormFinder, were used to normalize the selected 21 miRNAs. The Cq-values of the 21 miRNAs were normalized as mentioned earlier using the geometric mean of miR-17-5p & miR-431-5p.

Statistical methods
To identify miRNAs with a statistically significant difference in expression compared to the control group, the initial miRNA screening data (576 miRNAs and 12 samples) was analyzed by one-way ANOVA followed by post hoc Dunnett's test. ΔCq values were log transformed. Ordinary linear regression was used to identify any log-linear effect of miRNA levels (Log(Fold-Change) on log(dose-rate+1)). These tests were used to select the final set of miRNAs for qPCR analysis of the full set of samples.
Expression of the final set of 21 miRNAs was analyzed for significance using the limma package in R (Ritchie et al., 2015). MiRNAs were regarded as differentially expressed when the fold change N ± 1.5 and the adjusted p b 0.05 (Benjamini-Hochberg (BH) correction). Pearson's correlations and principal components were inspected for clustering of miRNAs.
Logistic regression was used to identify miRNAs that could classify mouse samples as γ exposed or not for each dose-rate level. In a multivariate approach to classify samples into the correct dose rate group, we applied a set of nine methods: The nearest shrunken centroids (PAM (Tibshirani et al., 2002)), random forest, bootstrap tree and generalized regression (elastic net, lasso, double lasso and their adaptive counterparts). K-fold cross validation with 5 folds was used for all methods except random forest and boosted tree. (Details on the prediction run in Supplementary material.) Random forest and boosted tree were run with different validation set ratios, and with a tuning design with 120 combinations of tuning parameters (including bootstrap options for random forest). No tuning design calibration could give a stable set of predictors with any validation sample size. Thus, only the results from the training sets are reported. Splitting the sample in training set and validation set according to mouse strain was attempted, but the sample size was too small to give any meaningful models. Partial contribution of each miRNA and the prediction accuracy with the different methods are presented. PAM, generalized regression and boosted tree also select the minimum number of miRNA predictors (variable reduction) by adding or removing predictor miRNAs from the model until misclassifications are at a minimum. The various prediction sets from the methods were compared and evaluated.
This set of prediction methods were selected because they could give reasonably consistent results with this small data set. Initial attempts to analyze the data set with PLS, K-means clustering, naïve Bayes were not successful. Neural networks would always classify the samples without misclassifications, but the solutions differed widely. The sample size is too small for the caret or tidymodels packages in R. Besides, JMP has all models easily available, while using the R packages requires a lot more effort.

miRNA profiling
To identify potentially modified miRNA expression levels following low dose-rate γ-irradiation, we investigated the expression profiles of 576 miRNAs in twelve randomly selected liver samples (i.e., 3 mice/ treatment group) from mice of both strains exposed to 0, 2.5, 10 and 100 mGy/h γ-irradiation, respectively. Using stringent selection criteria, i.e., removing miRNAs with inadequate measurements, miRNAs with multiple Tm-peaks, Cq-values N35 and miRNAs with expression pattern with b70% existing values, we obtained a miRNA informative matrix of 197 miRNAs × 12 samples (data not shown) that was used in the downstream analysis. The normalized data were used to calculate FC between γ-irradiated samples and control samples (SFigure 1).
From the results of the initial miRNA screening, miRNAs with expression level difference over ±2-fold or p b 0.05 between exposed and control samples were selected for single sample analysis in all mice: Based on these selection criteria, we identified a panel of 21 signature miRNAs (SFigure 1). The miRNA expression level of the identified 21 signature miRNA panel were subsequently analyzed in all samples (N = 39 samples; i.e., 4-5 mice/treatment group/mouse strain in liver samples from 39 mice) of the experiment. MiRNAs were called as statistically significantly differentially expressed when the fold change N ± 1.5 and the adjusted p b 0.05 (Benjamini-Hochberg (BH) correction for multiple testing). Seven miRNAs (miR-125b-5p, miR-130a-3p, miR-140-3p, miR-145a-5p, miR-181a-5p, miR-455-3p and miR-499-5p) were identified as statistically significantly differentially expressed (FC N ±1.5; adjusted-p b 0.05, BH-method) in all γ-irradiated groups (i.e., both mouse strains) compared to the control group (Table 1). Two other miRNAs miR-126a-3p and miR-133a-5p were also significantly differentially expressed (FC N ±1.5; p b 0.05) upon γirradiation, but they were not significant in all γ-irradiated groups after BH correction (Table 1). These seven miRNAs including miR-126a-3p and miR-133a-5p were all upregulated following γirradiation. Furthermore, when we searched for significantly differentially expressed miRNAs for each treatment group separately, we identified eleven miRNAs in the 2.5 mGy/h exposed group, seventeen miRNAs in the 10 mGy/h exposed group and nine miRNAs in the 100 mGy/h exposed group, following γ-irradiation (Table 1). Upon γirradiation, most of significantly identified miRNAs were upregulated compared to the control group, except one miRNA (miR-23a-5p) which was downregulated in the 10 mGy/h group (Table 1).

Functional enrichment analysis
To identify biochemical signaling pathways affected by γ-irradiation induced miRNA expression changes, we analyzed the predicted target genes of the seven significantly affected miRNAs in both mouse strains ( Table 1). The target genes of the seven identified miRNAs were extracted from miRWalk database (Dweep and Gretz, 2015;Dweep et al., 2011), and the intersection of identified target genes from at least seven prediction programs was chosen. To investigate which pathways may be affected by the seven dysregulated miRNAs, we evaluated the biological functions of their predicted target genes using WebGestalt tool (Wang et al., 2013;Wang et al., 2017), and correlated their predicted target genes with the KEGG (Kyoto Encyclopedia of Genes and Genomes) biochemical pathways (Kanehisa, 2009) in order to identify enriched pathways. The results from enrichment analysis represent a global picture of pathways that are significantly enriched with target genes for dysregulated miRNAs following γ-irradiation exposure. Significantly (FDR b 0.05) enriched top 20 KEGG pathways annotating the predicted target genes for the seven dysregulated miRNAs following γ-irradiation are presented in Fig. 2. The KEGG enrichment analysis indicated that the seven miRNAs targeted genes were involved in a range of signaling pathways. Of the top 20 enriched KEGG pathways, miR-125b-5p was the only miRNA whose target genes were enriched in all 20 KEGG categories. KEGG enrichment analysis showed a link between the γ-irradiation induced miRNA target genes and several cancer-related pathways (Fig. 2) such as pathway of cancer, Wnt signaling pathway, MAPK signaling pathway, ErbB signaling pathway, apoptosis, PI3K-Akt signaling pathway, and p53 signaling pathway.

Correlation analysis
Each miRNA can regulate numerous target genes and therefore has the potential to alter multiple biochemical pathways. MiRNA target gene sharing principle is based on that if two miRNAs share a common set of target genes, they may probably influence or co-regulate similar biological pathway(s). To investigate how miRNAs co-regulate biological processes together, we first conducted correlation analysis of the 21 miRNA panel and the result is presented in Fig. 3. Visual inspection of the correlation heat map (Fig. 3) reveals that the 21 miRNA panel clustered together into three main clusters: cluster 1 (seven miRNAs) and cluster 3 (two miRNAs) mainly consisting of miRNAs significantly expressed in only one or none of the γ-irradiated groups; while twelve miRNAs cluster together in cluster 2. The seven statistically differentially expressed miRNAs in all γ-irradiated groups (Table 1) were among the twelve miRNAs in cluster 2 (Fig. 3).
Correlated miRNAs may have similar target genes and miRNAs targeting the same genes may infer a broader range of target level alteration. To investigate shared target genes among the miRNAs in cluster 2, we focused on the identified seven significant miRNAs in all γirradiated groups in cluster 2 (miR-125b-5p, miR-130a-3p, miR-140-3p, miR-145a-5p, miR-181a-5p, miR-455-3p and miR-499-5p). The predicted target genes of the seven miRNAs were compared and the number of shared target genes for the seven miRNAs are presented in STable 1. Three miRNAs (miR-181a-5p, miR-125b-5p and miR-130a-3p) have N600 predicted target genes. The miRNA pairs sharing the highest number of target genes was miR-181a-5p and miR-499-5p, and these two miRNAs co-target N60 target genes (STable 1). This analysis indicated that some genes may be regulated by several miRNA following γ-irradiation and may explain the regulatory role of dysregulated miRNAs in biological process.

Prediction of dose rate
To identify potential miRNA predictive markers of γ-irradiation, we first conducted linear regression analysis, and the linear regression of dose rate on miRNAs revealed several significant linear effects (SFigure 2). However, the R 2 and beta coefficients were small which indicates that there is little general linear effect of dose rate on miRNAs. Because many of the miRNA transcript levels are correlated (Fig. 3), a principal component analysis of the relative miRNA transcript values was conducted to reveal similarities and contrasts between the miRNA predictors. The partial contribution plot (Fig. 4A) shows that the principal component 1 and 2 separates groups of miRNAs quite well into two clusters. Principal component 3 (blue bars) consists mainly of two miRNAs (23a-5p and 758-3p) that can also be seen close to the center of the biplot (Fig. 4B). The biplot also shows that the control group is separate from the three exposed groups indicating that using the full set of the 21 miRNA panel in prediction analysis may improve the prediction of the exposed group.
The next approach was to investigate whether the expression values of the 21 miRNA panel could separate the three dose groups and controls when considered as discrete categories by clustering techniques. PAM, K-means clustering, or naïve Bayes could not identify sets of miRNAs from the panel that could separate the dose rates without misclassifications (data not shown). Further, a linear prediction model was applied, considering the dose rates as a linear dose-response to the set of miRNAs with lasso and elastic net. None of these methods could come up with significant predictors. Finally, three contrasting categories were fitted separately: controls versus 2.5 mGy/h, controls versus 10 mGy/h, and controls versus 100 mGy/h. Initial logistic regression revealed up to 15 significant predictors (p b 0.05) for exposure at the three dose rate levels (SFigures 3, 4 and 5). With a few exceptions, the probability of samples being classified as irradiated increased with increased miRNA transcript level. However, none of the miRNAs could alone predict exposure with high sensitivity and specificity (i.e., without misclassifications).
We combined the whole set of 21 miRNAs as predictors of each contrast and analyzed them using several prediction methods (generalized logistic regression, random forest, boosted tree and PAM). The rsquares, p values and misclassification numbers for the prediction methods are presented in Table 2. The results from neural networks and partial least squares (PLS) are not shown because of the results changed for each repeat of the analyses, while K-means clustering and naïve Bayes could not present any prediction results. Several models revealed significant effects (Mean -log p b 0.05), and the predictions for the lowest dose rate are generally better than the highest dose rate. Note that the high significant boosted tree values are for the training sets only, cross validation weakens the predictions considerably. The contribution from each of the miRNAs in the prediction analysis is shown in Fig. 5, and fifteen miRNAs could predict one, two or all three γ-irradiated groups. Among the top 10 miRNA predictors which contribute highest for each three contrasting categories (unexposed controls vs 2.5 mGy/h; unexposed controls vs 10 mGy/h; and unexposed controls vs 100 mGy/h) included three miRNA predictors (miR-140-3p, miR-133a-5p and miR-145a-5p) common for all γ-radiated groups (Fig. 5). The average contribution across all models (Fig. 5B) illustrates the ranking of the miRNAs according to importance. The lower right panel sums up the average contribution of each miRNA across all contrasts.

miRNAs affected by X-ray exposure
Separate groups of mice (n = 10 mice) were exposed acutely to high dose rate X-rays (91,200 mGy/h, 3 Gy total dose) in order to compare and verify predictor miRNA signature panel. The miRNA expression levels of the identified 21 miRNA panel were analyzed in all samples (N = 10 liver samples; i.e., combined both mouse strains; 5 CBA/ CaOlaHsd and 5 C57BL/6NHsd mice). In the acutely exposed X-ray group, fifteen miRNAs were significantly differentially expressed: eleven upregulated and four downregulated miRNAs (Table 3 and SFigure 6). Of the identified significantly expressed miRNAs in the Xray group, six of them are among the seven common miRNAs identified in the γ-irradiated groups (Table 3 and SFigure 6). Then, the same prediction methods as those exposed to γ-radiation were used to identify candidate miRNAs that predict exposure to X-rays. Mice exposed to X- .5 are marked with green arrow or red arrow, bold and italic marked miRNAs (N = 6 miRNAs) in column 1 are statistically (p-value and adjusted p-value b0.05; BH-method, limma-method) significantly differentially expressed in all γ-irradiated groups. p-Values (pval) and adjusted p-values (adj.pval) b 0.05 are indicated with pink, green or yellow depending on dose rate group. *mmu-miR-126a-3p and mmu-miR-133a-5p were significantly differentially expressed (FC N ±1.5; p b 0.05) limma-method, but adjusted pvalue N0.05. MiRNA was regarded as statistically significantly differentially expressed when the fold change N ±1.5 and the adjusted p b 0.05 (BH correction for multiple testing). rays were correctly classified using the full set of the 21 miRNA panel (Fig. 6A). The average contribution across all models (Fig. 6B) illustrates the ranking of the miRNAs according to importance, and miRNA predictors which contributed highest included miR-122-3p and miR-125b-5p.

Discussion
Our study has identified a 21 miRNA signature panel of a total of 576 miRNAs in mice liver tissue samples, and the expression level of these miRNAs were analyzed by nine prediction methods to identify whether a mouse have been subjected to oxidative stress such as ionizing radiation, including low grade chronic exposure or not.
In this study, two different mouse strains were used to evade strainspecific results and generate miRNA predictor sets that are general for more than one specific mouse strain. This design was chosen to account for potential differences in strain-associated radiation susceptibility. To mimic a scenario where individuals are subjected to chronic low grade oxidative stress or low dose rates ionizing radiation for longer periods, mice were subjected to low dose rate γ-irradiation (2.5 mGy/h, which is below the threshold (b 6 mGy/h) defined by The United Nations Scientific Committee on the Effects of Atomic Radiation (UNSCEAR) (UNSCEAR, 2010) and others (Brooks et al., 2016), or to higher dose rates (10 and 100 mGy/h) in a continuous manner in the Figaro γfacility (Graupner et al., 2016). The mice obtained a total dose of 3 Gy (Fig. 1), demonstrated in previous experiments to give rise to genotoxic effects (Graupner et al., 2017). Dose-response-related differentially expressed miRNAs were analyzed in liver, since this tissue is metabolically active, mainly consist of one cell type (hepatocytes), it is capable of proliferation and probably contribute with a significant portion of the circulating miRNAs present in blood plasma. Moreover, most solid cancers including liver carcinomas are significantly associated with radiation (Grant et al., 2017;Preston et al., 2007;UNSCEAR, 2010).
miRNAs can be measured in many biological materials such as exosomes (Szatmari et al., 2019), plasma, saliva (Gai et al., 2018), tears and can also be extracted from paraffin-embedded materials (Anastasov et al., 2012), which opens up possibilities for assessing miRNA expression in historical materials. MiRNAs are quite stable in biological specimens after storage (Mraz et al., 2009). MiRNAs are also quite stable in the host organism after ionizing radiation compared to other biomarkers of radiation exposure (Tomasik et al., 2018). Note. Generalized regression methods with K-fold cross validation provide the best models. Random forest and Boosted tree were run without cross validation. High R-Square is better prediction, Mean -Log p indicates significance. Green bars indicate the level of R-square. Note. MiRNAs with FC above ±1.5 are marked with green arrow or red arrow, p-value and adjusted p-value b0.05; BH-method significantly differentially expressed in X-ray groups (limma-method). P-values (pval) and adjusted p-values (adj.pval) b 0.5 are indicated with pink, green or yellow depending.
We identified seven statistically significantly differentially expressed miRNAs present in all γ-irradiated groups (i.e., both mouse strains) compared to the control group (Table 1). The target genes of these seven miRNAs were involved in several cancer-related pathways. The seven miRNA target prediction and pathway enrichment analysis demonstrated that a single miRNA can affect a large number of genes and its targeted genes are involved in several signaling pathways. Furthermore, several miRNAs collectively work together to regulated set of genes, and these genes likely participate in common signaling pathways. Comparison of shared target genes for the seven miRNAs is presented in STable 1. MiRNAs targeting the same genes may infer a broader range of target level alteration. Regulation of a target gene may depend on the combined effect of multiple miRNAs, and therefore, some target genes may be regulated by several miRNA following γ-irradiation. The outcome of miRNA dysregulation upon γ-irradiation may have marked effects in many biological functions.
The miRNA transcriptome profile in irradiated TK6 cells given 2 Gy X-rays at six different times after exposure were assessed (Chaudhry et al., 2013), and our predictive miRNA set were among the miRNAs identified as differentially expressed at the latest time point in this study. Others have used gene expression as approach to identify signatures used as radiation biomarkers, following acute single exposures of cells (Macaeva et al., 2016).
We then evaluated the potential biological relevance and importance of the predicted target genes of the 21 miRNA panel. The KEGG enrichment analysis (presented in SFigure 7) indicated that the 21 miRNAs targeted genes were involved in a range of signaling pathways. The predicted target genes have been implicated having important roles in several cancer-related and radiation-induced cellular pathways (SFigure 7). Overall, there were several overlapping pathways enriched in both the seven significantly identified miRNAs and the 21 miRNA panel target genes ( Fig. 2 and SFigure 7).
Using nine prediction methods, we identified miRNA predictors that distinguished samples exposed to γ-radiation from controls. In contrast to differential expression analysis where each contrast is evaluated separately, prediction models will choose the most important miRNA from possible clusters of correlated miRNAs. If there are many highly correlated variables, the prediction models will select the one that explains most of the dependent variable (a probability in logistic models). The second variable to be included will be the one that is second best, given the previous choice. Although the algorithms for the methods used here are highly different, the outcome of the analysis is always that each method comes up with a combined set of variables that can separate controls from exposed mice with maximum confidence and preferably without misclassifications. Generalized regression methods frequently selected subsets of miRNAs in the model that correctly classified the mouse samples with no misclassifications (100% sensitivity and specificity). Neither boosted tree nor random forest presented stable predictions when cross validation and parameter tuning was applied. Hence, we present these predictions without cross validation. The nearest shrunken centroids method could not come up with a set of miRNAs that could predict without misclassification. The relative contribution of the predictor miRNAs differs between the estimation methods, but a few miRNAs appear consistently in all methods with similar contribution. Generalized regression seems to perform well when classifying radiated and non-irradiated mouse liver samples. The adaptive lasso and adaptive elastic net methods perform better than elastic net (and lasso, data not shown). Adaptive elastic net tends to select a higher number of predictors than adaptive lasso, but fewer than elastic net. It can be seen that most of the extra predictors selected by elastic net are also selected by random forest, boosted tree and PAM. The predictions are most consistent at the low dose rates, while the prediction at 100 mGy/h is weak (Fig. 5), including the generalized regression methods. Considering that this is a limited dataset with few and correlated predictors, it seems appropriate also to explore the potential role of those miRNAs called statistically non-significant (limma analysis) as biomarkers of radiation.
From the prediction analysis, miR-26a-5p appears as the most important predictor after exposure γ-irradiation (Fig. 5B) and it is the highest predictor for 2.5 mGy/h and 10 mGy/h. Interestingly, miR-26a-5p is just slightly significant in simple logistic regression at 2.5 mGy/h (not significant when using limma analysis) and much more significant at 10 mGy/h in both limma and logistic regression analysis. miR-26a-5p also appeared as a significant predictor at 100 mGy/h, but this was inconsistent (data not shown). This illustrates that even slightly dysregulated miRNA expression may have significant impact after environmental stress. MiR-26a-5p belongs to the mir-26 family together with miR-26b-5p, and the mir-26 family is implicated as potential tumor suppressor roles. These two miRNAs, miR-26a-5p and miR-26b-5p are located at different chromosomes; 9 and 1, respectively. MiR-26a has been associated with mouse liver function where upregulation of miR-26a reduce hepatocyte proliferation and induce Fig. 1. Experimental exposure design: Male mice of two different strains were continuously exposed to three dose rates of γ-radiation, or to acute high dose rate X-rays, to an accumulated total dose of 3 Gy.
signs of liver damage such as increased aspartate aminotransferase, alanine aminotransferase activities and increased levels of total bilirubin. Human miR-26a has been associated with numerous biological processes such as cell proliferation, apoptosis, tumorigenesis at different stages of non-tumor diseases, growth and development of normal tissues, and other biological processes, as reviewed (Gao and Liu, 2011;Fukumoto et al., 2015;Kato et al., 2015;Miyamoto et al., 2016). Increased expression of miR-26b-5p was associated with radiation exposure in a post-Chernobyl breast cancer epidemiology study (Wilke et al., 2018). MiR-26a has similar functions in humans and mice, among them inhibiting TGF-β (Koga et al., 2015). Thus, miR-26a-5p has tumor suppressor activities. MiR-26a has also been observed to be involved in the regulation of pro-inflammatory cytokines in microglia. Overexpression of miR-26a decreased the production of inflammatory cytokines (TNF α, IL-6) and the activating transcription factor (ATF) 2 was directly targeted by miR-26a in mice (Kumar et al., 2015). MiR-26a-5p is thus a very likely candidate as a predictor and biomarker of radiation, as indicated by previous studies showing links to cancer, liver function and radiation.
An enigma is that Miyamoto and co-workers (Miyamoto et al., 2016) assigns the same role to miRNA-26a-5p and miR-26b-5p, while we observed that these two miRNAs have opposite estimates in the 2.5 mGy/h in the generalized regression analysis, and their logistic regressions also go in opposite directions (SFigure 3). However, this is not observed at the higher dose rates, where both predict in the same direction. One can speculate that miR-26a and 26b may have opposing functions below a certain dose rate threshold.
MiR-140-3p only appears as a predictor at 2.5 mGy/h in generalized regression, while it is one of the top five most significant miRNAs in the logistic regression for all three dose rates. PAM and random forest make miR-140 a favorite pick at most dose rates. Logistic regression (SFigure 3) reveals that the expression level is increasing with dose rate. It seems to be a good predictor for exposure in general in this study, and even as predictor of dose-response. MiR-140 has been seen to reduce radiation fibrosis after radiation therapy through inhibiting TGF-β (Duru et al., 2016). Overexpression of miR-140-5p suppressed cell proliferation/SOX4 expression/invasion in colorectal cancer. Similar roles have been seen in many other cancer developments, including hepatocarcinomas, some actions also related to the Wnt1 pathway (Fang et al., 2017;Ji et al., 2018;Li et al., 2018;Lu et al., 2017;Song et al., 2009;Takata et al., 2013;Wolfson et al., 2014;Wu et al., 2019). It seems that miR-140 and miR-26a have similar roles, further supporting that this miRNA could be a general predictor of radiation exposure independent of dose rate.
MiR-133a-3p appears as a predictor of dose rates ≥10 mGy/h with all prediction methods. Although its expression level was low in all samples, the contrast is slightly higher than controls already at 2.5 mGy/h, higher at 10 mGy/h and higher still at 100 mGy/h. It ranks as number two in the average contribution to the predictors across dose rates (Fig. 5B). This indicates that MiR-133a is another predictor of dose-response (SFigure 2). MiR-133a has been seen as a conserved radiation predictor in mice and monkeys (Tomasik et al., 2018), and is also associated with cancer (Jia et al., 2013), inflammation and oxidative stress Sturrock et al., 2014). One report (Roderburg et al., 2013) describes a negative link between TGF-β and miR-133a in human hepatocytes, and how it suppresses collagen production and may prevent fibrosis. MiR-133a serum levels were higher in patients with chronic liver disease. It is involved in development and function of smooth and striated muscle cells (Liao et al., 2013), linking it to muscular and cardiovascular impairment  which is known to be one important noncancer effect of radiation.
MiR-145a is involved in the development from stem cells to mature cells, particularly smooth muscle cells, therefore indicating involvement in cardiovascular disease (Turczynska et al., 2012). MiR-145 directly target TGFβ receptor II (TβRII) (Ishii et al., 2018). Although the TGFβ pathway has a complicated role in cell differentiation, it has also been associated with liver fibrosis (Men et al., 2017). MiR-188-3p is the miRNA with the second highest estimate in generalized regression at 2.5 mGy/h, and it also contributes highly in the other methods, including the logistic regression. It does not, however, contribute to the principal components for the overall response to radiation, thus it appears to be the most specific predictor of the lowest dose rate. MiR-188-3p has no long record of proposed activities, however, it is reported to lower macrophage inflammatory responses (Zhang et al., 2018) connected to neurological and cardiovascular degeneration . This may be related to radiation influence, since low dose radiation may lower inflammatory reactions, while higher doses stimulate inflammation. The dose of 3 Gy is not considered as a low dose and there may be anti-inflammatory actions triggered at lowest dose rate, while not at the higher dose rates (Arenas et al., 2012;Large et al., 2015).
In our analyses miR-122-3p was non-significantly upregulated in the 2.5 and 10 mGy/h-groups and downregulated in the 100 mGy/hgroup, whereas it was the most important downregulated miRNA in the X-ray exposed group (Table 3). MiR-122 is liver-specific, and suppress the expression of the paternally expressed gene 10 (PEG10), which is implicated in the primary liver malignancy hepatocellular carcinoma (HCC) (Shyu et al., 2016).
The small sample size in this study limits the conclusions about what miRNAs and which mechanisms are involved in the biological processes after γand X-ray exposures. We identify several significantly A) expressed miRNAs, and we can predict radiation exposure at several dose rates without misclassifications. However, no clear pattern of dose-response to increasing dose rate was found. On the contrary, the best predictions were at the lowest dose rate. This may be due to different miRNA response patterns at the higher dose rates and our panel works best on the lower dose rates. It is also clear that the prediction methods applied were unstable with this small sample, although the generalized regression methods gave consistent results most of the time. Another limitation of the study was that these 21 miRNAs were analyzed only in liver samples, and miRNAs usually exhibit tissuespecific expression patterns. As discussed above several of the 21 miRNAs are reported to be expressed in other organs than liver, and there is no reason why these 21-miRNA panel cannot serve as candidate markers of chronic low grade stress exposure in other organs. However, in future study, the predictive capability of these 21-miRNA panel following chronic low grade stress exposure has to be demonstrated in other organs. As far as we know the only liver-specific miRNA among the 21-miRNA panel is miR-122-3p.
Overall, the biological response of radiation also includes conditions such as inflammation and DNA damage responses that may be inflicted by other stressors than radiation itself. Future studies will reveal the specificity of the currently established miRNA predictor signatures.
The identified miRNA predictor panel in our study needs to be confirmed with a larger number of samples and other bio-specimens.

Conclusions
This study demonstrates an approach to establish miRNA-based predictor signatures for stressors, such as radiation, with power to predict also low chronic stressor scenarios, in this case chronic low dose rate ionizing radiation. Our overall finding is that by selecting a small subset of miRNAs (21 of 576 investigated miRNAs), one can predict whether a sample has been exposed to γ-radiation or not with high accuracy. This includes long-term exposure to a low dose rate γ-radiation. However, the nine different estimation methods used for identifying the predictors often present different combinations in their prediction models. Partition models like random forest and boosted tree require larger sample sizes to run the validation part with reproducible results, hence, only training set results are presented. The nearest shrunken centroid method is fairly consistent but does not come up with significant predictors. Generalized regression showed superior performance with our data and presented reproducible results with high R 2 , although they differ somewhat in the choice of predictors. Signature sets of four to six miRNAs distinguish between mice exposed to 2.5 mGy/h, B) 10 mGy/h or 100 mGy/h γ-radiation were identified. Because of the variability in the predictor selection, the safest variable reduction is to omit the least informative miRNAs in a selection of methods. However, prediction of X-ray exposure required the use of the full set of 21 predictors. Using lager samples, our miRNA-based signature identification approach may be extended to other stressors and biospecimens, as well as species such as humans, and serve as powerful bio-predictors of a wide variety of stressors such as environmental, nutritional and lifestyle mediated stressors, also in chronic low stress scenarios.

Abbreviations
Cq quantification cycle threshold qPCR quantitative real-time PCR miRNA microRNA

Authors' contributions
Conceived and designed the experiments: AG, DME and AKO. Animal experiments: AG, DME, DAB. Laboratory analysis: AG, MLA and ND: Analyzed the data: MLA, ND, DME and AKO. Contributed reagents/materials/analysis tools: ND and DAB. Wrote the paper: ND, DME, AKO and DAB.

Declaration of competing interest
The authors declare no competing interests.

Acknowledgements
We would like to thank Arip Ikhsani and Yetneberk Kassaye for excellent animal work and handling of the often very inconvenient exposure schedules in this experiment.

Availability of data and materials
All relevant data of this article are included within the article and its additional files.

Consent to publish
Not applicable. No individual personal identifiable information is used in this work.