**Comparison of Serum Pharmacodynamic Biomarkers in Prednisone-Versus Deflazacort-Treated Duchenne Muscular Dystrophy Boys**

**Shefa Tawalbeh 1,2, Alison Samsel 2, Heather Gordish-Dressman 3, Yetrib Hathout 2,\*, CINRG-DNHS Investigators** † **and Utkarsh J. Dang 4,\***


Received: 7 August 2020; Accepted: 10 October 2020; Published: 12 October 2020

**Abstract:** Prednisone (Pred) and Deflazacort (Dfz) are commonly used glucocorticoids (GCs) for Duchenne muscular dystrophy (DMD) treatment and management. While GCs are known to delay the loss of ambulation and motor abilities, chronic use can result in onerous side effects, e.g., weight gain, growth stunting, loss of bone density, etc. Here, we use the CINRG Duchenne natural history study to gain insight into comparative safety of Pred versus Dfz treatment through GC-responsive pharmacodynamic (PD) biomarkers. Longitudinal trajectories of SOMAscan® protein data obtained on serum of DMD boys aged 4 to 10 (Pred: *n* = 7; Dfz: *n* = 8) were analyzed after accounting for age and time on treatment. Out of the pre-specified biomarkers, seventeen candidate proteins were differentially altered between the two drugs (*p* < 0.05). These include IGFBP-2 and AGER associated with diabetes complications, and MMP-3 associated with extracellular remodeling. As a follow-up, IGFBP-2, MMP-3, and IGF-I were quantified with an ELISA using a larger sample size of DMD biosamples (Dfz: *n* = 17, Pred: *n* = 12; up to 76 sera samples) over a longer treatment duration. MMP-3 and IGFBP-2 validated the SOMAscan® signal, however, IGF-I did not. This study identified GC-responsive biomarkers, some associated with safety, that highlight differential PD response between Dfz and Pred.

**Keywords:** Duchenne muscular dystrophy; pharmacodynamic biomarkers; prednisone; deflazacort; glucocorticoids; corticosteroids; safety

#### **1. Introduction**

Duchenne muscular dystrophy (DMD) is an X-linked recessive disorder affecting the expression of dystrophin protein, an essential protein that maintains muscle fiber integrity and function [1]. The lack of dystrophin expression leads to muscle inflammation at an early stage of the disease, followed by progressive muscle degeneration and wasting [2]. While there are promising gene therapy and exon-skipping treatments, some of which have received conditional approval from the EMA [3] and FDA [4–7], these are mutation-specific and only restore a partial amount of truncated dystrophin protein [8]. Hence, DMD patients will continue to need combination therapy using the current standard of care: glucocorticoids (GCs). GC use helps reduce muscle inflammation and delay the loss of motor

abilities [9], delay loss of independent ambulation, improve pulmonary function, and delay the onset of cardiomyopathy [10–12]. However, a GC regimen can have many side effects, such as weight gain, growth stunting, loss of bone density, hirsutism, Cushingoid features, osteoporosis, hypertension, diabetes, behavioral disturbances, and difficulty in sleeping [12–14]. Commonly used steroidal drugs for DMD are prednisone (Pred) and deflazacort (Dfz). Prednisone has been used to treat DMD patients in the USA since the seventies [15]. In February 2017, the FDA approved deflazacort to treat DMD patients aged five years and older [16] although this same drug has been widely used in Europe to treat DMD patients for years. Figure 1 shows structures of prednisone and deflazacort along with their respective active metabolites: prednisolone and 21-desacetyl deflazacort, respectively. Both drugs are given to patients in their prodrug forms, which are then metabolized to their active forms in the liver. The active metabolite prednisolone has a hydroxyl group at carbon 17, while the active 21-desacertyl deflazacort has an oxazoline structure at that same position (red arrow in Figure 1).

**Figure 1.** 2D structures of Prednisone and Deflazacort (top panels) and their active drugs: prednisolone and 21-desacetyl deflazacort. The similarities are clear. The red arrows point to the structural differences between the active drugs.

Many efforts have been undertaken to investigate the use of biomarkers in DMD [17–19] and the efficacy and safety of Pred versus Dfz [10–12,20,21]. Both drugs are known to be efficacious in prolonging ambulation [10,11,14,22]. For our purposes here, we provide a small review on safety comparisons from the literature. These studies were compared to clinical outcome measures like changes in height, weight, Cushingoid features, and erythema, among others to evaluate differences in safety. A comparison of high dose Dfz and Pred in 70 systemic lupus patients concluded that Pred was associated with a significant increase in weight gain, Cushingoid severity index, and hirsutism compared to Dfz [23]. In a double-blind, randomized study on 196 DMD subjects, Pred was found to be associated with more weight gain than Dfz [10]. In another randomized study (18 DMD patients), patients treated with Pred showed higher weight gain compared to Dfz-treated patients [13]. As compared to other studies [10,14], a comparison of prednisone/prednisolone versus Dfz treated patients using 340 participants from the Cooperative International Neuromuscular Research Group Duchenne Natural History Study (CINRG-DNHS) showed that participants treated with Dfz had increased frequency of Cushingoid appearance, cataracts, and growth delay. Another study [13] reported that the time of initiating, dosing, and duration of treatments were associated with side effects; longer duration and increased Dfz dosage predicted growth stunting and Dfz was reported to be associated with lighter weight and shorter heights compared to Pred [13]. To summarize, these studies suggest that

treatment with Pred is associated with relatively more weight gain, whereas Dfz treatment is associated with relatively more growth stunting. While these studies provide important comparisons between Dfz and Pred, the underlying molecular mechanism leading to these differences, remain unknown. Furthermore, predictive outcome measures, especially for adverse effects are highly desirable.

There are ongoing debates among families, clinicians, and regulatory agencies over which GC drug (Pred or Dfz) and regimen is better for DMD boys. GC responsive biomarkers (for both safety and efficacy) have been defined and confirmed in different diseases and multiple cohorts [17,24,25], using data with pooled corticosteroid drugs not differentiating between Pred and Dfz and different regimens [22]. However, much remains unknown about the comparative effects of Pred vs. Dfz on blood accessible biomarkers and how they can inform clinical decision-making and drug development. Adding pharmacodynamic (PD) biomarkers to the evidence from the aforementioned clinical outcome studies could bring insights into differences between Pred and Dfz at the molecular level. To the best of our knowledge, we are the first to compare serum accessible PD biomarkers in DMD patients (4–16 years old) treated with Dfz or Pred to gain insight into differences at the serum protein level. For this, we use data generated by a high-throughput SOMAScan® technique, followed by confirmation of a specific set of PD biomarkers by ELISAs. The objective here is to define PD biomarkers that differ in their response to Pred and Dfz, and investigate if these biomarkers are known to be associated with reported differences in term of side effects between these two drugs.

#### **2. Materials and Methods**

The study protocol was approved by Institutional Review Boards at all participating institutions that provided serum samples. These included, the Office of Research IRB administrations at the University of California Davis, Davis, CA, the University of Pittsburgh, Pittsburgh, PA, Children's National Health System, Washington DC, the Conjoint Health Research Ethics Board, University of Calgary and the Human Subject Research Review Committee at Binghamton University, NY. Informed written consent was obtained from the parents of the participants or their legal guardians for biomarker studies at each site.

#### *2.1. SOMAscan*® *Dataset*

The CINRG-DNHS is a natural history study of 440 DMD boys. The CINRG-DNHS investigators enrolled these subjects at 20 centers in nine countries [11]. In this prospective cohort study, subjects genetically confirmed to have DMD were followed for up to 10 years. Details about the study, including informed consent and others, have previously been published [26]. At entry, some DMD subjects were steroid-naïve and were then treated during follow up visits, some remained steroid-naïve throughout, while others were already on GCs [26]. Characteristics (height, weight, age, time on GC treatment, and type of GC) were recorded on follow-up visits as well as clinical outcome data on time to run/walk, time to stand, and time to climb until loss of ability [26]. Note that six-minute walk distance and the NorthStar Ambulatory Assessment score were measured for a small subset of patients enrolled later during the study [26]. Blood samples for biomarker investigations were collected by the CINRG-DNHS investigators during some visits on a subset of participants. GC-treated patients were on Pred or Dfz with variation in dose and dose schedules (intermittent or daily). In a recent study [17], we generated biomarker data on a subset of these DMD samples (*n* = 31 boys, 4 to 10 years old) and used steroid-naïve and steroid-treated paired samples to define 107 (false discovery rate-adjusted *p*-values < 0.05) PD biomarkers from 1310 serum proteins (SOMAscan® array).

Based on the above-mentioned study [17], we pre-specified the 107 PD biomarkers of interest and focused on an analysis of 35 longitudinal samples from 8 Pred-treated and 7 Dfz-treated ambulatory patients only (screening dataset A) to identify potential candidate serum proteins differentially altered between Pred and Dfz. These subjects were age-matched (see Table 1 for a summary of their characteristics). For these longitudinal data, only subjects with at least 2 visits were used. The name

of the study, type of treatment, duration of treatment, age, and some clinical measurements (height, weight, BMI) are tabulated in Supplementary Table S1.

**Table 1.** Characteristics of patients/samples for SOMAscan® screening dataset A. All samples were from ambulatory patients.


#### *2.2. Validation of Key Pharmacodynamic Biomarkers Using ELISA*

A larger data set (confirmation dataset B) contained up to 76 longitudinal samples on 29 ambulatory subjects (17 deflazacort-treated and 12 prednisone-treated) from CINRG DNHS. Note that all samples/ subjects from screening dataset A were used in confirmation dataset B as well. For both datasets, we focused on ambulatory patients to minimize any possible confounding issues with the stage of disease (after loss of ambulation). Again, these subjects were age-matched. Summary characteristics of the subjects are provided in Table 2.

**Table 2.** Characteristics of patients/samples for confirmation dataset B. These serum samples were used for ELISA confirmation of MMP3, IGFBP-2 and IGF-1. All the samples were from ambulatory patients.


A subset of PD serum protein biomarkers, including MMP3, IGFBP-2, and IGF-I, were selected for confirmation analysis using ELISA assays. An ELISA kit from Meso Scale Diagnostics (MSD) was used to measure levels of MMP3. Similarly, ELISA kits from Thermo Fisher Scientific and R&D Systems, Inc. were used to measure levels of IGFBP2 and IGF-I, respectively. These candidates were chosen based on the availability of a validated ELISA assay kit. All three assays were sandwich ELISA, which were highly specific for antigen detection. Dilution factors for serum samples were 1:10, 1:400, and 1:100 for MMP-3, IGFBP-2, and IGF-I, respectively. The three assays were performed following the manufacturer's instructions. Details about serum samples used in ELISA validation are shown in Table 2.

#### *2.3. Data Analysis and Statistical Methods*

As depicted in the schematic (Figure 2), we first pre-specified 107 steroidal PD biomarkers previously identified as being significantly altered over time in their levels between GC-naive and paired GC-treated DMD samples from the same subjects [17]. To compare the serum levels and trajectories of these pre-specified 107 PD biomarkers between DMD patients treated with Dfz and Pred, linear mixed effect models were used to analyze longitudinal measurements for 7 boys on Pred and 8 on Dfz (daily regimen). All statistical analyses were performed using R [27]. Linear mixed effect models were run using the lme4 and lmerTest packages [28,29]. Serum protein levels, which had previously been hybridization control and median signal normalized, were log-transformed (see [17] for details). Random intercept linear mixed effect models were used to investigate the association of (mean centered) time on drug, type of GC (Pred or Dfz), and interaction between time on drug and type of GC on log transformed protein RFUs. The samples were age-matched at baseline (average age of Pred-treated subject at baseline = 5.7 (min age = 4.3, max age = 7.3) years; average age of Dfz-treated subject at baseline = 6 (min age = 4.7, max age = 7.4) years) and over time (see Table 1 for age range at sample collection). Some models did not converge (7 out of 107) due to numerical optimization issues (likely due to small sample size). Here, the interaction coefficient represents the difference in estimated slopes of biomarker trajectory over time between Pred-treated and Dfz-treated DMD patients. If the interaction coefficient is significant, this indicates that treatment type (Pred or Dfz) is associated with different trajectories of biomarker over time. The coefficient for the type of GC represents the difference in mean protein levels between Pred-treated and Dfz-treated subjects for a patient with an average treatment duration with similar biomarker trajectories over time (this is considered after the interaction effect). The same model was also fit to log-transformed ELISA data. To study the effect of time on GC on weight, height, and BMI, data analyses were performed using linear mixed effect models on confirmation dataset B. For this, we used a random intercept linear mixed model adjusting for age at baseline, time on GC, and interaction between type of GC and time on GC. Concordance of SOMAscan® and ELISA measurements was investigated on the subset of samples common to screening dataset A and confirmation dataset B.

**Figure 2.** Schematic describing the workflow for the biomarker analyses. Proteinsidentified as differentially affected by prednisone vs. deflazacort were identified from the SomaScan® based screening dataset. This was followed by confirmation analysis using ELISA assays of three selected biomarkers on a larger set of subjects (confirmation dataset). n: number of subjects, Ns: number of serum samples used in this study.

#### **3. Results**

#### *3.1. Longitudinal Trajectory of Serum PD Biomarkers in Prednisone vs. Deflazacort Treated DMD boys*

We examined differences in the longitudinal trajectories of pre-specified 107 PD biomarkers in Pred-treated patients (*n* = 7) versus Dfz-treated DMD patients (*n* = 8) accounting for duration of treatment (time on GC) and interaction between treatment duration and type of GC (Pred or Dfz). In general, Pred and Dfz seem to engage the PD biomarkers similarly and altered their trajectory in the same direction in screening dataset A. This was not the case for 17 PD biomarkers that were found differentially altered in their average levels and/or longitudinal trajectories (unadjusted *p* value < 0.05) in Pred-treated vs. Dfz-treated DMD patients. These 17 differentially altered PD biomarkers are listed in Table 3 with fold change between non-GC-treated DMD subjects and healthy controls (from [17]), fold change between Pred and Dfz-treated subjects, p values for the difference in the mean levels and p values for the difference in longitudinal trajectories between Pred and Dfz treated DMD patients, along with potential significance and references to published literature.


*J. Pers. Med.* **2020** , *10*, 164

molecule L1; IGFBP-2 =

Pred-treated and Dfz-treated subjects, respectively.

Insulin-like growth

factor-binding

 protein 2. 2 The two p-values are for the difference in mean levels and difference in longitudinal

 trajectories between

Two major groups of differentially affected steroidal PD biomarkers were thus identified. The first major group consisted of PD biomarkers that were repressed by GC in general, but exhibited a significantly lower mean level in Pred treated group relative to Dfz treated group. These included leukocyte immunoglobulin-like receptor subfamily B member 1 (LILRB1), tumor necrosis factor receptor superfamily member 21 (TNFRSF21), chordin-like protein 1 (CHRDL1), soluble advanced glycosylation end product-specific receptor (sRAGE), annexin A2 (ANXA2), CD166 antigen (CD166), scavenger receptor cysteine-rich type 1 protein M130 (sCD163), induced myeloid leukemia cell differentiation protein (MCL-1), transmembrane glycoprotein NMB (GPNMB), mitogen-activated protein kinase 14 (MAPK14), and neural cell adhesion molecule L1 (NCAM-L1). The second (minor) group consisted of PD biomarkers that increased in their levels following GC treatment and consisted of two subgroups. The first subgroup consisted of those that were significantly elevated in serum samples of Dfz relative to the Pred-treated group such as hemojuvelin (HJV), stromelysin-1 (MMP3) while the second group consists of those biomarkers that were significantly elevated in serum samples of Pred relative to Dfz-treated group such as insulin-like growth factor I (IGF-I), ficolin-1 (FCN1), and cGMP-inhibited 3 ,5 -cyclic phosphodiesterase A (PDE3A). One interesting PD biomarker is insulin-like growth factor-binding protein 2 (IGFBP-2); this showed diverging longitudinal trajectories (*p* = 0.040) in the Dfz-treated group, as compared to the Pred-treated group. Similarly, longitudinal trajectories of Ficolin-1 (FCN1; *p* = 0.03) and Annexin A2 (ANXA2; 0.054) also seemingly diverged in addition to their differential levels between the two drugs. Figure 3 shows selected examples of these differentially affected steroidal PD biomarkers. For example, IGFBP-2 sharply decreased over time in Dfz-treated DMD patients, while it remained unchanged or slightly increased over time in Pred treated group (*p* = 0.040). While TNFRSF21 (also known as DR6) was not found to be differentially affected by the two different drugs in terms of longitudinal trajectory slopes but its mean sera levels was significantly lower in the Pred-treated group compared to the Dfz-treated group (*p* = 0.005). Similarly, mean levels of FCN-1 were lower in the Dfz-treated group relative to the Pred-treated group (*p* = 0.038), while MMP3 mean levels were significantly elevated in the Dfz-treated group relative to the Pred-treated group.

#### *3.2. E*ff*ect of Dfz and Pred on Height and Weight of DMD Boys*

We also investigated differences in growth stunting in Dfz treated patients relative to Pred treated patients in confirmation dataset B. DMD patients treated with Dfz had lower height growth rates (*p* = 0.006 for difference in trajectory slopes over time) compared to those treated with Pred (Figure 5). We did not find any significant difference in weight (*p* = 0.112) and BMI (*p* = 0.08) over time between patients treated with Dfz and Pred (Figure 5). Additionally, note that more variation is observed in longitudinal BMI measurements around the estimated model (Figure 5).

#### *3.3. Data Validation Using ELISA*

For the reported data above, we had a small sample size and we did not adjust for multiple testing, thus a follow-up confirmation analyses on a subset of PD biomarkers was carried out. Unfortunately, from the list of PD biomarkers identified in Table 3 above, a good validated ELISA assay that used low sera volume was available for only three PD biomarker candidates: MMP3, IGFBP2 and IGF-I. There were other ELISA assays for other candidates, but they were either not validated or required a larger volume of sera samples, which we did not have. To confirm the SomaScan® findings obtained for MMP3, IGFBP2, and IGF-1, we used a larger sample size of samples over a longer treatment duration and a wider age range of subjects. The findings from the ELISA runs are summarized in Table 4. Results for MMP-3 validated the SOMAscan® signal (same directionality; p-value for the difference in mean levels = 0.022), IGFBP-2 neared significance (p-value for difference in trajectory slopes = 0.051), while IGF-I did not validate the SOMAscan® signal. Figure 4 shows correlations between SOMAscan® data and ELISA data with the IGF-I measurements having a substantially worse Pearson correlation than MMP-3 and IGFBP-2.

**Figure 3.** Selected examples of serum protein PD biomarkers showing difference between prednisone (Pred) and deflazacort (Dfz)-treated samples. FCN-1 has relatively higher mean RFU levels in Dfztreated patients (*p* = 0.038 for mean levels; *p* = 0.03 for difference in trajectory slopes over time). TNFRSF21 has higher mean RFU levels in Pred- vs. Dfz-treated DMD patients (*p* = 0.005 mean levels; *p* = 0.934 for difference in trajectory slopes over time). The longitudinal trajectories of IGFBP-2 levels are different between the two drugs (*p* = 0.04). MMP-3 mean RFU level is elevated in DMD boys treated with Dfz, as compared to Pred (*p* = 0.008).

