High-resolution mass spectrometry-based non-targeted metabolomic discovery of disease and glucocorticoid biomarkers in an animal model of muscular dystrophy

Urine is increasingly being considered as a source of biomarker development in Duchenne Muscular Dystrophy (DMD), a severe, life-limiting disorder that affects approximately 1 in 4500 boys. In this study, we used the mdx mice—a murine model of DMD—to discover biomarkers of disease, as well as pharmacodynamic biomarkers responsive to prednisolone, commonly used to treat DMD. Longitudinal urine samples were analyzed from male age-matched mdx and wild-type mice randomized to prednisolone or vehicle. We used high-resolution mass spectrometry to discover metabolic biomarkers of both disease and glucocorticoid treatment. A large number of metabolites (869 out of 6,334) were found to be significantly different between mdx and wild-type mice at baseline (Bonferroni-adjusted p-value < 0.05), thus being associated with disease status. These included a peak with m/z=357 and creatine, which were also discovered in a previous human study looking at serum. Novel observations included biliverdin and hypusine. These four peaks were also significantly higher in mdx mice compared to wild-type, as well as significantly associated with time after the baseline. Creatine and biliverdin were also associated with treatment after the baseline, but the association with creatine may have been driven by an imbalance at baseline. In conclusion, our study reports a number of biomarkers, both known and novel, which may be related to either the mechanisms of muscle injury in DMD and/or prednisolone treatment.


Introduction
Oral corticosteroids are the standard of care in Duchenne Muscular Dystrophy (DMD)-an X-linked genetic disease that affects 1 in 4000 to 5000 boys worldwide [1,2] . DMD is caused by mutations in the dystrophin gene located on the short arm of the X chromosome [3,4] . Oral corticosteroid use in young boys with DMD changes the natural history of the disease. One of the largest natural history studies in DMD showed that corticosteroids increased survival in DMD, improved motor outcomes such as continued ambulation, as well as cardiopulmonary outcomes [5] . Despite these unambiguous clinical benefits, there is great variation in the dosage and regimen of corticosteroids, and with compliance to treatment in DMD [6] . Patients and families often discontinue corticosteroids mainly due to unpleasant side-effects including sleep disturbance, weight gain, increased fracture risk, and metabolic syndrome. Identification of non-invasive pharmacodynamic biomarkers associated with glucocorticoid side effects in boys may help predict those individuals who might be at a higher risk for developing corticosteroid-induced side-effects [7] .
We have previously identified a metabolic signature consisting of serum biomarkers in DMD [8] . Metabolites that were increased in individuals with DMD compared to healthy controls included creatine, arginine, and unknown compounds with m/z values of approximately 357 and 312 Da; metabolites that were decreased in DMD included creatinine, androgen derivatives, and other unknown, yet-to-be-identified metabolites. In order to better characterize these metabolites as potential biomarkers in DMD, we used the mdx mouse-a widely used pre-clinical model-to agnostically evaluate urine metabolites associated with dystrophin deficiency and/or response to corticosteroids. Urine is an easily accessible bodily fluid, can be collected non-invasively and on multiple time-points, and is becoming an attractive source of biomarker development in DMD [9][10][11][12] .
The identification of a palette of pharmacodynamic urine metabolites will allow the identification of metabolic pathways that are affected by the disease pathogenesis and that might be responsive to therapies such as corticosteroids and exon skipping. Such pharmacodynamic biomarkers may also help in drug development and dose finding.