**Figure 4.** SOMAscan® signal confirmation using ELISA. Upper panel shows longitudinal trajectories of selected DMD serum protein PD biomarkers from ELISA assays. Lower panel shows correlation plots between SOMA and ELISA data for the samples that overlapped between the screening and confirmation data sets.

**Figure 5.** Comparison of longitudinal trajectories between DMD patients treated with deflazacort and prednisone of BMI (*p* = 0.08 for difference in trajectory slopes), height (cm; *p* = 0.006), and weight (kg; *p* = 0.112).


**Table 4.** Summary for biomarker signal confirmation using ELISA.

<sup>1</sup> Two *p*-values are provided for difference in mean levels and difference in longitudinal trajectory slopes between Pred-treated and Dfz-treated subjects for both SOMAscan® and ELISA® analysis, respectively.

#### **4. Discussion**

GCs are and continue to be the standard of care for several chronic inflammatory and auto-immune diseases including DMD [11,21,25]. In this study, and due to the interest and ongoing debates between families and clinicians regarding which GC drug is more beneficial and safe for DMD boys, we focused on a subset of Dfz- an Pred-treated DMD patients to define the differential pharmacodynamic response to these two commonly prescribed drugs using blood circulating proteins. We previously defined a set of PD biomarkers that are responsive to GC treatment in DMD [17]. In that previous study, we showed that GC affected the levels of 107 circulating serum proteins. In general, use of GCs decreased the levels of several circulating pro-inflammatory and immune response associated proteins, but caused an increase in certain proteins associated with metabolism and extracellular remodeling [17]. In this current study, we find that, in general, both Dfz and Pred engaged the PD biomarkers in a similar fashion, however, among the 107 PD biomarkers, 17 exhibited a differential longitudinal response to Dfz relative to Pred, in directionality, mean levels, or both.

A close examination of the 17 differentially altered PD biomarkers (from screening dataset A) between Dfz and Pred in term of their levels, longitudinal directionality and their potential physiological function led to the following interpretation. Among these 17 PD biomarkers, 12 were significantly reduced in the Pred-treated group relative to the Dfz-treated group after adjusting for treatment duration (the samples were age matched). These 12 decreased PD biomarkers by Pred can be classified into subgroups with the first subgroup consisting of inflammatory and immune associated proteins such as LILRB1, TNFRSF21, sRAGE, CD166, and sCD163. Previous studies have suggested that Dfz is a stronger immune suppressant than Pred [43]. In contrast, our analysis showed that for the above-mentioned five serum proteins, Pred seemed to be more immunosuppressive than Dfz. However, further studies using a larger sample size and additional cellular biomarkers are needed to claim this finding. Another subgroup that was differentially decreased by Pred relative to Dfz included bone mineralization protein GPNMB and the cell adhesion protein NCAM-L1. GPNMB, also known as osteoactvin in rats, has been shown to be implicated in bone mineralization and bone regeneration [39,40] and the larger decrease in its mean level in the Pred treated group relative to the Dfz treated group could suggest that Pred treated patients might have a higher risk of losing bone density than Dfz treated patients. This agrees with earlier studies showing that decreases in bone density was markedly spared by Dfz, as compared to Pred in both adult and children cohorts [44]. However, further studies correlating low levels of circulating GPNMB to loss of bone density in Pred treated group relative to Dfz treated group are needed to support this hypothesis. The differential decrease in the circulating levels of NCAM-L1 by Pred could be indicative of another potential side effect of Pred over Dfz. Indeed, a recent study linked low levels of circulating NCAM-L1 to risk of developing type 2 diabetes [41], which is in agreement with an earlier study showing that Pred was more diabetogenic than Dfz in a pediatric population [45]. The differential decrease of sRAGE by Pred over Dfz could also be linked to the diabetogenic effect of Pred since sRAGE has been suggested to act as a decoy that dampens the advanced glycation end product signaling (e.g., RAGE signaling) and thus increases the risk of developing diabetes [33].

However, the differential decrease in the mean levels of MCL-1, MAPK14, and CHRDL1 by Pred relative to Dfz could be considered to be evidence of possible efficacy in DMD. Indeed, MCL-1 and MAPK14 were previously reported by us [17] to be elevated in blood of untreated DMD boys relative to age matched healthy controls then decreased by GC treatment toward the levels in healthy controls. CHRDL1, however, is known to bind to BMP4 and antagonizes its function [30]. BMP4 is a member of the member of the transforming growth factor-β (TGF-β), a well-known pathway reported to be involved in DMD pathogenesis [46]. Thus, a differential decrease of CRDL1 by Pred relative to Dfz could suggest a beneficial effect but further experiments using a CRDL1 knockout animal model are needed to verify this hypothesis.

Another relevant PD biomarker that was differentially altered over time between the two drugs is IGFBP-2. It sharply decreased in its longitudinal trajectory in Dfz-treated group compared to the Pred-treated group. This observation was further validated by ELISA assay (*p* = 0.0507; near significance). Note that, in our ELISA confirmations, we had the smallest number of additional samples available for IGFBP-2. IGFBP-2 binds to IGF-1 and regulates growth. While IGFBP-2 was decreased by Dfz treatment, IGF-1 was relatively increased by Dfz over Pred. Unfortunately, we were unable to validate the IGF-1 SomaScan® signal using ELISA; nevertheless, an increase in IGF-1 by GC treatment might be associated with antiinflammation efficacy [31]. However, the selective longitudinal decrease of IGFBP2 by Dfz could be associated with significant growth stunting by Dfz over the Pred reported by others [20,47,48] and confirmed by us in this study. Furthermore, a gene knockout study, conducted on mouse model comparing Igfbp2-/- mice colony with Igfbp2+/+ mice colony determined the role of IGFBP2 in bone turnover and showed that Igfbp2-/- males had shorter femurs and were heavier than controls but were not insulin resistant [42]. The decrease in IGFBP2 could have dual side effects on both growth stunting and risk of bone fracture. Indeed, a recent study compared a cohort of DMD boys treated with Dfz to a cohort treated with hydrocortisone and concluded that DMD boys receiving daily dose of Dfz had a higher incidence of bone fractures with greater risk of growth stunting [47].

The second group of steroidal PD biomarkers that were differentially affected by the two drugs are those that increased following treatment with GC. FCN1, IGF-I and PDE3A were relatively increased by Pred over Dfz, while MMP3 and HJV were relatively increased by Dfz over Pred. While FCN1 mean levels were higher in the Pred treated group relative to Dfz treated group, it decreased over time with Pred treatment while remaining relatively unchanged over time following Dfz treatment. This is intriguing and requires further validation using orthogonal methods such as ELISA. Unfortunately, the available ELISA assay required the use of a larger sample volume, which was a limitation of our study.

A comparison of the relative effects of Pred and dexamethasone (Dex), another glucocorticoid, on short-term growth and bone turnover confirmed that both drugs affected short-term bone turnover and growth [20]. However, Dex may be more potent in suppressing linear growth, simulating weight gain and bone turnover compared to prednisone. Dex was more potent at depressing IGF-I levels than prednisolone. This is consistent with our SOMAscan® data, IGF-1 was significantly more depressed or decreased in Dfz treated group relative to Pred treated group. However, we were not able to confirm this finding using ELISA. The discrepancy between the SOMAscan® data and ELISA data could be associated with an epitope effect. The two techniques recognize different epitopes on the IGF-I. Furthermore, IGF-I might co-exist as free form and bound form to circulating insulin growth factor binding proteins and these might further interfere with both these affinity-based assays that relies on epitope binding. Future analyses using denatured condition and targeted mass spectrometry analysis are needed to examine the true levels of circulating IGF-I.

MMP3, also known as stromelysin-1, was dramatically increased in the blood circulation following GC treatment, as previously shown by us [17] and others in both DMD patients [24] and inflammatory bowel disease patients [25] and anti-neutrophil cytoplasmic antibody-associated vasculitis, and juvenile dermatomyositis [25]. In this study, we further show that this increase is associated more with Dfz use

than Pred use. This was further confirmed by an ELISA assay on a confirmation dataset using a larger sample size with longer treatment duration. The mechanism by which GC induces levels of circulating MMP3 remains unclear. Further studies using cell culture expressing MMP3 are needed to explain the difference of action of Dfz and Pred on the expression of MMP3. Owing to the function of MMP3, i.e., degradation of extracellular matrix, this differential increase might result in more extracellular matrix remodeling in Dfz-treated DMD patients relative to Pred-treated patients. Whether this extracellular matrix remodeling is adverse or beneficial remain to be carefully examined.

#### **5. Conclusions**

Here, we investigated less invasive and objective PD biomarkers that might prove useful to monitor disease progression and response to GC therapies in DMD. We identified differences in serum levels of PD biomarkers between Pred- and Dfz-treated subjects, some of which may be associated with safety. Such blood accessible biomarkers in DMD could play an important role in clinical trials and decision making [49]. In our comparisons, we adjusted for duration of treatment, the comparison groups were age matched, and for a subset of proteins, we conducted ELISA validation testing of SOMAscan® signal. IGFB2 had a decreasing longitudinal trend associated with Dfz, at odds with the observation for Pred. This may be associated with differential growth stunting seen between the two drugs. The dramatic increase of MMP3 by Dfz relative to Pred remains to be interpreted. Further studies are needed to test the physiological significance of these Dfz and Pred differentially affected PD biomarkers. The study's limitations include the small sample size, no adjustments for multiple testing, and that the data come from a natural history study with lots of observed variability. We have tried to overcome these limitations by performing lab validation using ELISA and using a larger sample size and longer treatment duration. For biological validation, an external cohort is needed to validate our findings. Note also that we did not investigate direct associations of biomarkers with efficacy/safety outcomes; larger sample sizes are required for such studies and are the objective of ongoing research. While preliminary, this study identified serum proteins that were altered between the Dfz and Pred groups using SOMAscan® array data and we successfully validated two biomarkers using ELISA that may be associated with adverse effects. However, further studies using larger sample sizes collected from well controlled cohorts enrolled in ongoing and future clinical trials comparing the safety and efficacy of Dfz versus Pred are needed to validate and test these differentially altered pharmacodynamic biomarkers identified herein.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2075-4426/10/4/164/s1, Supplemental Table S1: A subset of DMD patients from the Cooperative International Neuromuscular Research Group Duchenne Natural History Study (CINRG-DNHS).

**Author Contributions:** Conceptualization, Y.H. and U.J.D.; Methodology, Y.H., U.J.D. and S.T.; Validation, S.T. and A.S.; Formal Analysis, S.T.; Data Curation, Y.H., U.J.D., H.G.-D. and S.T.; Writing—Original Draft Preparation, S.T.; Writing—Review & Editing, Y.H. and U.J.D.; Funding Acquisition, Y.H. The CINRG DNHS study was run by the CINRG DNHS investigators who also collected the bio-samples. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by DoD grant W81XWH-16-1-0557, NIH grant P50HD090254 and the Muscular Dystrophy Association USA (MDA353094).

**Acknowledgments:** The authors would like to thank the patients and their families for their participation in the study. The authors would also like to thank study team members for our participating CINRG sites: University of California Davis: C. McDonald, E. Henricson, M. Cregan, L. Johnson, J. Han, N. Joyce, A. Nicorici, D. Reddy; Sundaram Medical Foundation and Apollo Children's Hospital, Alberta Children's Hospital, Calgary: J. Mah, A. Chiu, T. Haig, M. Harris, M. Kornelsen, N. Rincon, K. Sanchez, L. Walker; Queen Silvia Children's Hospital: M. Tulinius, A. Alhander, A. Ekstrom, A. Gustafsson, A. Kroksmark, U. Sterky, L. Wahlgren; Children's National Health System: R. Leshner, N. Brody, B. Drogo, M. Leach, C. Tesi-Rocha, M. Birkmeier, B. Tadese, A. Toles and M. Thangarajh; Royal Children's Hospital: A. Kornberg, K. Carroll, K. DeValle, R. Kennedy, V. Rodriguez, D. Villano; Hadassah Hebrew University Hospital: Y. Nevo, R. Adani, A. Bar Leve, L. Chen-Joseph, M. Daana, V. Panteleyev-Yitshak, E. Simchovitz, D. Yaffe; Instituto de Neurosciencias Fundacion Favaloro: L. Andreone, F. Bonaudo, J. Corderi, L. Levi, L. Mesa, P. Marco; Children's Hospital of Pittsburgh of UPMC and the University of Pittsburgh: P. Clemens, H. Abdel-Hamid, R. Bendixen, C. Bise, A. Craig, K. Karnavas, C. Matthews, G. Niizawa, A. Smith, J. Weimer; Washington University: J. Anger, T. Christenson, J. Florence, R. Gadeken, P. Golumbak, B. Malkus, A. Pestronk, R. Renna, J. Schierbecker, C. Seiner, C. Wulf; Children's Hospital of Richmond at VCU: J. Teasley, S. Blair, B. Grillo, E. Monasterio; University of Tennessee: T. Bertorini, M. Barrett-Adair, C. Benzel, K. Carter, J. Clift, B. Gatlin, R. Henegar, J. Holloway, M. Igarashi, F. Kiphut, A. Parker, A. Phillips, R. Young; Children's Hospital of Westmead: K. North, K. Cornett, N. Gabriel, M. Harman, C. Miller, K. Rose, S. Wicks; University of Alberta: H. Kolski, L. Chen, C. Kennedy; Centro Clinico Nemo: M. Beneggi, L. Capone, A. Molteni, V. Morettini; Texas Children's Hospital: T. Lotze, A. Gupta, A. Knight, B. Lott, R. McNeil, G. Orozco, R. Schlosser; University of Minnesota: G. Chambers, J. Day, J. Dalton, A. Erickson, M. Margolis, J. Marsh, C. Naughton; Mayo Clinic: K. Coleman-Wood, A. Hoffman, W. Korn-Petersen, N. Kuntz; University of Puerto Rico: B. Deliz, S. Espada, P. Fuste, C. Luciano, J. Torres; and the CINRG Coordinating Center: Lauren Morgenroth, M. Ahmed, A. Arrieta, N. Bartley, T. Brown-Caines, C. Carty, T. Duong, J. Feng, F. Hu, L. Hunegs, Z. Sund, W. Tang, and A. Zimmerman. The CINRG-DNHS was funded by the U.S. Department of Education/NIDRR (#H133B031118, #H133B090001); U.S. Department of Defense (#W81XWH-12-1-0417); National Institutes of Health/NIAMS (#R01AR061875); Parent Project Muscular Dystrophy.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Multi-Omics Identifies Circulating miRNA and Protein Biomarkers for Facioscapulohumeral Dystrophy**

**Christopher R. Heier 1,\*, Aiping Zhang 2, Nhu Y Nguyen 2, Christopher B. Tully 2, Aswini Panigrahi 2, Heather Gordish-Dressman 1,2, Sachchida Nand Pandey 2, Michela Guglieri 3, Monique M. Ryan 4, Paula R. Clemens 5, Mathula Thangarajh 6, Richard Webster 7, Edward C. Smith 8, Anne M. Connolly 9, Craig M. McDonald 10, Peter Karachunski 11, Mar Tulinius 12, Amy Harper 13, Jean K. Mah 14, Alyson A. Fiorillo 1,2, Yi-Wen Chen 2,\* and Cooperative International Neuromuscular Research Group (CINRG) Investigators** †


Received: 25 September 2020; Accepted: 6 November 2020; Published: 19 November 2020

**Abstract:** The development of therapeutics for muscle diseases such as facioscapulohumeral dystrophy (FSHD) is impeded by a lack of objective, minimally invasive biomarkers. Here we identify circulating miRNAs and proteins that are dysregulated in early-onset FSHD patients to develop blood-based molecular biomarkers. Plasma samples from clinically characterized individuals with early-onset FSHD provide a discovery group and are compared to healthy control volunteers. Low-density quantitative polymerase chain reaction (PCR)-based arrays identify 19 candidate miRNAs, while mass spectrometry proteomic analysis identifies 13 candidate proteins. Bioinformatic analysis of chromatin immunoprecipitation (ChIP)-seq data shows that the FSHD-dysregulated DUX4 transcription factor

binds to regulatory regions of several candidate miRNAs. This panel of miRNAs also shows ChIP signatures consistent with regulation by additional transcription factors which are up-regulated in FSHD (FOS, EGR1, MYC, and YY1). Validation studies in a separate group of patients with FSHD show consistent up-regulation of miR-100, miR-103, miR-146b, miR-29b, miR-34a, miR-454, miR-505, and miR-576. An increase in the expression of S100A8 protein, an inflammatory regulatory factor and subunit of calprotectin, is validated by Enzyme-Linked Immunosorbent Assay (ELISA). Bioinformatic analyses of proteomics and miRNA data further support a model of calprotectin and toll-like receptor 4 (TLR4) pathway dysregulation in FSHD. Moving forward, this panel of miRNAs, along with S100A8 and calprotectin, merit further investigation as monitoring and pharmacodynamic biomarkers for FSHD.

**Keywords:** FSHD; biomarkers; miRNA; proteomics; calprotectin; dystrophy; muscle

#### **1. Introduction**

Facioscapulohumeral muscular dystrophy (FSHD) is an autosomal dominant muscle disorder with no current therapy, a variable prognosis, and complex genetic and molecular mechanisms. FSHD is caused by aberrant expression of *double homeobox 4* (*DUX4*) due to epigenetic changes of the *D4Z4* repeat region at chromosome 4q35 [1–3]. Roughly 95% of patients have Type 1 FSHD (FSHD1) due to contraction of the D4Z4 array; a small portion (~5%) of patients have Type 2 FSHD (FSHD2) caused by mutations in *the structural maintenance of chromosomes flexible hinge domain containing 1* (*SMCHD1*) gene, the *DNA methyltransferase 3B (DNMT3B)* gene, or the *ligand-dependent nuclear receptor-interacting factor 1 (LRIF1)* gene [4–6]. The aberrant expression of DUX4 protein causes mis-regulation of genes involved in germline function, oxidative stress responses, myogenesis, post-transcriptional regulation, and additional cellular functions [7–13]. These downstream molecular changes are believed to cause FSHD, although the exact mechanisms are not clear.

Although the onset of FSHD is generally around adolescent years, a small portion (~4%) of patients present with an early-onset or infantile form of FSHD [14]. Previous studies have shown that the disease severity of FSHD1 is negatively correlated with the size of *D4Z4* repeats [15,16]. Individuals with early-onset FSHD1 tend to have smaller *D4Z4* repeats and more severe disease phenotypes, including more profound muscle weakness, younger age at loss of independent ambulation, and extramuscular manifestations such as retinal vasculopathy or hearing loss [14,15,17,18].