Animals
The experiments using mice described in this paper were conducted according to the institutional guidelines regarding the humane treatment of animals. This study was done in compliance with the National Institutes of Health guidelines for pre-clinical studies as described by [13] and the standard operating procedure for pre-clinical studies [14] . The protocol was approved by the local Institutional Animal Care and Use Committee (#30424).
Age-matched 5-week-old male wild-type (C57BL/10ScSnJ) and mdx (C57BL/10ScSn-Dmd < mdx >/J) mice litters were purchased from Jackson Laboratory (Bar Harbor, ME). Mice were acclimatized for 4 days prior to commencement of experiments. A total of 48 mice (23 wild-type (WT) and 25 mdx ) were used in order to ensure a rigorous statistical analysis. WT and mdx mice matched for body weight were randomized into treatment with either vehicle control (cherry syrup) or prednisolone (dissolved in cherry syrup) for a total of 4 weeks as previously published [15] . Animals were treated with vehicle or prednisolone (5mg/kg/day) by mouth based on body weight. Animals had access to water and food ad libitum. The study team was blinded to the treatment groups. Unblinding of the different treatment groups occurred after data collection and data processing.

Urine collection
Urine from mice was collected by placing them in clean cages early in the morning. Urine was collected at three different time points: baseline (prior to treatment), 15 days following commencement of treatment, and at completion of treatment at 4 weeks; these are coded as times T0, T1, and T2 in the remainder of this manuscript. We note that at T0, the mice had not received any vehicle cherry syrup. The urine collection had the following caveats: 4 mdx and 2 WT mice died between times T0 and T1, and some mice did not urinate at every time point; one sample was not available for T2 and therefore not analyzed. In total, 48 mice were evaluated at time T0 and 43 unique mice at times T1 and T2, out of which 41 provided measurements at T1 but not T2 and 41 provided measurements at T2 but not T1.

LC-MS method
For metabolite extraction, 80 µl of a solution of 50% acetonitrile in water containing internal standards (10µl of 1mg/ml debrisoquine and 50µl of 1mg/ml 4-nitrobenzoic acid added to 10 ml of 50% acetonitrile in water) was added to 20 µl of each urine sample in an Eppendorf vial. The samples were centrifuged at 13,000 rpm for 20 minutes at 4°C and the supernatant transferred to fresh vials for UPLC-Qtof analysis. 2µl of each sample was injected onto a Waters Acquity BEH C18 1.7 µm, 2.1 × 50 mm column using an Acquity UPLC system by Waters Corporation, Milford, MA. The gradient mobile phase consisted of solvent A -100% water with 0.1% formic acid -and solvent B -100% acetonitrile with 0.1% formic acid. The column temperature was set to 40°C and flow rate to 0.5 ml/min. The gradient started with 95% of solvent A and shifted to 80% of solvent A at 4 minutes at a ramp of curve 6. Pooled quality controls samples were injected after every 10 injections.
Data processing LC-MS metabolomics samples were preprocessed as previously described [8] . In brief, the XCMS approach [16] , available through the xcms package on Bioconductor [17] was used to detect features and estimate intensities, perform retention time correction, and group peaks from different samples, followed by filling in missing peaks via integrating the signal in the peak region defined by the other samples. The resulting peaks were then annotated using the CAMERA package [18] , resulting in the prediction of possible isotopes.
Following these initial preprocessing steps, peak groups with fewer than 37 peaks were removed, along with predicted isotopes. This led to 3,674 peaks in the positive mode and 2,663 peaks in the negative mode.
Internal standard normalization was then performed, with the intensities in the each mode being divided by the corresponding standard and the standards being removed from the list of peaks, leading to a total number of 6,335 peaks. A small number of peak-sample combinations had intensities of 0 (51 for the negative mode and 3 for the positive mode); these intensities were replaced by the smallest non-zero intensity for that peak across all samples. The intensities were then log2-transformed and the peaks from both modes were quantile-normalized together [19] . All the quality control samples were removed from the downstream analysis.

Statistical analyses
Given that the mice had a different diet at time T0 compared to times T1 and T2, the analyses considered were: 1) a comparison of the genotypes at time T0 via a two-sample t-test and 2) a mixed effects linear model for times T1-T2 which included genotype, time, and treatment, as well as all their interactions, with a random intercept to account for the within-mouse correlation. We note that for analysis 1) we did not consider treatment given that all the mice were at baseline, with none being yet treated with prednisolone. Prior to analysis 1), a single peak that had the same value in all samples at T0 was removed, so that this analysis was performed on 6,334 peaks. The correlation between the p-values for t-tests assuming equal variances and t-tests assuming difference variances in the two groups was > 0.9999; moving forward, we used the results from the analysis with the equal variance assumption. For analysis 2), we focused on the top 50 most significant peaks from analysis 1), with the addition of creatine and creatinine, which are well-known to be important in DMD and were among the top metabolites associated with disease status in previous work [8] . In the mixed effects models, likelihood ratio tests were used to assess the association of genotype, time, and treatment with the transformed peak intensities. In order to account for multiple testing, statistical significance was determined to be reached if the Bonferroni-adjusted p-values were less than 0.05, when considering all 6,334 peaks. to the online available spectra from the databases Metlin [20] and HMDB [21] and manually. Further putative identification was performed as follows: The MZXML files were parsed using the pyteomics python library [22] .
For each sample, the retention time was used to access the scan of the targeted compound. The mzxml function in pyteomics was used to access the m/z and corresponding intensity information of the retention time.
The top 30 peaks with highest intensities were used for compound identification using MS/MS spectral libraries such as Metlin, HMDB, and MONA ( http://mona.fiehnlab.ucdavis.edu/ ). In addition to searching tools provided by these sources, we used MetaboSearch Pro [23] to search against these libraries and assign scores to the putative metabolite IDs.

Comparison of genotypes at baseline
Out of a total of 6,334 detected peaks, 868 were determined to show significant differences between the mdx and the WT mice (Bonferroni-adjusted p-value < 0.05). Given this large number, for the downstream analyses, we chose to focus on the top 50 most significant peaks, along with creatine (m/z = 132.08) and creatinine (m/z = 114.06), as described above. A subset of these 50 peaks was also selected for validation via MS/MS; the MS/MS spectra for the 11 selected peaks considered in this follow-up experiment are in Supplemental Figure   1. Creatine and creatinine are known to be important in DMD and were among the top peaks in our human serum study [8] . We note that the positive mode peak, with m/z = 357.25 (p = 5.11 x 10 -25 unadjusted), appears to be the same as the top peak in our previous work with DMD patients [8] . Interestingly, creatine was not among the top 50 peaks in the mouse model, but ranked 472 and was still deemed significant (p = 4.91 x 10 -8 unadjusted). Creatinine, however, was not significant (p = 0.98 unadjusted). Boxplots representing results for these 3 peaks for the comparison between genotype are shown in Figure 1a Table 2. For all peaks in Figure 1 besides creatinine, the intensities are significantly elevated in the urine of mdx mice compared to WT mice. In particular, for the peak with m/z = 357.25, this finding recapitulates our earlier findings from the serum of younger individuals with DMD. For creatine, older boys with DMD had higher values [8] . We note that the only match for this peak is a tripeptide, Pro Ile Gln; however, the MS/MS spectrum did not appear to have a perfect fragmentation pattern for a tripeptide perhaps due to overlapping contaminants. Upon matching this tripeptide against the Mus musculus proteome using only UniProt/SwissProt and not including isoforms at https://research.bioinformatics.udel.edu/peptidematch/index.jsp [25] , we obtained 770 proteins, including the mouse titin protein, in which it appears 7 times. Titin fragments have been previously detected in urine samples collected from DMD patients [10,[26][27][28] . Comparison of genotypes, time, and treatment for T1-T2 As discussed above, we focused our analysis at the later time points on the peaks found to have the lowest 50 p-values in the genotype comparison at T0, in addition to creatine and creatinine. This is because we are of the opinion that it is important to focus on peaks that are likely to be indicative of the disease process, with differences seen at an early age. Creatine and all the peaks selected for validation except for the one at m/z = 370.04 were among those associated with genotype at T1-T2. Creatine, biliverdin, and the peaks at m/z = 210.11, 484.16, 884.93, and 485.33 were significantly associated with treatment at T1-T2. Creatine and all the peaks selected for validation except for the one at m/z = 484.16 were among those associated with time.
Overall, 31 of these 52 peaks were significantly associated with genotype, 6 with treatment, and 17 with time.
The results from the HMDB query for the top peaks are in Supplemental Table 1. The overlaps between these categories are summarized in Figure 2. Figure 2. Overlaps between the peaks significantly associated with treatment, time, and genotype at times T1-T2 (Bonferroni-adjusted p-value < 0.05 for 6,334 peaks) out of the 52 peaks considered, consisting of the top 50 peaks most significantly associated with genotype at T0, plus creatine and creatinine. 19 peaks were not significantly associated with treatment, time, or genotype (bottom right corner).
In particular, we note that the peak with m/z = 357.25 was associated with genotype (p = 1.71 x 10 -40 and time (p = 1.99 x 10 -25 unadjusted), but not with treatment (p = 7.11 x 10 -4 unadjusted, p = 1 when adjusting for multiple comparisons). Creatine was associated with genotype, as well as with treatment and time. However, the treatment association for creatine may be driven at least in part by an imbalance at baseline between the wild-type mice later randomized to treatment 1 versus those randomized to treatment 2 (p = 0.009). Creatinine was not significantly associated with any variable. Figure 3 shows the trends for these five metabolites.
Hypusine and biliverdin were also strongly associated with genotype (p = 9.28 x 10 -9 unadjusted, p = 1.59 x 10 -16 respectively). Hypusine was also associated with time (p = 2.66 x 10 -8 unadjusted), but not with treatment (p = 0.16 unadjusted). Biliverdin was associated with both treatment (p = 2.24 x 10 -6 unadjusted) and time (p = 2.01 x 10 -14 ) and did not show a major imbalance at baseline. Once again, the results for the peak with m/z = 357.25, showing a decrease over time in the genotype with the disorder, recapitulate the results from the human serum and urine studies [8,10,29] . An example of a peak showing treatment differences (p = 6.87 x 10 -9 ) that did not have a large imbalance at baseline and that, interestingly, showed overall higher values in the mdx versus the WT group is the peak with m/z = 484.16, shown in Figure 4. Supplemental Figure 2 shows the baseline (T0) results for this peak.

Discussion
Our primary goal in this study was to identify urinary metabolite biomarkers associated with (i) muscle pathology and (ii) prednisone response in an animal model of DMD. We foresee that such information can guide future studies to evaluate biomarkers to predict disease progression and response to corticosteroids and other anti-inflammatory drugs. To our knowledge, this is the first study evaluating urine metabolites in the murine model of DMD. There is a growing body of literature that shows that urine is becoming an attractive biomarker in DMD [9][10][11][12]30] , including that changes in urinary biomarker signatures may inform regarding response to therapy [30] . It is also easily accessible and less invasive, especially to very young boys. injury in DMD. The best-known biomarkers of DMD is creatine kinase, which seeps from muscle into blood and is generally used to identify young boys with this disorder, being many times higher than in healthy individuals [31] . While this an appropriate diagnostic biomarker, we have shown that it decreases with disease progression and muscle loss [32] . We have also previously shown that creatine and creatinine are substantially different in the serum of DMD cases compared to controls, with creatine being higher in the cases and creatinine being lower; these metabolites are both in the creatine metabolism pathway, which includes creatine kinase [8] . Our current work shows an increase in creatine at baseline in the mdx mice over the WT mice, although it tends to decrease over time, whereas in the human serum study, it stayed constant or slightly increased over time in DMD cases and decreased in controls. However, creatinine was not significant in any of the analyses considered. Differences in these results could be due to the biospecimen (serum vs urine), model (human disease vs murine model), or type of study (observational natural history study vs randomized laboratory study).
Particularly interesting is the metabolite corresponding to the peak with m/z = 357.25 (positive mode). This metabolite was the top metabolite in our previous study [8] . Furthermore, it recapitulates the trend of being generally higher in cases than in controls, as well as decreasing over time in cases. While we have been unable to definitively identify it herein, we consider that could be an important key to the pathogenesis of DMD and thus worthy of further investigation. If this peak is indeed the tripeptide Pro Ile Gln and it is a titin fragment, this would confirm many prior findings of elevated titin levels in both the urine and serum of young DMD patients [10,29,33] , however this remains speculative at this point.
Two novel observations from our study concern biliverdin and hypusine, metabolites which were higher in mdx urine compared to WT mice. At times T1 and T2, both metabolites showed associations with genotype and time, with biliverdin also showing an association with treatment. Biliverdin is an end-product in the heme pathway. Heme oxygenase catabolizes heme to equimolar biliverdin, carbon monoxide, and ferrous iron.
Heme oxygenase is considered a regulator of oxidative stress and is anti-inflammatory. Second, its product carbon monoxide has vasorelaxant property. Recent studies have provided contradictory information on 13 whether heme oxygenase-1 is decreased [34,35] or increased [36,37] in the muscle of mdx mice compared to controls. We postulate that a pro-inflammatory state seen in mdx mice drives active breakdown of heme, resulting in an increase in biliverdin levels in urine [38,39] . Further investigation of this pathway could thus lead to both an improved understanding of the disease pathogenesis and a validated biomarker to monitor disease progression and response to treatment. Hypusine is an amino acid found exclusively in the eukaryotic translation factor 5A (eIF5A) family [40] . It is important in nonsense-mediated decay [41] but its role in muscle pathology has not yet been established. It is possible that the increase in hypusine is a reflection of ongoing extensive muscle regeneration in mdx mice, especially during early disease stages. It is also know that [42] eIF5A induces inflammatory cytokine cascade and the nitric oxide synthase iNOS and inhibition of eIF5A protects mice from developing diabetes. While the precise role of biliverdin and hypusine in vivo in skeletal muscle in mdx is unknown, we postulate that their levels in urine likely reflect inflammation and regeneration in dystrophin deficient skeletal muscle.