In clinical practice, particularly with pediatric-onset FSHD, there is a low use of serial histological assessments because they require painful biopsies of muscle tissue that typically reveal patchy or uneven pathology. Given this, many patients no longer undergo muscle biopsy once a genetic diagnosis is made. Functional motor scales provide a non-invasive alternative to study neuromuscular disease progression; however, they can show great variability, can be age- or disease stage-limited, and they can be subject to placebo or coaching effects in clinical trials [19,20]. Circulating molecular biomarkers provide a promising alternative to these clinical assessments because they are objective measurements that can be assayed repeatedly over time using minimally invasive methods. Blood-based miRNAs or proteins that measure the progression of disease or a patient response to therapy over time are known as a monitoring biomarker [21]. In clinical trials, monitoring biomarkers may also be used as pharmacodynamic biomarkers to identify patients who are early responders to therapy, to demonstrate exposure-response relationships, or to improve statistical power and modeling. As patient populations are sensitive and limited for this relatively rare pediatric disease, less invasive monitoring or pharmacodynamic biomarkers are important for early-onset FSHD, as frequent serial biopsies are especially problematic in this population.

Recently, circulating miRNAs have emerged as exciting potential diagnostic, prognostic, and drug-responsive biomarkers. This is a class of small non-coding ribonucleic acid (RNA) molecules

(~22 nucleotides in length) that can help to regulate gene expression [22], and which are highly stable in biofluids such as blood and urine [23,24]. In rare diseases with highly variable symptoms, such as multiple acyl-coenzyme A dehydrogenase deficiency (MADD), the serum-based detection of muscle-specific miRNAs termed myomiRs can signal the presence of underlying muscle-specific pathologies [25]. In Duchenne and Becker muscular dystrophies, myomiRs are up-regulated in serum from both patient populations, while detection of miR-206 up-regulation can be used to differentially diagnose severe Duchenne versus Becker patients [26–28]. In addition to myomiRs, inflammatory miRNAs such as miR-146a, miR-146b, miR-221 and miR-155 have been found to be dysregulated in multiple forms of muscular dystrophies [29–31]. These two classes of miRNA show potential as pharmacodynamic biomarkers, with myomiRs proposed for muscle-stabilizing treatments such as gene therapy [32,33], and inflammatory microRNAs proposed for current steroids [34,35] as well as newly emerging dissociative anti-inflammatory drugs such as vamorolone [36–38] or edasalonexent [39,40]. In parallel to development of miRNA monitoring biomarkers, new advances in whole exome sequencing are enabling clinicians to diagnose novel mutations in over 60 genes known to be responsible for muscular dystrophies such as FSHD and limb-girdle muscular dystrophy (LGMD) [41–44]. Together, these advances will help to improve the diagnosis, monitoring, and treatment of a diverse number of diseases affecting muscle.

The development of circulating biomarkers for FSHD has the potential to improve clinical management and to facilitate the development of new treatments. In this study, we test plasma samples from a cohort of individuals with early-onset FSHD1 using both miRNA and proteomic profiling approaches. Our goal is to identify molecules that can be used to monitor FSHD disease activity and that may ultimately facilitate future therapeutic trials. Initial analysis of a discovery group identifies a panel of miRNAs and proteins as biomarker candidates. Bioinformatic analyses of ChIP-seq data provide a rationale for the changes in candidate biomarkers, as their behavior is consistent with changes in transcription factor pathways that are disrupted in FSHD1. Subsequent characterization in separate, non-overlapping groups of FSHD1 patients provides validation of nine biomarkers whose expression can be conveniently assayed by qRT-PCR or Enzyme-Linked Immunosorbent Assay (ELISA), and are increased in early-onset FSHD.

#### **2. Materials and Methods**

#### *2.1. Ethics Statement*

We obtained institutional ethics and research review boards approval for these clinical studies from the Institutional Review Board of Children's National Hospital and at all participating Cooperative International Neuromuscular Research Group (CINRG) sites, in accordance with all requirements, as previously described in Mah et al. 2018 [45]. Written informed consent was obtained from all the participants before the study procedures. Where applicable, informed consent and/or assent was obtained from all patients or legal guardians before enrollment.

#### *2.2. Patients and Sample Collection*

Plasma samples were collected and biobanked from a previous early-onset FSHD study conducted by CINRG as described by Mah et al. [45]. For the discovery experiments, FSHD1 patients aged 10 to 51 years old were included (*n* = 16 for miRNA discovery, *n* = 25 for proteomics discovery), along with healthy control volunteers (*n* = 8 for miRNA discovery, *n* = 17 for proteomics discovery) aged 16 to 54 years old. All patients had Type 1 FSHD caused by epigenetic changes due to *D4Z4* contraction which results in up-regulation of *DUX4*.

#### *2.3. miRNA Profiling*

RNA was isolated and quantified from the discovery cohort of patients as described previously [34]. Briefly, RNA was isolated from 150 μL of plasma using Trizol liquid sample (LS) reagent (ThermoFisher,

Waltham, MA, USA), then converted to cDNA using the High Capacity Reverse Transcription Kit with multiplexed reverse transcription (RT) primers (ThermoFisher). Synthesized cDNA was then pre-amplified using PreAmp MasterMix with multiplexed TaqMan (TM) primers corresponding to the RT primers used in initial cDNA reaction. Quantitative analysis of miRNA was performed via TaqMan Low-Density Array Cards (TaqMan™ Array Human MicroRNA A Cards v2.0; ThermoFisher). The ThermoFisher Cloud software suite with the Relative quantification (Rq) application was used to perform statistical analysis and determine expression of miRNA in either mild or severe FSHD1 patient groups versus healthy controls. A value > 1 indicates an increase and a value < 1 indicates a decrease in miRNA expression in FSHD1 versus healthy controls, with *p*-values ≤ 0.05 considered significant. To reduce false-positive discovery in this setting, we used an evidence-based approach where candidate miRNAs that significantly increased in the discovery groups were cross-referenced to a separate set of non-overlapping CINRG patients used as a validation group.

#### *2.4. Bioinformatics of miRNA Regulation via DUX4 and FSHD-Associated Factors*

Surrounding DNA regulatory regions of candidate miRNA genes were queried in ChIP-seq datasets for binding by transcription factors known to be impacted by FSHD. These analyses were performed using the University of California Santa Cruz (UCSC) Genome Browser with alignment to the GECh37/hg19 genome build. For primary effects, due to the underlying mutation that causes FSHD, DUX4 binding was queried. For this, we uploaded a user-supplied DUX4 ChIP-seq track published by Geng et al. [9] to determine which candidate miRNAs displayed physical binding of DUX4 at potential regulatory regions within 100 kb of the gene for each miRNA.

To investigate secondary factors whose dysregulation is associated with FSHD-causing mutations, we investigated DNA binding by transcription factors shown to be significantly up-regulated in cultured human muscle cells using microarray data by Geng et al. [9]. For this, we used ChIP-seq data from the Encyclopedia of DNA Elements (ENCODE) [46,47]. From a master list of DUX4-regulated genes published in [9], we identified a list of 34 transcription factors with ChIP-seq data from ENCODE available within the UCSC Txn Factor ChIP Track and 47 transcription factors from the Txn Factor ChIP E3 Track [48–50]. After an initial survey of these full transcription factor lists for the 19 candidate miRNAs, we narrowed down to a shorter focus list of 9 transcription factors whose binding was most frequently associated with the candidate miRNAs. DNA binding by transcription factors was queried in datasets produced using ChIP-seq from all 9 available cell line tracks, including GM12878 (lymphoblasts), H1-hESC (embryonic stem cells), HeLa-S3 (cervical cancer cells), HepG2 (liver cancer cells), HSMM (skeletal muscle myoblasts), HUVEC (umbilical vein endothelial cells), K562 (immortalized myelogenous leukemia cells), NHEK (epidermal keratinocytes), and normal human lung fibroblasts (NHLF).

In addition to binding by DUX4 and the transcription factors described above, ChIP-seq data for histone modifications were queried to gain insight into potential promoter or enhancer regulatory functions for the identified transcription factor binding sites. For this, histone H3K4 tri-methylation (found near promoters), H3K4 mono-methylation (found near regulatory elements), and H3K27 acetylation (found near active regulatory elements) were included. These histone modifications were queried in ChIP-seq datasets using all 9 available cell line tracks.

Pathway analysis was performed using Ingenuity Pathway Analysis software version 52912811. Candidate miRNAs from these studies were uploaded along with transcription factors whose dysregulation is associated with FSHD. Defined network connections were identified using the Pathway Builder application. Molecules confirmed to have established relationships were used to visualize a novel network built from these FSHD expression data.

#### *2.5. Expression of Individual miRNAs in a Validation Sample Set*

Circulating miRNAs that were significantly up-regulated in individuals affected by FSHD1 were examined in a separate set of non-overlapping CINRG patients used as a validation group. For this group, FSHD patients had a confirmed diagnosis of FSHD1 (*n* = 12; 9 females, 3 males) and were compared to healthy volunteer control samples (*n* = 7; 4 females, 3 males). RNA was isolated from 150 μL of plasma using Trizol LS liquid extraction. Total RNA was converted to cDNA using a High Capacity Reverse Transcription Kit with multiplexed RT primers, pre-amplified using PreAmp MasterMix with multiplexed TM primers, and quantified with individual TaqMan assays on an ABI QuantStudio 7 real time PCR machine (Applied Biosystems; Foster City, CA, USA). Assay IDs used are: miR-32-002109, miR-103-000439, miR-505-002089, miR-146b-001097, miR-29b-000413, miR-34a-000426, miR-141-000463, miR-98-000577, miR-576-3p-002351, miR-9-000583, and miR-142-3p-000464. Expression levels of all miRNAs were normalized to the geometric mean of multiple control genes (miR-150 and miR-342-3p) determined previously to be stable circulating miRNA controls [35,51]. Expression was analyzed in FSHD1 versus healthy control patients via *t*-test analysis, including assessment of directionality. A *p*-value of ≤ 0.05 was considered significant. Data are presented as mean ± SEM unless otherwise noted.

#### *2.6. Proteomics Profiling*

Plasma samples were first processed using Pierce™ Top 12 Abundant Protein Depletion Spin Columns (Thermo Scientific) before mass spectrometry analyses using the Q Exactive HF mass spectrometer. Briefly, the 12 most abundant proteins from 5 μL of plasma sample were affinity depleted by incubating with Top 12 protein depletion resin. Following this, the unbound fraction was collected according to the manufacturer's protocol. Proteins were precipitated with pre-cooled acetone (1:5 vol) for 30 min at −20 ◦C and centrifuged at 4 ◦C for 15 min at max speed in a micro-centrifuge. The liquid was decanted and the pellet was air dried briefly and resuspended with 8 M Urea, followed by reduction and alkylation with 5 mM DDT and 15 mM idodoacetamide for 30 min at room temperature. Samples were diluted with 100 mM ammonia bicarbonate to final urea concentration of less than 2 M. Afterwards, the samples were digested with 1 μg of trypsin (Promega) at 37 ◦C overnight. Trypsin was inactivated by 0.1% TFA and samples were desalted by capturing the peptides onto C18 100 μL bed tips (Pierce®C18 tips, Thermo Scientific) following the manufacture's protocol. The bound peptides were eluted with 60% acetonitrile, 0.1% TFA, then dried using a SpeedVac, and resuspended in 20 μL buffer containing 2% acetonitrile with 0.1% acetic acid.

The peptide mixtures from each fraction were sequentially analyzed by liquid chromatography tandem mass spectrometry (LC-MS/MS) using Thermo Ultimate 3000 RSLCnano-Q Exactive mass spectrometry platform nano-LC system (Easy nLC1000) connected to Q Exactive HF mass spectrometer (Thermo Scientific). This platform is configured with nano-electrospray ion source (Easy-Spray, Thermo Scientific), Acclaim PepMap 100 C18 nanoViper trap column (3 μm particle size, 75 μm ID × 20 mm length), EASY-Spray C18 analytical column (2 μm particle size, 75 μm ID × 500 mm length). The data from each sample was collected in triplicate at 2 μL per injection, following which the peptides were eluted at a flow rate of 300 nL/min using linear gradients of 7–25% Acetonitrile (in aqueous phase and 0.1% Formic Acid) for 80 min, followed by to 45% for 25 min, and static flow at 90% for 15 min. The mass spectrometry data was collected in data-dependent manner switching between one full scan MS mode (*m*/*z* 380–1600, resolution 70,000, AGC 3e6) and 10 MS/MS mode (resolution 17,500); where MS/MS analysis of the top 10 target ions were performed once and dynamically excluded from the list for 30 s.

The MS raw data sets were searched against UniProt human database that included common contaminants using MaxQuant software (version 1.5.5.1) [52]. We used default parameters for the searches, first search peptide tolerance 20 ppm, main search peptide tolerance 4.5 ppm, maximum two missed cleavage; and the peptide and resulting protein assignments were allowed at 0.01 FDR (thus 99% confidence level). Protein levels were quantified in 25 FSHD1 patients and 17 healthy controls and reported for each protein as the number of unique peptides detected and the intensity measured. Proteins with altered abundance with greater than 2-fold were selected for further inquiry.

Several pre-processing steps were performed on the raw data values before statistical analysis. Each sample had either 2 or 3 replicates which were averaged to yield a single quantification for each subject for each protein. When a value of zero occurs, it can indicate either a true zero or an assay that did not detect that protein. To accurately reflect protein levels, we incorporated zeroes into our analysis in the following way. If one replicate yielded a zero value, that zero was left as is and treated as a true zero. If two replicates yielded a zero, all values for that protein/sample were set to missing as we cannot distinguish true zeroes from artificial ones. We then applied a normalization factor to the average values to account for differences in the amount assayed per sample. We summed the protein counts for all proteins for each sample and used the maximum value to normalize all other samples. This allowed us to ensure that the amount of proteins assayed were proportional for all samples.

All values were log-transformed for analysis. We assessed the relationship between protein levels and disease severity in the FSHD1 patients using a linear regression model where protein level was the dependent variable, severity was the independent variable, and age and gender were covariates. Regression models were performed only for proteins found in 5 or more samples. Model estimates were reported for each protein and included the coefficient and *p*-value for all terms in the model (severity, age and gender) along with an indication of the direction of each effect. This same method was used to assess the relationship between protein level and the number of D4Z4 repeats. We assessed the difference in protein expression between FSHD1 patients and healthy controls using a linear regression model where protein level was the dependent variable, a categorical indicator of disease was the independent variable, and age and gender were covariates. Again, regression models were performed only for proteins found in 5 or more samples. Model estimates were reported for each protein and included the coefficient and p-value for all terms in the model (disease status, age and gender), an indication of the direction of each effect, and age and gender adjusted means for each disease group. As this part of the analysis was discovery in nature, we did not adjust resulting *p*-values for multiple testing. Our intention was to find those proteins showing some evidence of an effect and to move those proteins forward for an additional evidence-based validation experiment. The significance level for all analyses was set at 0.05.

#### *2.7. Enzyme-Linked Immunosorbent Assay (ELISA)*

Five proteins were chosen for further validation in a separate set of patients via protein-specific ELISA assays. Human specific protein ELISA kits for human insulin-like growth factor-1 (IGF1) (R&D Systems, Minneapolis, MN, USA), profilin 1 (PFN1) (LSBio, Seattle, WA, USA), S100 Calcium-Binding Protein A8 (S100-A8) (Biotechne, Minneapolis, MN, USA), Proteoglycan 4 (PRG4) (AVIVA Systems Biology, San Diego, CA, USA), Human Tropomyosin alpha-4 chain (TPM4) (MyBioSource, San Diego, CA, USA) were performed to determine protein level in FSHD1 and unaffected controls. Plasma (20 μL) from FSHD1 patients (*n* = 19) and healthy volunteer (*n* = 13) controls (age and gender matched) were tested in duplicate following the manufacturer's recommended protocols. ELISA values were assessed for normality and a log-transformation applied where appropriate. We assessed the relationship between protein level and severity using, as described above, a linear regression model where protein level was the dependent variable, severity was the independent variable, and age and gender were covariates. We assessed the difference in protein expression between FSHD1 and healthy controls using a linear regression model where protein level was the dependent variable, a categorical indicator of disease was the independent variable, and age and gender were covariates. All analyses were performed at the 0.05 significance level.

#### **3. Results**

#### *3.1. Discovery of Novel Candidate miRNA Biomarkers Associated with FSHD*

Sixteen FSHD1 patients with pediatric-onset, matched for sex and age, were selected into two groups of a discovery sample set for circulating biomarker studies: one mild FSHD1 group (*n* = 8), and one severe FSHD1 group (*n* = 8), as determined by an FSHD disease severity score. These two groups were each compared to a group of healthy control volunteers (*n* = 8). Demographics are displayed in Table 1. Patients with severe FSHD1 showed a significantly higher FSHD severity score (12.25 ± 2.76; *p* ≤ 0.00001) than patients with mild FSHD1 (4.88 ± 1.46), with any value of nine or higher being classified as severe FSHD.


**Table 1.** Clinical characteristics of study group patients.

\*\* *p* ≤ 0.00001, *t*-test of mild FSHD versus severe FSHD severity score.

Ten miRNAs showed a significant change in expression level in mild FSHD1 plasma versus healthy controls, and 12 miRNAs showed a significant change in expression level in severe FSHD samples versus controls (Table 2). Of these, three miRNAs showed a significant increase in both mild and severe FSHD1 in comparison to healthy controls: miR-32, miR-505, and miR-29b. Each of these three miRNAs showed an approximately two-fold higher change in expression in severe FSHD1 patients than in mild FSHD1 patients versus healthy controls. Of the 19 unique miRNAs identified, several have been previously found to play a role in muscle disease pathways. miR-29b, which is associated with TGFβ-signaling and fibrosis, was up-regulated in both mild and severe FSHD1 patients. Both miR-146b and miR-142-3p, which are known to be up-regulated in inflammatory disease states, were up-regulated in mild FSHD1 patients and have previously been shown to be up-regulated in dystrophinopathy (Becker and Duchenne muscular dystrophy) patients and/or animal models [30,36]. miR-486 has previously been defined as a muscle-enriched microRNA or "myomiR" [53], and was found here to be down-regulated in mild FSHD1 patients (*p* < 0.005).




**Table 2.** *Cont.*

*Italics* = dysregulated in both mild and severe FSHD; ACAD = acute coronary artery disease, BMD = Becker muscular dystrophy, COPD = chronic obstructive pulmonary disease, DMD = Duchenne muscular dystrophy, IBD = inflammatory bowel disease, LMNA = Lamin A/C, TGFβ = Transforming Growth Factor β, VSMC = vascular smooth muscle cell. \* *p* < 0.005.

#### *3.2. Bioinformatic Analysis of miRNA Regulation and Pathways*

To examine their regulation by transcription factors which are dysregulated by the FSHD disease process, we next performed bioinformatic analyses of ChIP-seq data for DNA binding by transcription factors in proximity to each candidate miRNA's genomic locus. To gain insight into direct consequences of *DUX4*-up-regulating mutations that cause FSHD, we analyzed ChIP-seq data for DUX4. To do this, we analyzed DUX4 binding via a user-supplied DUX4 ChIP-seq track published by Geng et al. [9]. Genes for 16 of the candidate miRNAs had at least one binding site within distances capable of providing gene enhancer functions (Figure 1). Examination of the miR-100 home gene (*MIR100HG*) locus was particularly interesting. In total, we found 18 DUX4 binding sites in the area surrounding *MIR100HG*, and many of these clearly overlapped with histone modifications associated with active promoters (H3K4 tri-methylation) and regulatory elements (H3K27Ac). These data are consistent with regulation of miR-100 expression by DUX4 (Figure 1b).

To gain insight into additional pathways that may drive expression of candidate miRNAs and contribute to FSHD molecular pathophysiology, we performed bioinformatic analyses of ChIP-seq data for transcription factors that are dysregulated as a result of DUX4 mutations. For this, we obtained a list of transcription factors which are expressed at significantly different levels in human skeletal muscle cells as a result of DUX4 overexpression [9]. We then queried publicly available ChIP-seq datasets to identify which of these transcription factors had ChIP-seq datasets available through the ENCODE public research consortium [46,47]. Of the transcription factors in this dataset, 34 had ChIP-seq datasets available in the Factorbook repository and 47 had ChIP-seq datasets available in the ENCODE 3 repository [48–50]. Genomic binding by each of these transcription factors was surveyed for each of these transcription factors for all candidate miRNAs (Table S1). Transcription factors that increased in response to overexpression of toxic, full-length DUX4 but did not increase in response to a non-toxic, truncated isoform of DUX4 were considered to be of particular interest (Figure 2a). Of these factors, four showed a particularly high number of binding sites within regulatory distance of the candidate miRNAs: early growth response protein 1 (EGR1), FOS, MYC, and yin yang 1 (YY1). As an example of these findings, miR-576 was up-regulated in FSHD patients, has five DUX4 binding sites neighboring its home gene (SEC24B), and has a high number of binding sites for the secondary transcription factors described here (Figure 2b). EGR1, FOS, MYC and YY1 all showed a large number

of binding sites around miR-576, and these frequently overlapped with histone modifications which mark active promoter and enhancer regions, consistent with these four transcription factors driving gene expression signatures in FSHD.

**Figure 1.** DUX4 binding sites at loci surrounding miRNAs dysregulated in FSHD patients. The 19 miRNAs dysregulated in FSHD1 patient plasma samples were queried for potential regulation by the DUX4 transcription factor, which aberrantly expressed in FSHD, using a DUX4 ChIP-seq dataset [9]. (**a**) Overview of all DUX4 binding sites within regions capable of acting as regulatory elements (100 kb) of the 19 miRNAs and their home genes. (**b**) Schematic of DUX4 binding sites within the miR-100 locus and its surrounding home gene (*MIR100HG*) variants. Note, miR-100 is transcribed from right to left on this image. Corresponding epigenetic modification maps display the location of histone modifications associated with active promoters (H3K4me3) and poised/active enhancers (H3K4me1 and H3K27Ac, respectively).

**Figure 2.** Candidate miRNA loci are consistent with regulation via transcription factors dysregulated in FSHD. (**a**) Table listing a subset of transcription factors which are each increased in human skeletal muscle cells in response to DUX4 overexpression [9], along with the number (#) of binding sites they show within potential regulatory distance (100 kb) of the 19 candidate miRNAs. (**b**) The miR-576 locus shows binding consistent with regulation by FOS, EGR1, MYC, YY1, and DUX4. Corresponding epigenetic modification maps display the location of histone modifications associated with active promoters (H3K4me3) and poised/active enhancers (H3K4me1 and H3K27Ac) in the vicinity of the miR-576 locus and its surrounding home gene, *SEC24 homolog B* (*SEC24B*). (DUX4 binding sites identified using ChIP-seq data uploaded from Geng et al. [9]; binding sites for additional transcription factors identified using UCSC Genome Browser and respective ChIP-seq datasets accessed via the ENCODE3 regulation track [46–50]).

Additionally, we used Ingenuity Pathway Analysis software to perform a bioinformatic analysis on the candidate miRNAs identified in this study, together with transcription factors previously published to be dysregulated in FSHD [9], to see if there are defined signaling pathways or interactions shared by these factors. Interestingly, this analysis showed that there are previously established connections between many of the miRNAs and transcription factors examined, with 15 of the miRNAs and 18 of the transcription factors found to make up a network with previously defined interactions (Figure 3). For example, increased levels of miR-34a are known to decease cAMP response element-binding protein (CREB) to drive neuronal dysfunction in HIV-induced neurocognitive disorders, and to increase AMP-dependent transcription factor 3 (ATF3) levels in colon cancer [77,78]. MYC binds to ATF3 as well as to lysine-specific demethylase 5B (KDM5B) and YY1, all four of which are elevated in FSHD [9,79–81]; in addition, MYC is known to activate transcription of both enhancer of zeste homolog 2 (EZH2) and miR-9 [82,83], both of which are also increased in FSHD. Together, these bioinformatics data show our candidate miRNA markers are consistent with a change in transcriptional programming that results from FSHD-causing DUX4 overexpression mutations.

**Figure 3.** Pathway analysis of miRNAs and transcription factors dysregulated by FSHD mutations. Ingenuity Pathway Analysis software was used to identify established connections between candidate miRNAs from this study with transcription factors known to be dysregulated by FSHD-causing overexpression of DUX4 [9]. Red-shaded miRNAs and transcription factors were observed to increase, while those shaded blue were observed to decrease. Solid arrows denote direct relationships, while dashed arrows denote indirect relationships.

#### *3.3. Confirmation of miRNA Increases in FSHD1 Patients*

Next, we assayed expression of candidate miRNA biomarkers in samples from a separate and non-overlapping group of patients. Upon clinical examination, all patients in this validation group were determined to have FSHD1. We selected 14 miRNAs that significantly increased in the discovery experiments for follow-up study in the validation group. We found three of these miRNAs (miR-9, miR-32 and miR-329) were not expressed at consistently high enough levels for detection within plasma from the validation set of FSHD1 patients, leaving 11 miRNAs for validation. Here, these 11 individual candidate miRNAs were quantified in FSHD1 (*n* = 12; 9 females, 3 males) versus healthy volunteer control samples (*n* = 7; 4 females, 3 males).

Upon quantification, we found 8 of these 11 candidate miRNAs also showed a clear increase in samples from the FSHD validation group in comparison to healthy controls (Figure 4). miR-100, miR-103, miR-29b, miR-34a, miR-454, miR-505 and miR-576 were all expressed at significantly higher levels (*p* ≤ 0.05) in FSHD1 serum. miR-100, miR-29b, miR-34a, miR-505, and miR-576 were the most highly up-regulated in FSHD1, showing up-regulation from approximately 4- to 20-fold higher than healthy controls. miR-146b was also expressed at an approximately 2-fold higher level in this set of FSHD1 patients; however it did not reach significance (*p* = 0.06). Of the remaining three miRNA candidates, miR-98 showed no apparent change, while miR-141 and miR-142-3p showed an approximately 50% increase that did not reach significance. As most candidate miRNAs showed consistent behavior in this separate validation set of FSHD samples, this panel of miRNAs merits further investigation as biomarkers moving forward.

**Figure 4.** Expression of candidate miRNAs in a validation group of FSHD patients. Candidate miRNAs that increased in the FSHD discovery experiment were assayed via individual qRT-PCR assay in a separate validation group of FSHD1 patient plasma samples. Expression levels of each miRNA are expressed as fold change versus healthy control volunteers. (values are mean ± SEM, \* *p* ≤ 0.05, \*\* *p* ≤ 0.01, one-tailed *t*-test comparing FSHD1 to control in direction of Discovery experiment; one outlier removed from miR-34a and miR-576 after significant Grubb's outlier test; *n* = 7 healthy control volunteers, 12 FSHD1).

#### *3.4. Proteomics Profiling*

To identify protein candidate biomarkers, we performed LC-MS/MS-based proteomic profiling of samples from a discovery group of FSHD patients (Table 3). For this, plasma from FSHD1 patients (*n* = 25) was compared to healthy volunteer controls (*n* = 17), with a roughly even mix of males and females, and an average age of early- to mid-twenties for each group. All FSHD patients were confirmed to have FSHD1 resulting from D4Z4 contraction mutations that alter epigenetic regulation of DUX4.


**Table 3.** Clinical characteristics of patients in proteomics discovery group.

Based on signal intensity, we identified 32 proteins that were significantly different between FSHD1 and healthy control samples (Table S2). To further filter the protein list, we used unique peptide count data to identify proteins that had significantly different counts between FSHD1 and control samples. This narrowed the candidates down to 13 proteins (Table 4). Within these protein markers, fibulin-1 (FBLN1) and insulin-like growth factor 1 (IGF1) showed potential effects of sex and age, while keratin 16 (KRT16) displayed a potential age effect and profilin-1 (PFN1) showed a potential sex effect (Table S3). Among the 13 total protein biomarker candidates, 11 proteins were higher in FSHD1 samples versus healthy controls, while two proteins were lower in the FSHD1 samples versus healthy controls.

**Table 4.** Thirteen circulating proteins identified as dysregulated in FSHD plasma via LC-MS/MS.


CFL1 = Cofilin 1, EFEMP1 = EGF-containing fibulin-like extracellular matrix protein 1, F13A1 = Coagulation factor XIII A chain, FBLN1 = fibulin-1, IBD = inflammatory bowel disease, IGFI = insulin-like growth factor 1, KRT16 = Keratin 16, PFN1 = Profilin-1, PRG4 = Proteoglycan 4 or lubricin, PROC = Protein C, S100A8 = S100 calcium-binding protein A8, SPP2 = Secreted phosphoprotein 24, TLR4 = Toll-like receptor 4, TMSB4X = Thymosin beta-4, TPM4 = Tropomyosin alpha-4 chain.

We selected five candidate protein markers for subsequent quantification via protein-specific ELISA analysis of a non-overlapping validation group of FSHD1 samples. These included insulin-like growth factor 1 (IGF1), proteoglycan 4 (PRG4), profilin 1 (PFN1), tropomyosin 4 (TPM4), and S100 calcium-binding protein A8 (S100A8). Of these candidate proteins, S100A8 showed a significant increase in FSHD1 plasma of approximately 4.5-fold over healthy controls in the validation group (Figure 5a), consistent with its behavior in the discovery experiment. To determine if elevated S100A8 signaling was consistent with the overall proteomic and miRNA profiling results, we performed bioinformatic pathway analyses focused on the S100A8 pathway along with the full list of candidate protein (Figure 5b) and miRNA (Figure 5c) markers. Nine proteins and 13 miRNAs were shown to have previously established connections to the toll-like receptor 4 (TLR4) signaling pathway, which is activated by S100A8 and drives increased inflammatory (NF-κB and AP-1) gene expression. As miRNAs can reflect a direct readout of transcription factor activity, we also surveyed ChIP-seq data to analyze DNA regions encoding miRNAs elevated in FSHD for binding by the NF-κB and AP-1 transcription factors activated by S100A8 (Figure 5d). All miRNAs except for one (miR-329) showed binding by NF-κB and/or AP-1 subunits at DNA regions capable of acting as regulatory promoter or enhancer elements. As S100A8 is a well-established biomarker of inflammatory disease processes (reviewed in [86]) and these can be up-regulated in the muscular dystrophies, this protein merits further investigation as a biomarker for FSHD moving forward.

**Figure 5.** Validation and pathway analysis of elevated S100A8 protein in FSHD. (**a**) ELISA of S100A8 protein in plasma from a separate validation set of FSHD1 patients. (**b**) Bioinformatic pathway analysis was used to identify known connections between candidate protein markers with S100A8 pathway proteins involved in TLR4 signaling. (**c**) Bioinformatic pathway analysis was used to identify established connections between candidate miRNAs with S100A8 pathway proteins involved in TLR4 signaling. (**d**) Bioinformatic analysis of ChIP-seq defined binding sites for the key S100A8 pathway transcription factors AP-1 (FOS and JUN) and NF-κB (RELA), at potential regulatory regions of the candidate miRNAs that were found to increase in FSHD plasma. Binding sites represent the combined number of potential promoter (within 2 kb of promoter) and enhancer (within 10 kb) regulatory regions with ChIP-seq-confirmed transcription factor binding for each miRNA home gene. (\*\* *p* ≤ 0.01; *n* = 13 healthy control volunteers, 19 FSHD1; panels (**b**,**c**) produced using Ingenuity Pathway Analysis software, red = increased, blue = decreased; data for panel (**d**) produced using the Factorbook ChIP-seq data repository from ENCODE and the UCSC genome browser).

#### **4. Discussion**

There is currently no effective treatment available for FSHD. However, research advances in FSHD are now beginning to yield promising and novel therapeutic strategies that will require well-designed clinical trials to evaluate effectiveness. Potential therapeutic strategies including antisense oligonucleotides (AON) and small molecules have been reported or are being actively pursued [12,97–100]. Changes in biomarkers following a treatment can be a powerful tool for evaluating the efficacy and safety of the treatment. Previous studies seeking to identify circulating miRNA biomarkers in muscular dystrophy have focused exclusively on assaying myomiRs, which are a defined group of miRNAs with muscle-specific or muscle-enhanced expression [101,102]. Previously, a study by Statland et al. identified 7 potential protein biomarkers in 22 FSHD serum samples, using a commercial multiplex assay [103]. A multi-site study using aptamer-based SomaScan proteomics to assay two FSHD populations identified a total of 115 proteins that were dysregulated, four of which behaved consistently between the two independent cohorts (creatine kinase MM, creatine kinase MB, carbonic anhydrase III, and troponin I type 2) [104]. In this study, we used -omics approaches to identify additional circulating miRNA and protein biomarker candidates using samples collected from individuals with early-onset FSHD.

There is an intriguing potential for developing miRNAs as biomarkers in diseases affecting muscle, as they are stable in biofluids, objective, minimally invasive, and well-conserved between human patients and preclinical animal models [23,24]. Recently, the utility of serum miRNAs to detect muscle involvement in complex diseases with highly variable symptoms has been demonstrated, as in patients with MADD [25]. Muscle-specific miRNAs are also elevated in Duchenne and Becker muscular dystrophy, along with a set of inflammatory miRNAs reflecting the chronic inflammatory pathology of these diseases [29,30,105]. Here we identify eight circulating miRNAs that are associated with FSHD in patient plasma samples. The prevalence of DNA binding by DUX4 and FSHD-associated transcription factors, within regions capable of regulating the candidate miRNAs, provides a molecular rationale for their up-regulation in FSHD. Several of the markers have also been previously shown to play a role in muscle diseases and associated pathological pathways. These candidate biomarkers hold potential as monitoring biomarkers in early-onset FSHD.

Several candidate miRNAs we identified have previously been proposed as circulating biomarkers and have shown similar behavior in other diseases. Plasma miR-454 has been identified as a biomarker of myotonic dystrophy [72,73]. Serum miR-146b is a pharmacodynamic biomarker in inflammatory bowel disease (IBD) [34,35]. Intriguingly, miR-146b is also known to down-regulate dystrophin in multiple muscle diseases, is increased in dystrophinopathies and in myositis, and is also drug-responsive in the *mdx* mouse model of DMD [30,31]. Urinary miR-141 provides a promising diagnostic biomarker for the identification of both prostate and bladder cancers [69]; it will be interesting to determine if this or other candidate miRNAs are also dysregulated in urine from dystrophic patients, as this sampling method could provide a completely non-invasive biomarker.

Increases in circulating S100A8, a subunit of calprotectin, are consistent with an inflammatory signature playing a role in FSHD. The inflammatory calprotectin protein consists of a heterodimer (S100A8/S100A9) which binds to toll-like receptor 4 (TLR4) to activate pro-inflammatory gene expression pathways through the NF-κB and AP-1 transcription factors. Consistent with such an inflammatory gene signature in FSHD, bioinformatic analyses here show five of the candidate miRNAs have established connections with TLR4 signaling, are increased in FSHD patients, and have gene promoters that are bound by AP-1 and/or NF-κB. Outside of FSHD, calprotectin is already a well-established biomarker across rheumatic diseases. Fecal calprotectin is a widely used diagnostic, monitoring and pharmacodynamic biomarker for IBD, and recent studies indicate serum calprotectin levels are also well-correlated with IBD disease state [87,106]. Serum calprotectin is used as a monitoring and pharmacodynamic biomarker for rheumatoid arthritis, and intriguingly S100A8/S100A9 may have further utility in arthritis as a molecular imaging marker of inflammatory activity [84,107,108]. Of particular relevance to the present study, calprotectin in both muscle and serum is a biomarker

for disease activity in juvenile dermatomyositis [109]. Moving forward, it will be interesting to see if S100A8 or calprotectin can show further utility as completely non-invasive or local biomarker for FSHD and other muscle diseases such as myositis.

Several of the molecular markers we identified here as elevated in FSHD may provide a new therapeutic target. In various states of muscle atrophy miR-29b is also up-regulated, while preventing its expression shows efficacy in mouse models of muscle atrophy [64,65]. In myositis and Becker muscular dystrophy, the inflammatory marker miR-146b is known to down-regulate dystrophin expression, whereas the reduction of miR-146b via anti-inflammatory drugs or via miRNA-targeting oligos is proposed as a method to increase dystrophin levels to help improve muscle health [30,31]. In various rheumatological disease states, the inhibition of S100A8 or calprotectin via small molecule inhibitors or antibodies is a very attractive therapeutic strategy; early studies of such inhibitors are already showing therapeutic efficacy in both human trials and/or in mouse models, including in studies for arthritis, asthma, IBD, and multiple sclerosis (reviewed in [86]). Similarly, decreases in PROC seen here in FSHD are also seen in several rheumatological disorders, where treatment with PROC activators are already being pursued as a therapeutic option (reviewed in [94]).

Bioinformatic analyses of the -omics results support muscle and inflammatory gene expression pathways as being dysregulated in FSHD. As expected, several muscle pathology-associated miRNAs are dysregulated in FSHD patients: miR-486 is a defined myomiR, miR-29b up-regulation promotes muscle atrophy, miR-146b is dysregulated in dystrophinopathies and myositis, miR-329 counteracts muscle hypertrophy, and three others are known to be dysregulated in myotonic dystrophy, lamin A (LMNA) dystrophy, and/or FSHD (miR-34a, miR-140-3p, miR-100, and miR-454). Consistent with these findings, several of the proteins that were dysregulated are known to function in muscle contraction, actin filament organization and/or muscle regeneration (TOM4, PFN1, CFL1, TMSB4X, and IGF1).

S100A8 and its associated inflammatory signaling pathway (TLR4, NF-κB and AP-1) appear to be a substantial hub for dysregulated expression of the candidate markers we identified. Nine of the candidate miRNAs have previously established connections to this TLR4-centered pathway. ChIP-seq analysis of the miRNAs up-regulated in FSHD shows all but one have promoters bound by NF-κB or AP-1, which are activated by S100A8-induced TLR4. In the proteomics data, several of the proteins that increased are pro-inflammatory (S100A8, KRT16 and SPP2) while in contrast the two proteins that decreased have anti-inflammatory (PROC and PRG4) roles. Consistent with our FSHD findings, KRT16 and S100A8 are also up-regulated together in inflammatory skin disorders; additionally, the pattern of increased S100A8 with decreased PROC is seen here in FSHD as well as in IBD and several other chronic inflammatory disorders [93,94]. Pathway analysis further establishes a link between the protein markers, as nine out of 14 have established connections to the S100A8 and TLR4 signaling pathway. Together these data confirm that circulating FSHD biomarkers reflect muscle pathogenesis, and suggest inflammatory S100A8/TLR4 signaling plays a role in pediatric-onset FSHD as well.

#### **5. Conclusions**

FSHD is chronic genetic muscle disease with a variable prognosis. There is no cure, and no pharmaceuticals for FSHD have shown efficacy in altering the disease course. Development of objective biomarkers will facilitate the clinical and preclinical development of novel therapies, as well as our ability to monitor disease activity. We identified eight circulating miRNAs (miR-100, miR-103, miR-146b, miR-29b, miR-34a, miR-454, miR-505, and miR-576) which may be developed as biomarkers for FSHD. Additionally, we identified the S100A8 subunit of calprotectin as a primary protein marker of interest for FSHD, consistent with its utility in numerous rheumatic diseases. These molecular markers warrant further investigation in additional cohorts, preclinical drug testing, and prospective clinical trials.

*J. Pers. Med.* **2020**, *10*, 236

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2075-4426/10/4/236/s1, Table S1: ChIP-seq of FSHD-disrupted transcription factors at candidate miRNA loci, Table S2: Total proteomic changes in FSHD plasma as detected via LC-MS/MS, Table S3: Proteins affected by sex or age.

**Author Contributions:** Conceptualization, C.R.H., A.A.F., Y.-W.C.; methodology, C.R.H., A.Z., N.Y.N., C.B.T., A.P., H.G.-D., A.A.F., M.G., M.M.R., P.R.C., M.T. (Mathula Thangarajh), R.W., E.C.S., A.M.C., C.M.M., P.K., M.T. (Mar Tulinius), A.H., J.K.M., A.A.F., Y.-W.C.; data analysis, C.R.H., A.Z., N.Y.N., H.G.-D., A.A.F., Y.-W.C.; experimental procedure carried out by, C.R.H., A.Z., N.Y.N., C.B.T., S.N.P., A.P.; writing-original draft preparation, C.R.H., A.Z., N.Y.N., C.B.T., A.P., H.G.-D., A.A.F., Y.-W.C.; writing-review and editing, C.R.H., A.Z., N.Y.N., C.B.T., H.G.-D., S.N.P., A.P., M.G., M.M.R., P.R.C., M.T. (Mathula Thangarajh), R.W., E.C.S., A.M.C., C.M.M., P.K., M.T. (Mar Tulinius), A.H., J.K.M., A.A.F., Y.-W.C.; project administration, C.R.H., M.G., M.M.R., P.R.C., M.T. (Mathula Thangarajh), R.W., E.C.S., A.M.C., C.M.M., P.K., M.T. (Mar Tulinius), A.H., J.K.M., A.A.F., Y.-W.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study is supported by funding from the Friends of FSH Research, the FSHD Society, Muscular Dystrophy Canada, and the FSHD Global Research Foundation. C.R.H. and this work is supported by the National Institutes of Health (NIH) grants R01HL153054, P50HD090254, and R00HL130035, as well as the Foundation to Eradicate Duchenne and the A. James & Alice B. Clark Foundation. Y.W.C. is partially supported by NIH/NINDS 1R03NS116444, Muscular Dystrophy Association and NIH/NCATS UL1TR001876.

**Acknowledgments:** The authors would like to thank participants and families for their support, and the FSH Society for assistance with patient travel and care. The authors would like to thank the Cooperative International Neuromuscular Research Group (CINRG) for patient recruitment and care, as well as the CINRG Early-onset FSHD investigators who provided patient samples for this study. The authors would also like to thank the ENCODE Consortium, the ENCODE ChIP-seq production laboratories, and the ENCODE Data Coordination Center for generating and processing ChIP-seq datasets used here.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Clinical and Laboratory Associations with Methotrexate Metabolism Gene Polymorphisms in Rheumatoid Arthritis**

**Leon G. D'Cruz 1,2, Kevin G. McEleney 1, Kyle B. C. Tan 1, Priyank Shukla 1, Philip V. Gardiner 3, Patricia Connolly 4, Caroline Conway 5, Diego Cobice <sup>5</sup> and David S. Gibson 1,\***


Received: 20 August 2020; Accepted: 24 September 2020; Published: 26 September 2020

**Abstract:** Rheumatoid arthritis (RA) is a chronic systemic autoimmune disease that causes loss of joint function and significantly reduces quality of life. Plasma metabolite concentrations of disease-modifying anti-rheumatic drugs (DMARDs) can influence treatment efficacy and toxicity. This study explored the relationship between DMARD-metabolising gene variants and plasma metabolite levels in RA patients. DMARD metabolite concentrations were determined by tandem mass-spectrometry in plasma samples from 100 RA patients with actively flaring disease collected at two intervals. Taqman probes were used to discriminate single-nucleotide polymorphism (SNP) genotypes in cohort genomic DNA: rs246240 (*ABCC1*), rs1476413 (*MTHFR*), rs2231142 (*ABCG2*), rs3740065 (*ABCC2*), rs4149081 (*SLCO1B1*), rs4846051 (*MTHFR*), rs10280623 (*ABCB1*), rs16853826 (*ATIC*), rs17421511 (*MTHFR*) and rs717620 (*ABCC2*). Mean plasma concentrations of methotrexate (MTX) and MTX-7-OH metabolites were higher (*p* < 0.05) at baseline in rs4149081 GA genotype patients. Patients with rs1476413 SNP TT or CT alleles have significantly higher (*p* < 0.001) plasma poly-glutamate metabolites at both study time points and correspondingly elevated disease activity scores. Patients with the rs17421511 SNP AA allele reported significantly lower pain scores (*p* < 0.05) at both study intervals. Genotyping strategies could help prioritise treatments to RA patients most likely to gain clinical benefit whilst minimizing toxicity.

**Keywords:** rheumatoid arthritis; SNP; DMARD; methotrexate; pharmacogenomics

#### **1. Introduction**

Rheumatoid arthritis (RA) is the most common chronic autoimmune inflammatory arthritis, affecting approximately 0.3–1% of the world's population [1,2]. The disease primarily affects the articular joints, causing swelling, stiffness, joint destruction [3], loss of function in joints [4], disability and a significantly lower quality of life. To prevent irreversible joint damage resulting in substantial

disability, it is important to introduce disease-modifying anti-rheumatic drugs (DMARDs) early after onset and failure of non-steroidal anti-inflammatory treatment.

Conventional synthetic disease-modifying anti-rheumatic drugs (csDMARDs) such as methotrexate (MTX), hydroxychloroquine (HCQ), cyclosporin, sulfasalazine (SSZ) and leflunomide are commonly used mainstays of the disease; however, it is widely known that a significant proportion of patients with RA often show poor or inadequate therapeutic response to csDMARDs [5]. The anti-folate MTX is the cheapest drug in treatment of RA and is often the first-line treatment [6]; however, only 55% of patients remain on this drug for more than 2 years due to a build-up of non-response or the accumulation of various adverse side effects [6,7]. MTX is subject to significant metabolic activity in the body; the polyglutamated derivatives of MTX are selectively retained in cells, therefore lengthening the activity of the drug which complicates treatment management, since patients would continue taking their daily drug dosage oblivious to the fact that their circulating drug levels are still high, potentially contributing to undesirable cytotoxic effects [8,9]. MTX is converted in hepatic parenchymal cells resulting in the 2- through 4-glutamate residues derivatives or the drug is catabolised to the 7-hydroxy-methotrexate (MTX-7-OH) form [10]. More than 10% of a dose of methotrexate is oxidised to MTX-7-OH, irrespective of the route of administration [11]. The MTX-7-OH metabolite is extensively (91 to 93%) bound to plasma proteins, in contrast to the parent drug (only 35 to 50% bound) and contributes to inactivity of the drug or poor response to treatment [11].

When non-response has been confirmed, NICE clinical guidelines recommend switching to the more costly biological disease-modifying anti-rheumatic drugs (bDMARDs) [5,12,13]. Various studies indicate that treat-to-target strategies which aim to reduce disease activity shortly after diagnosis result in better long term outcomes and can minimise permanent joint damage, thus there is a genuine need for earlier identification of patients who do not respond well to csDMARDs treatments [6].

It is estimated that 15–30% of variation in drug responses are attributable to genetic or single-nucleotide polymorphisms or SNPs [14]. Not all SNPs are functional; some are in non-coding areas (introns) and there is a variety of ways that a SNP can affect or inhibit downstream transcription factor, gene or protein function [15]. The promise of pharmacogenomics is that identification of SNPs and associated risk alleles could identify patients who may be susceptible to accumulating cytotoxic levels of a drug during therapy (as in the polyglutamate derivatives of MTX) or when certain metabolite levels accumulate rendering a drug as inactive (as in MTX-7-OH).

In this study, we sought to determine the metabolite levels in RA patients taking DMARDs. We then carried out genotyping of 10 SNPs known to influence the metabolic pathways of DMARDs in arthritis. Our aim was to determine if genetic variations or polymorphisms associate with metabolite levels. This could help design studies to improve clinical management, which risk stratify patients at greater predisposition of forming ineffective or potentially harmful metabolite levels, by adequately planning ahead the appropriate drug and dosage.

#### **2. Materials and Methods**

#### *2.1. Participant Recruitment*

The research team at Ulster University collaborated with rheumatologists from the Western Health and Social Care Trust (WHSCT) to design, conduct and recruit patients to the study. Informed consent to participate was obtained from all RA patients enrolled to the study. One hundred patients identified using following inclusion/exclusion criteria were recruited into the prospective observational cohort study: Remote Arthritis Disease Activity MonitoR (RADAR); ClinicalTrials.gov Identifier: NCT02809547. Inclusion criteria: aged between 18–90 years, diagnosed with RA (according to American College of Rheumatology criteria [5,16]), diagnosed with RA for a minimum of 1 year and maximum 10 year duration, active disease flares on a regular basis, and receiving a disease-modifying anti-rheumatic drug (DMARD). Exclusion criteria: any other inflammatory conditions, any infections or trauma during study period, and have restricted hand function (determined by clinical team). Office

for Research Ethics Committees Northern Ireland (ORECNI) (16/NI/0039), Ulster University Research Ethics Committee (UREC) (REC/16/0019) and WHSCT (WT/14/27) approvals were obtained for the study. Informed consent was obtained from all RA patients enrolled to the study.

#### *2.2. Whole Blood and Dried Blood Spot Sample Collection*

Venepuncture whole blood samples as part of the normal routine care pathway were forwarded to the hospital laboratories for multiple tests including CRP, ESR, Bilirubin, Liver enzymes and full blood count. An additional 5-mL EDTA tube of blood was collected from each of the 100 participants for DNA genotype and drug metabolite analyses within the RADAR study. Samples were collected at both study baseline and at a 6-week follow-up appointment at an outpatient rheumatology clinic for all participants. Additionally, a sub cohort of 30 of the above participants were supplied with a kit containing sufficient dried blood spot (DBS) cards [17], finger lancets and pre-paid and addressed postal envelopes with desiccant and biohazard sealable pouch to send a weekly samples (approximately 3–5 droplets of blood ~20 μL each) from home to the NICSM laboratory for drug metabolite analysis. Finger lancet blood droplets were deposited onto dried blood spot (DBS) Protein Saver 903TM cards (Whatman, GE Healthcare Life Sciences, Buckinghamshire, UK), pre-treated with a protein stabiliser coating.

#### *2.3. Nucleic Acid Isolation from Peripheral Blood*

Total DNA was isolated from peripheral blood samples using TRIzol reagent (TRIzol LS Reagent, Thermo Fisher Scientific, Basingstoke, UK. cat. No 10296-010) according to manufacturer's directions. Total DNA concentration was estimated by spectrophotometry (NanoVue Plus—GE Health Care, Buckinghamshire, UK).

#### *2.4. Mass-Spectrometry Analysis*

Determination of methotrexate (MTX) metabolites was performed using a liquid extraction surface analysis (LESA) coupled with nanoESI-triple quadruple mass spectrometer (QQQ) using Triversa nanomate (Advion, New York, NY, USA) and API 4000 QQQ Mass Spectrometer (AB Sciex, Cheshire, UK). Control MTX metabolites and internal standards were from Schircks Laboratories (Jona, Switzerland).

Quantitation for MTX and MTX metabolites was performed by the matrix-matched standards approach using an intensity ratio (ISTD/MTXs) calibration (10–2000 nM). Signal for each metabolite was the average of *n* = 2 (duplicate injection). A total of 5 nM was selected as LLOD (S/N ~ 3) for MTX, 7-OH MTX and MTX-PG2 and 8 nM was selected for MTX-PG3 to PG5; 10 nM was selected as LLOQ for MTX and all metabolites (S/N ~ 10). Intra- and inter-day precision was assessed at both 50 and 500 nM and coefficient of variation (CV) for MTX metabolites ranged from 2.0–7.2%. Linear regression coefficient (R2) of the back-calculated concentration against the nominal concentration for MTX and its metabolites was above 0.995.

Determination of sulfasalazine metabolites and teriflunomide, analysis was performed by liquid chromatography tandem mass spectrometry (LC-MS/MS) using a HP 1200 LC (Agilent, Palo Alto, CA, USA) and a Quattro micro mass spectrometer (Micromass, Manchester, UK). Control sulfasalazine metabolites, teriflunomide and internal standard were obtained from Sigma Aldrich, Gillingham, UK. Quantitation of sulfasalazine, its metabolites and teriflunomide was performed by the matrix-matched standards approach using an intensity ratio (ISTD/Analyte) calibration (5–500 μg/L). Signal for each metabolite was the average of *n* = 2 (duplicate injection). 5 μg/L was selected as LLOD (S/N > 5) for metabolites and 10 μg/L was selected as LLOQ for all metabolites (S/N > 10).

Intra- and inter-day precision was assessed at both 20 and 100 μg/L and coefficient of variation (CV) for sulfasalazine metabolites ranged between 1.4–5.8%. Linear regression coefficient (R2) of the back-calculated concentration against the nominal concentration for sulfasalazine, its metabolites and teriflunomide was above 0.992.

#### *J. Pers. Med.* **2020**, *10*, 149

Separation of targeted analytes was carried out by reverse phase chromatography using a C18 column in gradient mode. Quantitation of all analytes were performed in positive ion mode multiple reaction monitoring (MRM) using matrix-matched standards and stable isotope ratios. All Mass Spectrometry methods were validated according to ICH Guidance for selectivity/specificity, limit of detection/quantitation (LLOD/LLOQ), linearity and precision [18]

#### *2.5. Endpoint-Genotyping Using Taqman Assay*

Endpoint genotyping analysis was carried out using the LightCycler 480 real-time PCR system (ROCHE). The assay is based on the competition during annealing between probes detecting the wild type and the mutant allele. The 5 -exonuclease activity of DNA polymerase cleaves the doubly labelled Taqman probe hybridised to the SNP-containing sequence, once cleaved, the 5 -fluorophore is separated from a 3 -quencher. Two allele-specific probes carrying different fluorophores (VIC®, emission: 554 nm and FAMTM, emission: 518 nm) permits SNP determination in a single well without any post-PCR processing. Genotype is determined from the ratio of intensities of the two fluorescent probes at the end of amplification (endpoint instead of the entire cycle in conventional PCR).

#### *2.6. Taqman Probes Used for Single-Nucleotide Polymorphism Genotyping*

The concentration and integrity of the genomic DNA were assessed by microvolume spectrophotometer (NanoDrop 2000; Thermo Fisher Scientific). DNA samples were genotyped by the following TaqMan SNP genotyping assays [*MTHFR*-rs1476413, Assay Identification (ID) C\_8861304\_10; RS1-rs 2231142, Assay Identification (ID) C\_354526997\_10; *ABCC2*-rs3740065, C\_22271640\_10; *SLCO1B1*-rs4149081, Assay Identification (ID) C\_1901759\_20; *MTHFR*-rs4846051, Assay Identification (ID) C\_25763411\_10; *ABCB1*-rs10280623, Assay Identification (ID) C\_30537012\_10; *ATIC*-rs16853826, C\_33295728\_10; *MTHFR*-rs17421511, Assay Identification (ID) C\_32800189\_20; *ABCC2*-rs717620, C\_2814642\_10; *ABCC1*-rs246240, Assay Identification (ID) C\_1003698\_10; Life Technologies Ltd.).

#### *2.7. Validation of Polymorphisms by Pyrosequencing*

Validation of SNP genotyping results from the Taqman assays was performed on a subset of samples using pyrosequencing. Due to the high number of SNPs to cover, a method using a universal biotinylated primer was employed [19]. Briefly, this method involves the use of standard target specific primer pairs with a universal M13 sequence at the 5 end of one of the primers. A third, biotinylated M13-targeting primer is included in the PCR amplification reaction, leading to incorporation of biotin into the PCR product without the need for individual biotin labelling of each individual primer pair and thus lowering the cost of pyrosequencing considerably. A list of primers used for pyrosequencing are shown in Table S3.

PyroMark Assay Design Software 2.0 (Qiagen) was used for primer design in the SNP calling assay design format. PCR amplification was carried out using the Pyromark PCR kit (Qiagen) in 25 μL total volumes with 10–20 ng DNA and final concentrations of 0.2 mM for each primer. Standard PCR cycling conditions were used as per manufacturer's instructions and were consistent for all samples and targets. PCR products were checked by agarose gel electrophoresis and those with a positive single band of the expected size were taken forward into pyrosequencing on the Pyromark Q48 (Qiagen) using standard manufacturers protocols and the instrument run setting for SNP calling.

#### *2.8. Statistics and Sample Size Calculations*

Statistical analysis was carried out with SPSS ver.25 (IBM Corp, NY, USA), SciPy module (ver. 1.3) for Python (version 3.7.2) and R (version 3.60) with *p* < 0.05 considered as statistically significant, all within 95% confidence intervals. Descriptive statistics were used to characterise the variability in mean MTX, MTX-7-OH and MTX2PG–5PG and MTXtotal concentrations between different genotype groups of patients. Normality of data was determined using the Shapiro-Wilks test in SPSS (ver. 25) prior to employing the Kruskal–Wallis non-parametric (distribution free) one-way ANOVA, with

Dunn–Bonferroni post hoc test to assess differences between genotype group means using GraphPad Prism version 8.0.0 for Windows (GraphPad Software, San Diego, CA, USA, www.graphpad.com).

The Hardy-Weinberg equilibrium was assessed for SNPs with significant clinical associations in the methotrexate treated cohort. A Chi-square test with Benjamini-Hochberg adjusted *p* values was used to assess if there were significant differences between the genotype frequencies expected from dbGAP European population and those observed. Power was calculated for the same SNPs using the GENPWR package [20] within R (v 4.0.2), using the linear regression model (with alpha at 0.05) since the goal was to calculate power in a continuous outcome (metabolite levels) between genotypes.

#### **3. Results**

#### *3.1. Study Population*

A total of 100 participants, *n* = 68 female and *n* = 32 male, with active rheumatoid arthritis were enrolled to this study (Table 1). The mean age of study participants was 59.5 years with a mean disease duration of 6 years and mean baseline disease activity score (DAS28ESR) of 3.6.

**Table 1.** RADAR Study Cohort Demographics. Clinical and laboratory feature summary across *n* = 100 participants. yr: years; Anti-CCP: anti-cyclic citrullinated peptide; ESR: erythrocyte sedimentation rate; DAS28: disease activity score across 28 joints; RBC, red blood cell (count); Hb: haemoglobin; WBC: white blood cells; ALT: alanine aminotransferase; AST: aspartate aminotransferase; ALP: alkaline phosphatase; SD: Standard deviation.


A subgroup of *n* = 66 participants (*n* = 46 female) were treated with weekly methotrexate at baseline was identified for subsequent genotype association analyses. Baseline and 6-week follow-up drug dose information is summarized in Table S2A,B for this main subgroup. A smaller, partially overlapping, subgroup of *n* = 27 participants (*n* = 20 female) being treated with daily sulfasalazine at baseline was also identified for subsequent analyses (Table S2B).

#### *3.2. Single-Nucleotide Polymorphisms Analysed in this Study*

Ten SNPs were analysed in this study, previous studies have linked these SNPs to various clinical consequences observed in RA patients being treated with DMARDs, documented in the PharmGKB database [21,22]. SNP genotypes were determined by endpoint PCR assay and allele specific probes (Figure 1). The SNPs characteristics and frequencies in the study cohort are summarised in Table 2.



*J. Pers. Med.* **2020** , *10*, 149

**Figure 1.** Endpoint polymorphism genotype assay. Genotype were determined for the listed polymorphisms in the RADAR study cohort (*n* = 100) from the ratio of fluorescent intensities (nm) of the two-allele specific Taqman probes (VIC and FAM) at the end of PCR amplification. Clusters in upper and lower quadrants represent groups of individuals with a homozygous genotype for either allele; the middle quadrant represents individuals with heterozygous genotype. (**A**) rs246240, (**B**) rs717620, (**C**) rs1476413, (**D**) rs17421511, (**E**) rs2231142. (**F**) rs3740065, (**G**) rs4149081, (**H**) rs4846051, (**I**) rs10280623, (**J**) rs16853826. Genotypes were also confirmed by pyrosequencing in selected individuals. Genotype frequency is summarized in Table 1.

#### *3.3. Methotrexate and Sulfasalazine Metabolite Polymorphism Associations*

The data strongly suggest plasma concentrations of methotrexate and sulfasalzine metabolites are associated with the allelic genotype for 2 particular polymorphisms, rs4149081 and rs1476413 (see Figure 2). Table 3 indicates that within the *n* = 66 methotrexate treated subgroup, *n* = 17 participants with the minor homozygote genotype AA in rs4149081 have a significantly lower mean plasma MTX-7-OH concentration compared to the GA genotype group (*p* = 0.002) at baseline. Although a similar trend is observed at the 6-week follow-up, this was not statistically significant. The GA genotype group mean concentrations of MTX and MTX-7-OH are significantly higher than those observed in the GG genotype group (*p* = 0.01 and *p* = 0.038, respectively; Figure 2A,B). The baseline mean blood bilirubin concentration was the only feature observed at significantly higher levels in the rs4149081 AA genotype, relative to the GG genotype group (*p* = 0.020; Figure 2C).

A total of *n* = 8 participants with the rs1476413 homozygous major allele genotype CC have significantly lower (*p* = 0.012) group mean plasma MTX-7-OH concentration, compared to the CT genotype group (Table 3) at the 6-week follow-up sessions. No significant difference was observed at baseline. The mean plasma concentration of tetraglutamate MTX metabolites are also significantly lower in the rs1476413 CC genotype group at both baseline (*p* = 0.02) and 6-week follow-up appointments (*p* = 0.008; Figure 2F,G and Table 3).

A total of n = 6 participants with the minor allele genotype AA in the SNP rs17421511 (Supplementary Figure S1A,B) show significantly higher mean plasma concentration (*p* = 0.013) of sulfapyridine at the 6-week follow-up appointment period only, relative to GG and GA genotype groups.

**Figure 2.** rs4149081 and rs1476413 genotype associations. Statistically significant associations between two polymorphisms, plasma drug metabolite concentration, other laboratory and clinical outcome measures are shown for individual genotypes in methotrexate treated participants (*n* = 66). Each symbol represents an individual participant of the genotype indicated on the x axes. Data grouped by rs4149081 genotypes: (**A**) baseline plasma concentration of unmetabolised methotrexate, (**B**) baseline plasma concentration of 7-hydroxy-methotrexate, (**C**) baseline bilirubin blood concentration, (**D**) weekly plasma concentrations (log scale) of listed methotrexate metabolites (PGs: polyglutamate subtypes) of a rs4149081 GA genotype participant. Data grouped by rs1476413 genotypes: (**E**) 6-week follow-up plasma concentration of 7-hydroxy-methotrexate, (**F**) baseline and (**G**) 6-week follow-up plasma concentrations of long-chain methotrexate 4-glutamate, (**H**) baseline and (**I**) 6-week follow-up disease activity (DAS28ESR) scores. Statistically significant differences between genotype group means are indicated by horizontal bars and an asterisk used to summarise *p* values adjusted by Bonferroni's multiple comparison test: (\*) *p* < 0.05; (\*\*) *p* < 0.005 (descriptive statistics data shown in Table 3). Red horizontal bar represents genotype group mean; error bars represent standard deviation. MTX: methotrexate; BL: baseline; 6w: 6-week follow-up.



#### *3.4. Clinical and Laboratory Feature Polymorphism Associations*

While the average red blood cell counts (RBC) remain within reference ranges in both men and women in this study (Table 1), there is a modest but statistically significant decrease in mean RBC levels in the *n* = 15 participants with the rs2231142 heterozygous genotype GT compared to the TT genotype group at both study time points (*p* = 0.044; Table 3, Figure 3A,B). Mean alkaline phosphatase concentration is also significantly lower at baseline in the rs2231142 GT genotype group, relative to the TT genotype participants (*p* = 0.019; Figure 3C). In the smaller group of *n* = 6 rs2231142 GG genotype participants, mean lymphocyte counts are significantly higher than the GT genotype group, again at both time points (*p* = 0.043, *p* = 0.019; Figure 3D,E).

**Figure 3.** rs2231142 and rs17421511 genotype associations. Statistically significant associations between two polymorphisms and other laboratory and clinical outcome measures are shown for individual genotypes in methotrexate treated participants (*n* = 66). Each symbol represents an individual participant of the genotype indicated on the x axes. Data grouped by rs2231142 genotypes: (**A**) baseline and (**B**) 6-week follow-up red blood cell (RBC) count, (**C**) baseline blood alkaline phosphatase (ALP) concentration, (**D**) baseline and (**E**) 6-week follow-up lymphocyte count. Data grouped by rs17421511 genotypes: (**F**) baseline and (**G**) 6-week follow-up patient assessed pain (PgPain) levels, (**H**) baseline disease activity (DAS28ESR) scores, (**I**) baseline blood alanine aminotransferase (ALT) concentration. Statistically significant differences between genotype group means are indicated by horizontal bars and an asterisk used to summarise *p* values adjusted by Bonferroni's multiple comparison test: (\*) *p* < 0.05; (\*\*) *p* < 0.005 (descriptive statistics data shown in Table 3). Red horizontal bar represents genotype group mean; error bars represent standard deviation. BL: baseline; 6w: 6 weeks follow-up.

Mean patient-reported general pain (PgPain) scores are significantly lower in *n* = 9 participants carrying the rs17421511 AA genotype at baseline (*p* = 0.033) and at the 6-week follow-up appointment (*p* = 0.013) compared to those with the GA genotype (Figure 3F,G). This trend is also reflected in significantly lower mean baseline DAS28ESR scores in participants carrying the rs17421511 AA genotype, relative to the GA and GG genotype groups (*p* = 0.002 and *p* = 0.005 respectively; Table 3, Figure 3H), though scores even-out at the 6-week follow-up period among all three genotypes. The baseline mean alanine aminotransferase (ALT) was recorded at significantly higher blood concentrations in the rs17421511 AA genotype group, relative to the GG genotype participants (*p* = 0.048; Figure 3I).

The Benjamini–Hochberg-adjusted Chi-square test p-values showed no statistically significant differences between the genotype frequencies observed and those expected from dbGAP European population: rs4149081, p*adj* = 0.515; rs 1476413, p*adj* = 0.945; rs17421511, p*adj* = 0.711; rs2231142, p*adj* = 0.9454. The power calculated for each of these SNPs at *p* < 0.05 was: rs4149081, 0.96; rs 1476413, 0.95; rs17421511, 0.97; rs2231142, 0.94.

#### **4. Discussion**

This study investigates the influence of ten well-characterised SNPs in RA and we have tried to correlate this with the appearance and accumulation of metabolites measured in the plasma of patients taking DMARDs such as methotrexate and or sulfasalazine. The in vivo pharmacotherapy of DMARDS and potential response biomarkers in RA have been previously described [6,23–25], however there studies of potential associations between circulating csDMARD levels and specific genetic variants remain limited in RA patients.

Typically, methotrexate treatment may cause elevations in serum AST and ALT, long term therapy has also been linked to development of fatty liver disease, fibrosis, cirrhosis, nephrotoxicity, and renal failure [26]. However, under active consultant-led clinical management, these effects are largely minimised. The mean values of all clinical biomarkers, liver enzyme and blood component cell-counts (Table 1) are within recommended normal reference ranges when viewed across all of the study participants. However, mean ALT was significantly higher at baseline in the methotrexate treated subgroup of participants with the rs17421511 AA genotype, albeit the potential effect of multiple drug combinations was not investigated in this subgroup.

Although the average RBC count when taken from all participants appear to be within normal ranges (Table 1), participants with the GT allele in the rs2231142 SNP have significantly decreased erythrocyte counts in their circulation compared to those with the homozygous alleles methotrexate affects folic acid metabolism, thus patients taking MTX may show variations in their mean corpuscular volume (MCV) of red blood cells (RBC), therefore resulting in megaloblastic anaemia.

RBCs retain MTX as the polyglutamate derivatives throughout their lifespan [27,28]. While normal RBC levels are between 4.7 to 6.1 million cells per microlitre (mill.c/μL) for men and between 4.2 to 5.4 mill.c/μL in women, the slight decrease shown in heterozygous rs2231142 SNP patients is statistically significant. However, the lower mean haemoglobin levels observed in the rs2231142 GT genotype is not statistically significant and there is no correlation with disease activity score as may have been anticipated in anaemia of chronic disease.

Apart from the impact of sex-linked genes in RA, the diversity in our genomes are partially accountable for the heterogeneity in the clinical presentation of synovitis among patients [29]. The genetic influence in RA is particularly strong, the heritability in RA is estimated to be around 60% [29] and with the high diversity of clinical presentations observed in RA, the goal in treatment would be to stratify patients according to their genetic profile and clinical outcome, eventually formulating a genetic-risk based personalised treatment management protocol.

MTX is an anti-folate drug, with anti-proliferative and anti-inflammatory effects, by inhibition of folate and adenosine pathways and inhibition of purines and pyrimidines synthesis [30–32]. Approximately 80–90% of methotrexate is primarily excreted by the kidneys [33]. MTX is converted in hepatic parenchymal cells of some patients resulting in the 2- through 4-glutamate residues derivatives or the drug is catabolised to the MTX-7-OH form.

Though not observed consistently on both study time points, participants with the AA allele in rs4149081 and CC allele in rs1476413 can have significantly lower mean plasma levels of MTX-7-OH in their plasma circulation. Since some genotype groups are modest in size, the potential for differences in the mean MTX dose between genotype groups was analysed, though no significant difference was observed (Table S2). Furthermore, only a weak correlation exists between MTX dose and circulating MTX-7-OH (r2 = 0.08213). Increasing levels of MTX-7-OH is known to inhibit the clinical responsiveness of RA patients to the MTX drug and therefore, reduced levels of this metabolite could signify a better clinical response to MTX [34]. Thus, with genetic profiling of expanded csDMARD naïve RA cohorts, it would be interesting to further investigate clinical responsiveness to MTX in these genotypes.

The 2- through 4-polyglutamate MTX metabolite-derivatives are selectively retained in cells and participants with the TT or CT alleles in the rs1476413 SNP tend to show significantly higher mean plasma levels of the tetraglutamate metabolite. A significantly higher mean DAS28 is observed only at study baseline in the rs1476413 CT genotype group relative to CC genotype, though due to modest numbers in the latter group this would require independent verification. It is likely that folic acid supplementation in the study cohort to mitigate toxicity of MTX has reduced the frequency of observable side effects.

Sulfasalazine is metabolized by intestinal bacteria, resulting in the release of sulfapyridine (SPY) and 5-aminosalicylate or 5-ASA (SPY and 5ASA are linked by an azo bond) [35]. Sulfapyridine is almost completely absorbed by the colon, metabolized by the liver, and renally excreted [1]. Commonly reported side-effects of sulfapyridine are minor gastrointestinal (GI) and central nervous system (CNS) abnormalities, and uncommon serious haematological and hepatic side-effects [36,37]. Although study participants with the AA and GA alleles of the rs17421511 SNP indicate higher mean plasma levels of sulfapyridine compared to those with the GG allele (Figure S1), no significant adverse phenotypic effects were observed in these subgroups.

The modest number of patients with the AA genotype of the rs17421511 SNP in our study report significantly lower levels of pain and disease activity, relative to the remaining methotrexate treated cohort. In future research with expanded patient cohorts, it would be pertinent to see if this phenomenon is observed in other patient groups carrying this particular genotype. As a general observation, a limitation of the current study is the low number of participants in particular genotype groups and the smoking status was not recorded, which may impact upon methotrexate metabolism. Furthermore, the findings for the methotrexate treated cohort are only generalizable to the European population, as no significant differences in Hardy–Weinberg equilibrium were found by Chi-square test between the dbGAP frequencies and those observed in this study.

While it is challenging to find a clear-cut relationship between genotype and circulating drug levels which translates through to a clear prediction of phenotypic consequence, useful leads are presented in the current study. The rs1476413 and rs17421511 *MTHFR* variants and the rs2231142 *ABCG2* variant display significant changes which are consistent at both study time points.

With further carefully powered studies of variability in both csDMARD response and predisposition to side effects, there is considerable potential to personalise effective treatments whilst avoiding any toxicity.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2075-4426/10/4/149/s1. Table S1: Methotrexate treated cohort SNP association with laboratory and clinical features; Table S2A: Methotrexate treated cohort prescription data; Table S2B: Methotrexate and sulfasalazine treated subgroup co-prescription data; Table S3: Pyrosequencing SNP primers; Figure S1: Sulfasalzine treated cohort SNP associations.

**Author Contributions:** Conceptualization, D.S.G.; methodology, L.G.D., K.G.M., C.C., D.C., K.B.C.T.; patient recruitment and data collection, P.V.G., P.C.; formal analysis, L.G.D.; data curation, K.G.M.; writing—original draft preparation, D.S.G., L.G.D.; writing—review and editing, P.S., P.V.G., D.C., C.C.; visualization, D.S.G.; supervision, D.S.G.; project administration, D.S.G.; funding acquisition, D.S.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Research and Development Office, Health and Social Care NI (RADAR study), The Northern Ireland Rheumatism Trust, and Invest Northern Ireland (PoC 631/731) in support of this work.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **Modelling Neuromuscular Diseases in the Age of Precision Medicine**

#### **Alfina A. Speciale, Ruth Ellerington, Thomas Goedert and Carlo Rinaldi \***

Department of Paediatrics, University of Oxford, Oxford OX1 3QX, UK; ambra.speciale@paediatrics.ox.ac.uk (A.A.S.); ruth.ellerington@paediatrics.ox.ac.uk (R.E.); tg5g18@soton.ac.uk (T.G.)

**\*** Correspondence: carlo.rinaldi@paediatrics.ox.ac.uk; Tel.: +44-(0)1865-272148; Fax: +44-(0)1865-282840

Received: 23 September 2020; Accepted: 15 October 2020; Published: 17 October 2020

**Abstract:** Advances in knowledge resulting from the sequencing of the human genome, coupled with technological developments and a deeper understanding of disease mechanisms of pathogenesis are paving the way for a growing role of precision medicine in the treatment of a number of human conditions. The goal of precision medicine is to identify and deliver effective therapeutic approaches based on patients' genetic, environmental, and lifestyle factors. With the exception of cancer, neurological diseases provide the most promising opportunity to achieve treatment personalisation, mainly because of accelerated progress in gene discovery, deep clinical phenotyping, and biomarker availability. Developing reproducible, predictable and reliable disease models will be key to the rapid delivery of the anticipated benefits of precision medicine. Here we summarize the current state of the art of preclinical models for neuromuscular diseases, with particular focus on their use and limitations to predict safety and efficacy treatment outcomes in clinical trials.

**Keywords:** neuromuscular diseases; translational research; disease models; precision medicine

#### **1. Introduction**

Neuromuscular diseases are a broad and heterogeneous group of conditions characterized by an impairment in one or more components of the motor unit, defined as the motor neuron and the muscle fibres it innervates. Whilst most are individually rare, collectively neuromuscular diseases are significantly prevalent, with a cumulative prevalence of approximately 100–200 cases per 100,000 individuals worldwide [1], accounting for a substantial proportion of population-wide health care costs [2]. Very few treatments currently exist to treat these diseases. Nevertheless, as research progressively disentangles their pathogenic mechanisms, many opportunities are finally starting to land in the clinic.

Precision medicine refers to a treatment approach wherein the most appropriate treatment for an individual is chosen based on their specific disease manifestation, alongside their genetic/epigenetic information and other features such as their microbiome, age, nutrition, and lifestyle. The clinical and genetic heterogeneity of neuromuscular diseases make them ideal candidates for personalized therapeutic approaches, with many individuals suffering from rare or ultrarare diseases that cannot be treated by conventional blanket approach treatment. One example is Duchenne muscular dystrophy (DMD), the most prevalent childhood-onset muscular dystrophy, where progressive muscle degeneration and weakness is caused by mutations in the *DMD* gene, leading to loss of dystrophin protein production [3]. The vast majority of DMD patients carry an exon deletion (~65%) or a duplication (~10%) of one or multiple exons and these mutations tend to manifest in regions of vulnerability between exons 2 and 20 and exons 45 and 55 [4–6]. In addition, small mutations (insertions, deletions, nonsense mutations and splice site mutations) account for the remaining ~25% mutations and occur throughout the length of the gene [4]. Excision of specific exons, or exon skipping, by use of

antisense oligonucleotides (AON) to allow restoration of the disrupted reading frame and therefore production of a shortened but functional dystrophin protein, has surfaced as a promising therapy for DMD [7]. Therefore, diagnosis by genetic sequencing has become a crucial tool in determining eligibility for these treatments, as multiple AON products need to address the large series of mutations carried by DMD subjects.

While presenting new challenges for researchers, precision medicine is rapidly taking the lead in the pursuit of radically transforming health care. Choosing the appropriate disease model that recapitulates the complexity and heterogeneity of patients is therefore paramount to understand disease mechanisms and increase the chances of success of translating a treatment opportunity into a safe and effective marketed drug.

In this review, we aim to discuss the currently available tools used to model neuromuscular diseases and to evaluate their utility and applicability to personalized medical research and therapeutic development (Table 1).

#### **2. Cellular Models**

#### *2.1. Myoblasts*

Primary myoblasts (activated satellite cells) obtained from human subjects or animal models typically go through multiple rounds of cell division until reaching confluence in growth media, followed by iterations of cellular fusions to form multinuclear myotubes and eventually terminal differentiation [8]. Due to several inherent traits of human-derived muscle cells, including the slower growth rate as well as the flattened morphology, primary human myotubes typically exhibit poorer contractile activity than their mouse counterparts in response to electric stimulation [9]. Obtaining a substantial number of satellite cells from skeletal muscle biopsies of patients is markedly limited by the restricted proliferative capability of activated satellite cells in culture. In order to overcome this limitation, myogenic conversion of non-muscle primary cells, such as primary human and murine fibroblasts from skin, has been widely employed, mainly using transduction of *MyoD* gene (myogenic differentiation), a master regulator of skeletal muscle differentiation [10]. In order to increase proliferative capacity, transduction with both telomerase-expressing and cyclin-dependent kinase 4-expressing vectors has been used to produce immortalized human muscle stem-cell lines from patients with different muscle diseases such as DMD, limb-girdle muscular dystrophy type 2B, facioscapulohumeral muscular dystrophy, oculopharyngeal muscular dystrophy and congenital muscular dystrophy [11]. These immortalized cultures have been extensively used both to study disease mechanism and to test treatment strategies.

#### *2.2. Induced Pluripotent Stem Cells (iPSCs)*

The development of induced pluripotent stem cell (iPSC) technology has brought a great paradigm shift in the field of precision medicine [12] and now they have a prominent role as a tool for disease modelling and drug screening. Moreover, they are highly expandable, are free from the ethical issues linked to the use of embryonic stem cells (ESCs), and their source of cells easily accessible.

Two major strategies have been recently developed to differentiate PSCs into satellite-like cells. The first involved overexpressing PAX7, the master transcription factor for satellite cells, in an inducible fashion [13]. After being generated from human embryonic stem cells and iPSCs, these cells showed capability for in vitro expansion and differentiation, as well as engraftment and myofibre formation in immunodeficient mice [13,14]. The second strategy involved the use of a small molecule, and consists of glycogen synthase kinase 3 beta (GSK3beta) inhibition, in order to activate the Wnt pathway, as well as treatment with fibroblast growth factor 2 (FGF2) in a minimal medium [15–20]. Alternative protocols have used bone morphogenic protein 4 (BMP4) inhibition to promote differentiation into the myogenic lineage [21–23], or Notch signalling inhibitor DAPT [24]. Purified by fluorescence-activated cell sorting (FACS) [15,19,24], partially purified, or unpurified [16,17,20,21,23], cell mixtures are then plated.


Keyfeaturesofthevariousmodelsusedforneuromuscular

*J. Pers. Med.* **2020**, *10*, 178

By generating an in vitro DMD model from patient-derived iPS cells, Shoji et al. noted excess Ca2<sup>+</sup> influx in DMD myocytes when compared to control myocytes in response to stimulation via electricity. This was alleviated by restoring dystrophin expression via exon skipping, therefore establishing a model that recapitulates early DMD pathogenesis and is appropriate for assessing the efficacy of exon-skipping drugs by phenotypic assay [25]. IPSC models of several other neuromuscular diseases are currently available, including Miyoshi myopathy, a muscle disease caused by the mutation in dysferlin [26], Pompe disease, a paediatric disease caused by lysosomal glycogen accumulation in skeletal muscle that leads to muscle weakness [27], and myotonic dystrophy type 1, a multisystem disorder that affects skeletal and smooth muscle caused by a CTG trinucleotide repeat expansion in the non-coding region of the *DMPK* gene [28]. Overall, the introduction of iPSC technology has allowed scientists to model diseases directly from patients' cells, this being a cornerstone for personalized medicine. However, if they are planned to be used for personalized cell therapy, several issues remain to be addressed, including alterations in the differentiation efficiency, line-to-line variability, and risk of tumorigenicity.

#### *2.3. Urine-Derived Stem Cells*

In addition to representing an ideal source of cells for generating iPSCs, with a reprogramming efficiency approximately 100-fold higher than that of fibroblasts [29], urine stem cells (USCs) can also be induced into myogenic lineage by direct MyoD1 reprogramming [30]. Muscle differentiation can be further enhanced by adding 3-deazaneplanocin A hydrochloride [31]. These cells carry pluripotency markers such as CD29, CD105, CD166, CD90, and CD13 [32], and are able to self-renew and differentiate into the mesodermal, endodermal and ectodermal lineage [33]. Direct reprogramming of these cells, which can be easily isolated by centrifugation method and standard cell culture, has been recently shown to efficiently and reproducibly establish human myogenic cells from patients with DMD and limb-girdle muscular dystrophy (LGMD) type 2 [30]. Upon further molecular characterisation, this cost-effective and efficient in vitro model system shows great potential for more efficient drug development and targeted therapies development for neuromuscular diseases.

#### *2.4. Skeletal Muscle Organoids*

As the use of human iPSCs for tissue engineering and disease modelling expands, iPSC-derived organoids are rapidly becoming a powerful tool for modelling human organogenesis, homeostasis, injury repair and disease aetiology [34]. These miniature 3D tissues are generated using a combination of signposted differentiation, morphogenetic processes, and the embryonic organogenesis mimicking intrinsically driven self-assembly of cells, resulting in architecture and function remarkably similar to their in vivo counterparts. By using natural or synthetic scaffolds to create the artificial tissue [35], these models account for the cell–cell and cell–extracellular matrix interactions as well as the mechanical and/or chemical cues [36,37]. The development of physiologically relevant 3D in vitro models holds great promise to provide more economic, scalable and reproducible means of testing drugs and therapies for successful clinical translation. Few studies have reported methods to engineer human skeletal muscle tissue [38–43]. Induced myogenic progenitor cells derived from multiple human iPSC lines have been shown to form functional skeletal muscle tissues and are able to survive, progressively vascularize, and maintain functionality when implanted into the hindlimb muscle or dorsal window chamber in immunocompromised mice [44]. Isogenic human iPSC-derived 3D artificial muscles from patients affected by DMD, limb-girdle type 2D, and lamin A/C (LMNA)-related muscular dystrophies have been recently generated, recapitulating several pathogenic hallmarks in these diseases and also showing potential for muscle engraftment [45]. These studies have indicated that generation of fully functional artificial muscles require the contribution from other cellular lineages, for example vascular cells and motor neurons [45–49]. The major challenges the field is currently facing are mainly related to improving organoids' scalability as well as their complexity and maturity. Recent success in growing brain organoids using multiwell spinning bioreactors represents a significant step towards

high-throughput drug screening via large-scale organoid generation [50]. These models resemble more closely foetal than adult tissue, therefore optimisation of protocols is essential before being able to advance these tissues into replacement therapy. Bearing in mind the speed at which the field has advanced over the past few years, the range of possible future applications of this platform in the study of human diseases and in regenerative medicine is expected to rapidly expand.

#### *2.5. Muscle on Chip*

Advancement in culturing models with mixed culture capabilities, together with the latest developments in 3D printing, microfluidics and microfabrication engineering, has led to the rapid expansion of organ-on-chip technologies. These platforms have recently attracted substantial interest due to their potential to be informative at multiple stages of the drug discovery process, while offering new ways to model disease states and perform mechanistic investigations in vitro. The critical and defining features of these platforms are the 3D structure, the possibility of integration of multiple cell types to reflect tissue physiology, and the presence of relevant biomechanical forces [51]. Organ on chips have been adapted for the human gut [52], heart [53], blood–brain barrier [54], and kidney [55]. Human primary myogenic cells have been engineered to form 3D myobundles, which respond to electrical stimuli and undergo dose-dependent hypertrophy or myopathy in response to pharmacological stimulation [40]. The decreased muscle regeneration capacity and weakness observed in DMD patients have been recapitulated in a human dystrophic skeletal muscle on a chip [56]. Using a 3D photo-patterning approach, other researchers have developed a skeletal muscle platform by confining a cell-laden gelatin network around two hydrogel pillars, which serve as anchoring sites for the cells, as the muscle tissues form and mature [57]. In other instances, neurons and rhabdomyocytes, both originating from mouse embryonic cells, have been differentiated in a 3D hydrogel culture, to effectively constitute a neuromuscular unit on a chip [58].

Tissue engineering requires a deep understanding of the functional interplay of cell types and the effect of the scaffold on cellular architecture, as well as careful characterising and validation of the model for the purpose of study. Additionally, due to safety concerns around the potential for unexpected toxic side effects, the biocompatibility of the materials to be used must be well profiled [51].

As iPSCs or adult stem cells taken from mass production of tissue organoids are increasingly employed as a source of cells for these platforms, organ on a chip represents an ideal tool for precision medicine.

#### *2.6. Other*

Sources in addition to the muscle-derived cells or reprogrammed cells can be employed to model muscle diseases. For example, melanocytes from DMD patients show the same morphological alterations as DMD muscle-derived cells [59]. Cultured melanocytes from skin biopsies have been shown to be a useful alternative to muscle biopsies for the mRNA-based molecular diagnosis of DMD [60]. Additionally, in the case of Ullrich congenital muscular dystrophy (UCMD) and Bethlem myopathy (BM), diseases caused by mutations in collagen VI genes [61], patients' derived melanocytes recapitulated the mitochondrial dysfunction and ultrastructural alterations that are found in patient myoblasts [62].

#### **3. Animal Models**

#### *3.1. Mouse Models*

A large fraction of currently available therapies have been developed with the help of animal models, especially mice, mainly due to the high similarity in sequence homology and organ physiology to humans, as well as cost-effective husbandry. Additionally, the external environment in mice studies can be well controlled and monitored and studies using inbred mice allow resampling isogenic individuals, therefore minimising variability.

Nevertheless, many differences remain: mice are smaller in size, have a markedly reduced lifespan and an increased heart rate, just to name a few. Approximately 1% of human genes are not present in the mouse genome [63], while the differences in the promoter regions, non-coding sequences, and RNA splicing might be even more marked, accounting for species-specific disparities in gene expression that in some cases can affect disease phenotype [64,65]. Overall these considerations, together with the realisation that treatments in mice have frequently resulted in disappointing outcomes in clinical trials, have recently called into question the translational potential of findings in mouse models [66].

One way of making mouse models for studying human diseases more suitable is to follow approaches pioneered over 30 years ago, which comprise incorporating human DNA into the mouse genome (genetic humanisation) and/or engrafting human cells and tissue into mouse tissues (cellular humanisation) [67–70]. Genetic humanisation can be achieved through a variety of methods, most commonly by injection of plasmids or artificial chromosome vectors into the mouse zygotes. Transgenic models have substantially contributed to advancing the understanding of human disease and have helped develop treatment strategies. One notorious major breakthrough in biomedical research using transgenic mice carrying the human *SMN2* gene led to the recent clinical approval of an AON, able to block an intronic splicing silencer in human *SMN2* [71], increasing full-length *SMN2* isoform expression, which compensates for the loss of *SMN1* that causes spinal muscular atrophy [72–75].

However, some key features must be considered: the cDNA or genomic DNA used to generate the transgenic mice tend to integrate randomly in multiple copies and thus overexpress the protein of interest. Overexpression of wild-type proteins may give a dose-dependent phenotype not related to the disease mutation, like in the case of the androgen receptor [76], and RNA binding proteins, such as TAR DNA-binding protein 43 [77]. The rise of genome engineering technology has revolutionized the field of molecular biology by allowing the generation of physiological, humanized knock-in mice models by precise editing [78,79]. Most DMD preclinical studies have been carried out in the mdx mouse that carries a nonsense point mutation in DMD exon 23 [80], which is only one out of the thousands of possible variations in this gene present in DMD patients. Despite a lack of dystrophin expression, these mice do not exhibit dilated cardiomyopathy or a shortened lifespan. To improve upon this model, a number of double knock-out mouse models have been created, such as mice deficient in both dystrophin and its homolog utrophin, which show decreased cardiac function and survival [81]. In recent years by using clustered regularly interspaced short palindromic repeat (CRISPR)-based editing, many new DMD mouse models carrying deletions, frameshifting mutations, a point mutation, and a mutant version of the human *DMD* gene have been generated [82–88], making testing of exon skipping strategies targeting different parts of the DMD transcript possible. It is worth considering that recent studies to assess the effects of disease-causing mutations or environmental stimuli in different mouse strains found a strong influence of the genetic background on phenotypic responses [89], highlighting the importance of genetic diversity of animal models in biomedical research.

It is becoming more and more evident that choosing the right model is critical. Depending on the specific research question, often combining different strains is the most appropriate way to minimize the risks of a lack of reproducibility of translational research. Despite the obvious differences between mice and humans, genetic mouse models have allowed us to look at the effects of a mutation at a system level. Combining genetic engineering, which has made genetic modifications of endogenous targets possible, with the use of genetic with cellular humanisation, we now have powerful tools to study human pathophysiology in vivo, in cell-autonomous and non-cell-autonomous contexts [90], as well as excellent preclinical models to identify and test the pharmacodynamic and pharmacokinetic properties of a treatment strategy, from gene therapy to small-molecule and cell replacement [91]. Overall, these considerations further support the use of 'mouse precision medicine' as a better prototype for future mouse studies.

#### *3.2. Drosophila Melanogaster*

Drosophila melanogaster can serve as a useful model of human neuromuscular disease, since flies have a neural circuitry, albeit much simpler than in humans, as well as multinucleated muscle cells and neuromuscular junctions (NMJ). The mechanisms of synaptic transmission seen at the NMJ in humans are conserved in Drosophila, with a key difference being that Drosophila uses glutamate, not acetylcholine, as the neurotransmitter. The ability to genetically manipulate Drosophila is useful when trying to better understand how certain myopathies occur. Moreover, their short life span and large progeny make flies a good system for carrying out large-scale genetic screens. Drosophila has helped us understand more about the NMJ, and in particular, the role that the dystrophin–glycoprotein complex plays (DGC). Like in mammals, the Drosophila gene of dystrophin also encodes multiple isoforms, which contain highly conserved domains and are mainly expressed in the muscle and the nervous system [92–94]. Studies into DGC function at the NMJ of Drosophila have shown that it plays an important role in the retrograde control of neurotransmitter release, neuronal migration and muscle stability and thus may help explain how neuromuscular pathology can occur. Removal of a dystrophin isoform (DLP2) in Drosophila, which is normally located at the post-synapse, has been shown to lead to an increase in presynaptic neurotransmitter release, causing increased muscle depolarisation, thus indicating a role of dystrophin in regulating presynaptic neurotransmitter release [95]. Previous work has shown that by studying sensory neurons (photoreceptor cells) in Drosophila [96], a lot can be learnt about axon guidance and target recognition. Perturbation of dystrophin and dystroglycan in photoreceptor cells led to disrupted axon guidance, similar to neuronal defects seen in human muscular dystrophy patients. Drosophila not only aids us in understanding the role that certain proteins play at the synapse of the NMJ, but also serves as a good model for studying age-dependent progression of muscular dystrophy. The reduction in levels of expression of dystrophin isoforms in Drosophila using RNAi led to muscle degeneration in larval and adult flies [95], thus potentially providing a useful model to help us understand Duchenne muscular dystrophy pathogenesis in humans.

#### *3.3. Zebrafish*

The zebrafish (Danio rerio) has become a useful organism for studying neuromuscular genetic disorders [97]. Comparison to the human reference genome has shown that approximately 70% of human genes have at least one zebrafish orthologue [98], and dozens of mutant zebrafish lines have already been generated to model the most common human myopathies [99–101]. As vertebrates, they possess desirable attributes, including small size, rapid development, and genetic tractability [97]. Zebrafish embryos are transparent, develop externally and can be easily genetically manipulated [102], making this model ideal for phenotypic high-throughput screening platform to investigate drug efficacy in a whole-organism context. The most commonly adopted screening criteria for assessing neuromuscular phenotype are spontaneous coiling, ability to hatch on time, swimming behaviour, and birefringence assay [103]. Compared to target-based drug discovery, a phenotype-driven approach offers several key advantages [104], such as rapid identification of compounds that have poor bioavailability, exhibit toxicity or off-target effects. By screening small-molecule libraries in the dystrophin-null zebrafish (sapje model), aminophylline, a non-selective phosphodiesterase inhibitor, was found to improve survival rate in animals, restore normal muscle structure and up-regulate the cAMP-dependent PKA pathway without affecting dystrophin expression [105]. In the sapje model, the mitochondrial defects present in DMD patients were recapitulated, making it an optimal model for the disease, and it was used to assess the effect of the cyclophilin inhibitor alisporivir treatment in vivo, resulting in an improvement in the morphology of mitochondria and myofibrils, and in mitochondrial respiration [106]. A zebrafish model showing severe myopathy has also been generated for UCMD via a deletion in the *col6a1* gene through the injection of an antisense morpholino [107]. Here, defects in the mitochondria permeability transition pore (mPTP) were corrected with the cyclophilin inhibitor NIM811 treatment [108]. In another study, the zebrafish model was used to test mitochondrial

respiratory capacity after treatment with stable analogues of mPTP inhibitors [109]. Additionally, the zebrafish model has also provided insight into functional aspects of disease pathogenesis for several muscle conditions: for example, studies in zebrafish relatively relaxed (ryr) mutant, a model of RYR1-related myopathies [110], have contributed to identifying oxidative stress as an important disease mechanism in RYR1-related myopathies [111].

#### *3.4. Caenorhabditis Elegans*

With 40% of human disease genes having a nematode ortholog [112], and a fully sequenced genome [113], C. elegans is a valuable model to investigate several human physiological and pathological mechanisms. Studies of sarcomere maintenance and function in striated muscle led to the first identification of many conserved proteins, including twitchin, unc-89 (obscurin), unc-112 (kindlin), unc-45 (myosin chaperone) and unc-78 (AIP1) [114]. Using a large-scale screens in a C. elegans model of muscular dystrophy, carrying mutations in the dys-1 and the hlh-1 genes, which are respectively the homolog for the mammalian dystrophin and *MyoD* gene [115], compounds such as prednisone and serotonin have been shown to be effective in reducing muscle degeneration [116,117]. The obvious advantages of using this scalable and high-throughput model are counterbalanced by the limited phenotypic analyses, such as counting the number of times a worm bends in a C-shaped fashion in liquid in one minute, although new automated methods of quantifying muscle contraction and relaxation kinetics are emerging [118].

#### **4. Computational Models**

In silico models are becoming an increasingly useful tool for investigating muscle function and in helping us to understand which key players cause muscle pathology. These models integrate published experimental data, thus allowing us to encompass the many variables linked to pathology in a single model, enabling the study of multifaceted diseases. In doing these studies, one may understand better the underlying interactions between different disease mechanisms that lead to pathology, which may prove harder to do in live experiments. Over the last twenty years, big steps have been made in the computational modelling of muscle. A recent development has been the creation of agent-based models (ABMs), which allow us to assess what roles different biological agents play in muscle pathology, both at cellular and systems levels. For example, the use of ABMs for DMD has indicated a link between low satellite stem cell counts and impaired muscle regeneration symptom [119]. ABMs can also be used to predict the outcomes of given scenarios based on the rules derived from the literature, as well as having certain parameters that cannot be measured experimentally. This system can even add software agents that mimic certain biological cells into the simulation, with the aim of helping us to better understand their cellular interactions. This has been carried out in studies showing that fibroblasts can affect a muscle's susceptibility to disuse-induced atrophy [120].

However, these models do have their limitations: the simulated model is not a full replicate of the muscle cell and its microenvironment, as it only accounts for the contribution of known variables, which renders this model system not fully translatable to the in vivo situation.

#### **5. Conclusions**

The increasing availability of genetic and phenotypic information on patients with neuromuscular diseases, coupled with the unprecedented opportunity to manipulate eukaryotic genomes to generate disease models to study these diseases, has the potential to accelerate the translation of new therapeutic opportunities from preclinical settings into medical practice. Among the models available to researchers, 3D cultures and muscle on chips are best suited for precision medicine applications, due to their structural complexity and opportunity for genetic and environmental manipulation. However, as it becomes increasingly evident that we need to abandon the concept of 'one drug fits all', modelling every disease-associated variant for preclinical applications is likely to be unattainable and in many cases unnecessary. Achieving model precision is critical in translational research as long as it provides

predictive validity, which is the ultimate goal of preclinical work, and may further be enhanced by using multiple models to capture the spectrum of mechanisms and testing therapies in diverse genetic backgrounds that more closely reflect the human population as a whole. This may be particularly true in complex diseases, where multiple risk loci concur to the development of a specific condition or to the treatment response.

**Author Contributions:** Literature review, writing and original draft preparation, A.A.S.; writing, review and editing, A.A.S., R.E., T.G., and C.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **An Omics View of Emery–Dreifuss Muscular Dystrophy**

#### **Nicolas Vignier and Antoine Muchir \***

INSERM, Center of Research in Myology, Institute of Myology, Sorbonne University, 75013 Paris, France; n.vignier@institut-myologie.org

**\*** Correspondence: a.muchir@institut-myologie.org; Tel.: +33-1-4216-5705

Received: 13 May 2020; Accepted: 12 June 2020; Published: 15 June 2020

**Abstract:** Recent progress in Omics technologies has started to empower personalized healthcare development at a thorough biomolecular level. Omics have subsidized medical breakthroughs that have started to enter clinical proceedings. The use of this scientific know-how has surfaced as a way to provide a more far-reaching view of the biological mechanisms behind diseases. This review will focus on the discoveries made using Omics and the utility of these approaches for Emery–Dreifuss muscular dystrophy.

**Keywords:** *LMNA*; Emery–Dreifuss muscular dystrophy; Omics

#### **1. Introduction**

To understand the complexity of systems biology, Omics' technologies adopt a holistic view. In this, in opposition to hypothesis-generating experiments, no rationale is known, but instead, biological inputs are acquired and analyzed to delineate a hypothesis that can be then tested. Omics technology can be used not only to decipher physiological conditions but also in disease states, where they have a key role in diagnosis, as well as promoting our knowledge of the development of diseases [1]. Omics approaches to conditions, such as muscular dystrophies, are being used for small molecule therapy discovery by isolating innovative targets for drug development [2]. The scope of this review is to provide an overview of the Omics approaches and their application in Emery–Dreifuss muscular dystrophy research.

#### **2. Emery–Dreifuss Muscular Dystrophy**

Muscular dystrophies are characterized by the progressive weakness and degeneration of the skeletal muscle system, which may or may not be associated with cardiac impairment, leading to loss of mobility, and swallowing and respiratory difficulties. Death originates from respiratory defects or heart failure. Muscular dystrophies are a heterogeneous group of inherited disorders, and they differ in the distribution of affected muscles, the rate of muscle weakness progression and the age of onset [3]. The development of molecular genetic mapping techniques has shown that these disorders are genetically heterogeneous, and more than 50 genes have been identified as causing muscular dystrophies [4].

In the 1960s, an X-linked muscular dystrophy associated with contractures, which was first diagnosed as a benign variant of Duchenne muscular dystrophy, was reported [5,6]. In the 1980s, Alan Emery re-investigated the original family and reported that cardiomyopathy was a significant feature of the disease, which was thereafter called Emery–Dreifuss muscular dystrophy (EDMD) [7]. In EDMD, the onset of symptoms occurs within the first decade of life [7]. Contractures of the elbows, neck extensor muscles and Achilles' tendons appear to be the first symptoms of the disease, and occur before muscle weakness and wasting. The progressive muscle degeneration begins during the end of the second decade of life, in a humeroperoneal distribution. Cardiac alteration begins during the teenage

years, with no link to the severity of the muscular dystrophy [7–10]. Over time, dilated cardiomyopathy develops, and is associated with severe ventricular tachydysrrhythmias. Sudden cardiac death is frequent, and an implantable defibrillator can be lifesaving [10,11].

#### **3. Application of Omics Approaches for Emery-Dreifuss Muscular Dystrophy**

#### *3.1. Omics and Diagnosis*

In the 1990s, a positional cloning study—a technique for the positioning of a trait-associated gene within the genome—showed that mutations in *EMD*, and *LMNA* cause the X-linked [12] and autosomal dominant [13] forms of EDMD, respectively [14]. EMD encodes emerin, which is a transmembrane protein of the nuclear envelope. *LMNA* encodes nuclear lamins A and C, which are intermediate filament proteins associated with the nuclear envelope. Until recently, genetic screening for EDMD was performed with Sanger sequencing of the exons and intron-exon regions of both the *EMD* and *LMNA* genes. Since then, diagnosis of EDMD has been made available if the clinical signs are suggestive or if a family member is known to have EDMD. However, classical DNA sequencing methods have many shortcomings (time and cost) hampering their use for diagnosis. By contrast, as introduced in 2005, next-generation sequencing is a potent novel technology and a cost-effective method that has completely revolutionized the field [15,16]. Full exome sequencing is routinely used in clinical diagnostic laboratories to identify pathogenic variants in a given patient at a reasonable cost [17–20]. A new next-generation sequencing approach to identify potential new candidate EDMD genes has also recently been tested [21].

Subsequently, *LMNA* mutations have been shown to cause other striated muscle diseases, i.e., dilated cardiomyopathy [22], limb-girdle muscular dystrophy type 1B [23] and congenital muscular dystrophy [24]. Limb-girdle muscular dystrophy type 1B and dilated cardiomyopathy can occur in the same families as subjects with EDMD, and can be therefore be considered variants of the same clinical entity. This supports the concept that modifier genes arouse the severity and the peculiar symptoms. To map the modifier locus, microsatellite markers were genotyped in a large French family, where patients carrying the same *LMNA* mutation exhibited phenotypic variability [25]. The linked DNA region harbors two candidate modifier genes, *DES* and *MYL1*, encoding desmin and light chain of myosin, respectively, thus providing insights for the natural history and the physiopathology of EDMD.

#### *3.2. Omics and Abnormal Cellular Signaling*

The mechanisms by which mutations in *EMD* and *LMNA* cause muscular dystrophies are poorly understood. A few models have been proposed to explain the physiopathology of EDMD [26]. The model called the 'mechanical stress hypothesis' relies on the premise that striated muscle is steadily subjected to mechanical strains. Abnormalities in nuclear envelope composition may imply a weakening of the nucleus, which could represent an initial step in the chain of events leading to EDMD.

The ambition of transcriptomics studies is to pinpoint genome-wide changes and to expose coordinately organized gene networks [27,28]. Gene expression profiles associated with EDMD have been studied in several experimental models using various technologies. Mouse models have been helpful in deciphering mechanisms involved in the pathogenesis of EDMD, as well as for opening pharmacological therapies perspectives. The development of *Lmna*-/- mice by Sullivan and colleagues was the first animal model of the disease [29]. The mice expressing a truncated peptide, lamin A delta8-11 [30], develop cardiomyopathy and skeletal muscle wasting reminiscent of human EDMD. Then, other mouse models of EDMD were generated. These are knock-in mice that express A-type lamins with the p.H222P [31], p.deltaK32 [32] and p.N195K [33] residue substitutions. Arimura et al. developed *Lmna* knock-in mice carrying the p.H222P mutation that was identified in the human *LMNA* gene in an EDMD family [31]. This mutation was also chosen because it putatively dramatically altered the coil-coiled organization of A-type lamins, based on in silico analysis. This was the first *Lmna* mouse

model mimicking human EDMD from the gene mutation to the clinical symptoms. Two separate groups have generated *Emd* null mice [34,35]. These mice have subtle motor coordination abnormalities, with a prolongation of atrioventricular conduction time in *Emd* null mice [35]. Notwithstanding, these animal models may not reflect the "natural" human condition in terms of physiological mechanism and genetic outlook. Induced pluripotent stem cell (iPSC) technology represents a means to surmount these shortcomings, allowing the generation of any cell type through peculiar differentiation protocols [36]. Tissue-specific in vitro models of EDMD have been created from iPSCs that recapitulate traits of the disease [37–42].

To explore the pathogenesis of cardiomyopathy associated with EDMD, we carried out a genome-wide RNA expression analysis of hearts from *Lmna*p.H222P/H222P mice and *Emd* knockout mice, two mouse models of EDMD. This analysis revealed changes in the expression of genes encoding proteins in mitogen-activated protein (MAP) kinases, Wnt/β-catenin, AKT/mTOR and transforming growth factor (Tgf)-β signaling pathways in the mutated models [43–48]. Using RNA-sequencing technology on *Lmna*-/- mice, Auguste and colleagues showed that the FOXO signaling pathway impacted different signaling pathways, i.e., NFκB, TNFα, P53 and OxPHOS signaling pathways and biological processes, i.e., apoptosis, sustaining the cardiac phenotype associated with EDMD [49]. Furthermore, whole genome expression analysis of the primary cells of EDMD patients showed aberrant activity of unfolded protein response signaling [50]. In a cardiac specific expression of *Lmna* p.D300N, the Marian group showed using bulk RNA sequencing strategy that an increase of DNA damage response/TP53 pathway was contributing to the pathogenesis of cardiomyopathy associated with EDMD [51]. All these datasets led to the hypothesis of a model of how abnormalities of A-type lamins and emerin may lead to EDMD [52]. Using a transcriptomic approach from regenerating skeletal muscle from emerin-deficient mice, Melcon et al. [34] have shown delayed myogenic differentiation, which is regulated by *Rb* and *MyoD* genes. Since then, small molecules have been used by others to rescue the impaired myogenic differentiation in emerin deficiency, which could represent a potential strategy for improving the muscle wasting phenotype seen in EDMD [53]. Hence, abnormalities in satellite cell behavior may be responsible in part for the skeletal muscle disease in EDMD.

The mechanisms that bridge the *LMNA* genetic defects to malignant arrhythmias [54] are unknown. To better understand this phenotype in EDMD, Dr. Wu's group modeled the disease in vitro using patient-specific iPSC-derived cardiomyocytes (iPSC-CMs) carrying *LMNA* frameshift p.K117fs mutation [42]. They showed an abnormal activation of the platelet-derived growth factor (PDGF) signaling in *LMNA* p.K117fs iPSC-CMs. The inhibition of the PDGF signaling improves the arrhythmic phenotype of mutated iPSC-CMs, opening novel therapeutic perspectives for the treatment of EDMD.

All these studies contributed to a better knowledge of the functional and molecular mechanisms of the disease. We can expect that new findings will design applications of iPSC-models to pharmacological testing in striated muscle-specific contexts [55], making the technology available to patients.

#### *3.3. Omics and Chromatin Regulation*

Among the models proposed to explain how EDMD phenotypes arise, the 'gene-expression model', posits an effect of mutated lamin A/C on the transcription activity of genes and/or pathways that could impact striated muscle-homeostasis process. According to this model, it has been described that A-type lamins interact with heterochromatin regions called Lamin-Associated Domains (LADs). These LADs play a role for chromatin organization and gene expression regulation [56–59]. It has been shown that the LADs are re-organized in EDMD steering modifications of the epigenetic program, ultimately driving the loss of myogenic differentiation [60]. This has been recently shown in an elegant manner by Bianchi and colleagues on a murine model of EDMD [61]. The authors described an abnormal positioning of polycomb proteins, which are epigenetic repressors involved in cell identity. This causes impairment in self-renewal, deficiency of cell identity and the early exhaustion of the quiescent satellite cell pool, demonstrating that muscular dystrophy in EDMD can be partially caused by epigenetic dysfunctions of muscle stem cells [61]. Mewborn and colleagues showed that

the *LMNA* p.E161K mutation perturbed the positioning and compaction of chromosomal domains in primary fibroblasts, resulting in an altered gene expression profile [62]. Another study focused on the organization of LADs in EDMD using a multi-Omics approach (Chip-Seq/RNA-sequencing) in explanted hearts from five patients carrying *LMNA* mutations [62]. LADs were redistributed, ensuing in a functional chromatin state in mutated hearts, suggesting the loss of specific functional chromatin binding. This aberrant distribution impacted both the gene expression profile and CpG methylation. An integrated analysis showed the combined role of LADs and CpG methylation in the regulation of gene expression, and identified numerous transcription factors involved in biological processes such as cell death/survival, cell cycle and metabolism [63]. Bertero and colleagues showed at the same time, using genome-wide chromosome conformation capture (HiC) analyses on iPSC-CMs carrying *LMNA* p.R225X mutation [64], a slight chromatin compartment dysregulation (around 1% of the genome). RNA-sequencing of these altered chromatin domains revealed the abnormal up-regulation of only a handful of genes, among which was *CACNA1A*. This latter encodes a subunit of a calcium channel, whose abnormal expression might partially explain the electrophysiological and contractile aberrations observed in the mutant hiPSC-CMs. The authors showed that pharmacological treatments to prevent both the electrophysiological and contractile alterations helped improve the abnormal phenotypes described in the mutated iPSC-CMs. This suggests that *CACNA1A* may be a good therapeutic target to reverse the cardiac abnormalities in EDMD patients. However, the work by Bertero et al. showed only minor alterations in chromatin compartmentalization in mutated hiPSC-CMs, a finding that challenges the aberrant gene-expression model [64].

#### *3.4. Omics and Metabolism*

Contracting incessantly, the heart demands a lot of energy to ensure optimal contractile function. Research has demonstrated that the high requirements of the heart are satisfied by a preference for the oxidation of fatty acids. Studies have demonstrated that the failing heart deviates from its inherent profile and relies heavily on glucose metabolism, primarily achieved by acceleration in glycolysis. To gain further insight into the molecular mechanism ruling this disease, scientists have studied the cardiac metabolic rates in *Lmna*-/- mice. West and colleagues described a targeted metabolomics assay that quantifies metabolites relevant to cardiac metabolism [65]. The assay demonstrates that the *Lmna*-/- mouse heart has decreased metabolites associated with the citric acid cycle and fatty acid oxidation [65]. This corroborated another study, which showed that activated AKT/mTOR signaling reduces tolerance to energy deficits in hearts from *Lmna*p.H222P/H222P mice [45]. A rapamycin analog that blocks AKT/mTOR activity has been used to prevent the progression of cardiomyopathy in *Lmna*p.H222P/H222P mice [45]. These works highlighted that the heart is unable to compensate for increased or fluctuating energy demand and, over time, develops dilated cardiomyopathy in EDMD. This metabolic remodeling probably represents an adaptive cardio-protective mechanism that can help improve contractile function, thus slowing the progression of EDMD and improving prognosis. As such, metabolic modulators, which have the potential to shift myocardial substrate utilization from fatty acids toward glucose metabolism, may have a place in the management of patients. Some of these modulators have already been investigated as treatment for cardiomyopathy, with some beneficial effects [66,67]. It would be relevant to test their efficacy in EDMD.

Moreover, West and colleagues showed increased responses to oxidative stress and reactive oxygen species (ROS) exposure in *Lmna*-/- mouse hearts [65]. ROS are small, short-lived signaling molecules that mediate various cellular responses. Based on this, we recently showed that N-acetyl cysteine treatment reduces cardiac oxidative stress injury and ameliorates contractile dysfunction in *Lmna*p.H222P/H222P mice [68].

#### *3.5. Omics and Biomarkers*

Omics are powerful tools to identify diseases' molecular biomarkers. Molecular biomarkers are molecules with particular biophysical properties, the quantities of which are measured in biological samples, and are of critical importance in the support of the clinical diagnosis of pathology, the monitoring of its time course and the evaluation of the impact of therapeutic approaches. Molecular biomarkers have to be quantified from patients in the least invasive manner possible. Circulating molecular biomarkers in body fluids, i.e., blood, plasma, serum or urine, are thus the main interest. Proteomic, transcriptomic and metabolomic analysis have been driven in many diseases to identify such biomarkers [69,70].

To identify circulating microRNA as molecular biomarkers for EDMD, the microRNA transcriptome (miRnome) from plasma of *Lmna*p.H222P/H222P mice was screened [71]. A specific and distinctive microRNA expression profile was identified, and three of the dystromirs (mir-133b, mir-133a and mir-1) were downregulated in these mice. Furthermore, two microRNAs were upregulated (mir-146b and mir-200a) and six microRNAs were downregulated (mir-130a, mir-133a, mir-133b, mir-1, mir-151-3p, mir-339-3p) in *Lmna*p.H222P/H222P mice compared with wild type animals [71].

#### **4. Conclusions**

Omics technologies have hastened the identification of genetic mutations associated with EDMD, and have unveiled the existence of rare variants and modifier genes that might establish the phenotypical heterogeneity of the disease. Transcriptomic studies have uncovered alterations in signaling mechanisms causing some of the symptoms in EDMD, but are restricted to experimental models (in vitro and in vivo) with technical and theoretical shortcomings. Furthermore, although the number of parameters being measured has increased with Omics technologies, the number of biological and methodological replicates has not. In addition, because of the large number of measurements and the limited number of subjects, unique problems arise in Omics studies involving statistics and bias. A single Omics technique will only capture changes in a subset of biological cascades; it cannot provide a systemic understanding of the complexity of systems biology. The integration of multiple Omics data sets promises a substantial improvement, through an increase of information and, especially, systemic understanding. Therefore, much work is needed before using this research in the clinic. All of this will depend on integrative collaborations among physicians and scientists that will be essential for major breakthroughs for both the diagnosis and treatment of EDMD.

**Author Contributions:** Writing—Original Version, A.M.; Writing—Review & Editing, A.M. and N.V.; Funding Acquisition, A.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Association Française contre les Myopathies, the Institut National de la Santé et de la Recherche Médicale and Sorbonne Université.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Review*
