Skip to main content
PMC home page
European Journal of Human Genetics logoLink to European Journal of Human Genetics
. 2023 Mar 20;31(6):663–673. doi: 10.1038/s41431-023-01329-5

A genome-wide association analysis of loss of ambulation in dystrophinopathy patients suggests multiple candidate modifiers of disease severity

Kevin M Flanigan 1,2,3,✉, Megan A Waldrop 1,2,3, Paul T Martin 1,2, Roxane Alles 1, Diane M Dunn 4, Lindsay N Alfano 1,2, Tabatha R Simmons 1, Melissa Moore-Clingenpeel 5,6, John Burian 5, Sang-Cheol Seok 5, Robert B Weiss 4,#, Veronica J Vieland 2,5,6,7,#
PMCID: PMC10250491  PMID: 36935420

Abstract

The major determinant of disease severity in Duchenne muscular dystrophy (DMD) or milder Becker muscular dystrophy (BMD) is whether the dystrophin gene (DMD) mutation truncates the mRNA reading frame or allows expression of a partially functional protein. However, even in the complete absence of dystrophin, variability in disease severity is observed, and candidate gene studies have implicated several genes as modifiers. Here we present the largest genome-wide search to date for loci influencing severity in N = 419 DMD patients. Availability of subjects for such studies is quite limited, leading to modest sample sizes, which present a challenge for GWAS design. We have therefore taken special steps to minimize heterogeneity within our dataset at the DMD locus itself, taking a novel approach to mutation classification to effectively exclude the possibility of residual dystrophin expression, and utilized statistical methods that are well adapted to smaller sample sizes, including the use of a novel linear regression-like residual for time to ambulatory loss and the application of evidential statistics for the GWAS approach. Finally, we applied an unbiased in silico pipeline, utilizing functional genomic datasets to explore the potential impact of the best supported SNPs. In all, we obtained eight SNPs (out of 1,385,356 total) with posterior probability of trait-marker association (PPLD) ≥ 0.4, representing six distinct loci. Our analysis prioritized likely non-coding SNP regulatory effects on six genes (ETAA1, PARD6G, GALNTL6, MAN1A1, ADAMTS19, and NCALD), each with plausibility as a DMD modifier. These results support both recurrent and potentially new pathways for intervention in the dystrophinopathies.

Subject terms: Genome-wide association studies, Neuromuscular disease

Introduction

Mutations in DMD, encoding the dystrophin protein, result in the X-linked dystrophinopathies. The most common of these are the severe Duchenne muscular dystrophy (DMD) and the milder Becker muscular dystrophy (BMD), with some patients demonstrating an intermediate phenotype. The primary determinant of disease severity depends upon whether the mutation truncates the DMD open reading frame, ablating dystrophin expression, or maintains an open reading frame, allowing expression of a partially functional muscle-specific (Dp427m) dystrophin isoform. Mutation analysis is often accompanied by interpretations of the predicted reading frame, which is used by clinicians to predict phenotype. Exceptions to this “reading frame rule” as interpreted from genomic DNA frequently occur, with most of them explained by molecular mechanisms that restore an open reading frame in the mRNA [1].

Although treatment with corticosteroids may alter the age of loss of ambulation (LOA) in patients who lack dystrophin expression, increasing evidence implicates other genetic modifiers of severity, including proteins involved in inflammation, fibrosis, regeneration, and sarcolemma repair [2–4]. Protein polymorphisms in latent TGF-β binding protein 4 (Ltbp4) are associated with increased levels of TGFβ-mediated fibrosis in mice [5], and a common human LTBP4 haplotype has been associated with disease severity [6]. The matricellular protein osteopontin (Spp1) contributes to the balance between inflammatory, regenerative, and fibrotic processes in muscle, and common regulatory SNPs in the proximal promoter of SPP1 are associated with variations in age at LOA in DMD patients [7]. Cohort effects of SPP1 and LTBP4 SNPs are evident with variable success at replication [3], but the effects of the “protective” and “risk” human LTBP4 polymorphisms on the dystrophic phenotype in mice have been confirmed, with epistatic interactions between Spp1 and Ltbp4 in modulating fibrosis via pathological TGFβ signaling [8]. The sarcolemmal repair protein and downstream effector of TGFβ signaling annexin A6 (Anxa6) modulates severity in mice, with alleles of Ltbp4 and Anxa6 acting jointly to modify the dystrophic phenotype [8]. The Notch signaling pathway gene Jagged1 modifies severity in the golden retriever muscular dystrophy (GRMD) dog model, implicating a role for satellite cell self-renewal and quiescence in severity [9].

As an extension of candidate modifier gene studies in DMD patients, there have been three previous reports of GWAS for age at loss of ambulation. The first was based on an exome chip (27 K markers) and 109 ancestrally homogeneous samples [2]. While no loci reached traditional genome-wide thresholds for statistical significance, subsequent filtering for genes in candidate pathways identified a haplotype that included a known functional SNP, rs1883832, in the Kozak sequence of the CD40 [2, 7] gene, which encodes a costimulatory receptor mediating immune and inflammatory responses in a variety of cell types. More recently, we reported a preliminary GWAS using 253 non-ambulant subjects [4] and observed two loci above the genome-wide significance threshold using the recessive model suggested by the prior candidate gene study. One locus was LTBP4 itself, where fine mapping revealed cis-eQTL variants in linkage disequilibrium with “protective” IAAM LTBP4 protein isoform and associated with LTBP4 mRNA expression. The other locus was a distal enhancer containing cis-eQTL variants associated with thrombospondin-1 (THBS1) expression, a matricellular protein with multifunctional properties including promoting TGFβ signaling through latent protein complex activation [10]. Finally, a recent two-stage design using whole exome sequencing of extreme phenotypes identified TCTEX1D1 as an LOA modifier in cohorts of 301 and 109 DMD patients [11]. This study highlighted the importance of novel study designs in searching for modifiers of DMD, but also emphasized the need for more complete knowledge of mutation-specific effects that may increase residual dystrophin levels.

The current work is distinct from our previous GWAS study in several critical ways. First, we have substantially increased the cohort to N = 419 dystrophinopathy patients. We originally published a candidate gene study in 253 non-ambulant patients [6]; we then published a genome scan in those same 253 subjects [4]. Here, using the same genotyping platform, we include all 174 of the original cohort who met current inclusion criteria for this study (see below), as well as 245 newly collected subjects, while utilizing survival analysis to permit inclusion of subjects who were still ambulatory at last assessment. Second, we have used a more conservative approach to excluding subjects whose DMD mutations may produce residual dystrophin, to minimize confounding of modifier gene effects with LOA variation due to mutations within the DMD gene itself, as described in detail below. This step is designed to increase the power of the analysis by minimizing heterogeneity in the causes of LOA variation within our data set. (Of note, 54 subjects from the original N = 253 cohort were dropped from the current cohort based on this new approach to DMD mutation classification.) Third, instead of the usual (frequentist) approach to assessing evidence for genetic association, we have utilized an evidential statistical paradigm that may be better adapted to analysis in smaller data sets than conventional approaches [12]. This is critical for the study of a rare disease like DMD, for which the accrual of data sets on the order of typical GWAS designs (thousands or even tens of thousands of subjects) is not feasible. Finally, we uniformly apply recent genome-wide datasets to the candidate SNP data to reveal their functional impact, utilizing a systematic in silico pipeline integrating ENCODE genomic annotations, chromatin interaction experiments, expression quantitative trait loci (eQTL), and in silico functional predictions, allowing us to implicate several novel regulatory loci as potential modifiers of DMD severity and to consider the biological plausibility of these candidate modifiers as a first step toward further studies to determine pathogenesis.

Materials and methods

Clinical classification

All subjects have a known DMD mutation confirmed by either MLPA or DNA sequencing, as previously described [13, 14]. The study was conducted under research protocols approved by the Nationwide Children’s Hospital Institutional Review Board (Approval #0502HSE046) and the University of Utah Institutional Review Board (Approval #IRB_00094017). Written informed consent was obtained from the parent/guardian of each minor subject, or from the subject themselves if 18 years of age or older. Genomic DNA was extracted from blood samples and phenotypic data were extracted from clinical records; most subjects underwent prospective clinical examinations as well.

Ambulatory loss was defined by full-time wheelchair use (requiring wheelchair within the home) and LOA was recorded to the nearest month when available, and otherwise to the nearest half-year of age. Glucocorticoid treatment (Steroid Y) was defined as use of any steroid regimen for more than 6 months that began at least 6 months prior to loss of ambulation. Untreated (Steroid N) was defined as never exposed, exposed to any steroid regimen for less than 6 months, or onset of treatment less than 6 months prior to loss of ambulation. Data was first entered in the UDP in the first decade of the 2000s, when steroid use was less widespread, reflected in the relatively low number of subjects considered to have been treated with steroids. However, all genotyped subjects had complete clinical data, including information on steroid use. Only one subject was on a potentially disease-modifying drug (eteplirsen) for more than 6 months prior to loss of ambulation.

Inclusion criteria for subjects

In previous studies, we used physician-assigned diagnoses of DMD and BMD, which are themselves defined clinically by age at LOA, and therefore excluded individuals who still ambulated at age 20. Here we include patients regardless of their age at LOA. However, we excluded any individuals still ambulant with last assessment at ≤6 years of age (the minimum observed LOA in our cohort), as these individuals are uninformative for the association analysis.

The primary criterion for inclusion was a conservative definition of loss-of-function (LoF) DMD mutations. The impact of trace dystrophin protein expression [15] was minimized by classifying DMD mutations using a conservative LoF definition from predicted muscle-isoform transcription, splicing, and translation of their DMD mutation (S1 and S2 Table). We excluded typical Becker mutations producing moderate to high levels of partially functional dystrophin, for example, “in-frame” deletions in the central rod domain, nonsense mutations in exon 1 that use a downstream translation initiation site [16], and “out-of-frame” exon 3 to 7 deletions [17]. However, we also excluded patients with mutations that typically have been reported as “out-of-frame” or nonsense mutations on clinical reports but may occasionally result in low-level dystrophin expression based upon a known or suspected mechanism for bypassing frame-truncating DMD mutations. Excluded mutations included “out-of-frame” exon 45 deletions, which sometimes can provoke low-level exon 44 skipping that restores the reading frame, as well as exon 8 “skippable” mutations [18]. We and others have previously demonstrated that exons 23–42 constitute a “zero-frame context,” in which deletion of any exon(s) results in an open reading frame. Nonsense mutations and splice site mutations within these exons can be associated with exon skipping due to poorly-predicted damage to exon definition elements recognized by the spliceosome that restores the open reading frame by excluding the mutated exon from the mature mRNA [1, 19, 20]. Because of this strict LoF mutation criteria, we included 311 non-ambulant LoF subjects, but excluded 166 non-ambulant partial LoF subjects, of which 132 had lost ambulation by 20 years (Fig. 1A).

Fig. 1. DMD mutation class and LOA cohort selection.

Fig. 1

Open in a new tab

A Each panel contains UDP cohort patients in a dystrophin mutation class, with the upper left panel showing loss-of-function (LoF) subjects included in the GWAS (N = 419, yellow highlight), and the remaining seven panels showing excluded partial LoF subjects (N = 414, gray highlight; see Tables S1B and S2B). The y-axis is age at last visit for ambulant patients or age at loss of ambulation. Non-ambulant patients are grouped by new subjects (NonAmbA) or subjects from a previous GWAS (NonAmbB) that used loss of ambulation by 20 years of age as an inclusion criterion (horizontal blue line). Individual patients are colored by their clinical diagnosis: DMD (red), IMD (orange), and BMD (green). B Flow diagram for cohort inclusion criteria, phenotyping, GWAS survival analysis, and candidate gene identification. C Kaplan-Meier curves and violin plot of the OTE residuals in the censored LoF cohort (N = 419 patients). The OTE residuals on the y-axis are deviations in age at LOA from expectation, after adjusting for steroid use, on the standard normal scale. Smaller (negative) values represent older-than-expected LOA, while larger (positive) values represent younger-than-expected LOA.

Genotyping methods and cleaning, methods of statistical analysis (based on the posterior probability of (trait-marker) linkage disequilibrium (PPLD), computed using the software package KELVIN [21]), and methods of functional annotation and eQTL analysis of candidate loci are all discussed in detail in Supplemental Methods.

Results

Description of the cohort

Using the LoF classifications described above, our cohort of N = 419 subjects comprised 201 distinct LoF DMD mutations, encompassing 71 multi-exon deletions, 52 nonsense, 23 multi-exon duplications, 21 frameshifts, 14 single exon deletions, 11 single exon duplications, 6 splice donor, and 3 splice acceptor mutations (see Tables S1A and S2A), while N = 414 subjects with partial LoF mutations were excluded (Fig. 1A). Recent observations of patients with exon 2 duplications cataloged within the UDP database suggests that many of these patients may have delay in LOA in relation to other mutations, perhaps due to utilization of a downstream internal ribosome entry site to initiate translation of a highly functional dystrophin; Dup2 patients were not excluded from the LOF group, but are discussed below.

Kaplan-Meier estimates of the survival function S (time to LOA) based on the N = 419 LoF cohort are shown in Fig. 1B, separately for the two steroid treatment groups. The estimated mean (standard deviation) of LOA is 10.1 (2.6) and 12.3 (3.4) for the steroid N and Y groups, respectively. As can be seen, there is considerable variability in LOA, with no individuals losing ambulation prior to age 6, but a substantial number (44) still walking beyond age 15. The distribution of the Ordinary Time-to-Event (OTE) residuals [22] as the phenotypic values used as input for PPLD analysis are also shown in Fig. 1C. The proportion of subjects still ambulatory (censored) as of last examination was 26.7. 69.8% of subjects had been treated with steroids, and 30.2% were untreated.

GWAS for loss of ambulation (LOA)

Figure 2A shows the primary genome-wide scan results of the PPLD analysis. 94.4% of SNPs returned a PPLD < 0.0004, indicating evidence against association with LOA. There were 74 (0.0053%) SNPs with PPLD ≥ 0.04, 28 (0.0020%) with PPLD ≥ 0.10, and 8 (0.0006%) with PPLD ≥ 0.40; this last group represented 6 distinct loci. (For reference, in 1,000,000 independent replicates of size N = 400 simulated under the null hypothesis of no association, we observed 15, 7 and 1 PPLDs ≥ 0.04, 0.10 and 0.40, respectively.) Note that the relatively “sparse” appearance of the PPLD Manhattan plot is largely a scaling issue. The PPLD drops off rapidly with small changes in SNP-to-SNP LD. Combined with the very small prior probability of trait-marker association, this leads to a very different visual impression than what is seen on typical p-value plots, as illustrated in the LocusZoom plots of the rs12657665 region on chr5 (Fig. 2B). A PPLD of 0.04 is still 100 times greater than the prior probability, hence, still supporting association, but on the Manhattan plot it will appear visually indistinct from 0. The PPLD signals reported here are consistently and substantially above the prior probability of association at those SNPs in high LD with the peak SNP in each region. Despite the unusual appearance of the Manhattan plot in this regard, the PPLD produces fewer false positive signals than conventional regression analysis [12].

Fig. 2. Manhattan plot of genome-wide scan.

Fig. 2

Open in a new tab

A The PPLD is on the probability scale and represents the estimated posterior probability that a SNP affects the phenotype, either directly or via linkage disequilibrium with one or more causal variants. Horizontal lines correspond to the heuristic cutoffs used to prioritize SNPs for further investigation in what follows. B LocusZoom plots of the lead SNP (rs10077875) from chr5 using PPLD values (upper panel) or p-values from a Cox Proportional Hazards regression with steroid exposure and SNP genotypes as covariates (lower panel).

Table 1 shows the list of 19 distinct loci with PPLD ≥ 0.10, indicating the peak genotyped SNP as well as the maximum PPLD value in each region including imputed genotypes (Table 1). While imputation increased the maximum PPLD at some loci, it did not change the set of loci with PPLD ≥ 0.40; it did however somewhat shuffle the rank-ordering within this set. For comparison, Table 1 also shows the p value and rank for these top SNPs from a Cox proportional-hazards regression (CPH) using the OTE residuals and an additive model. Although this is an “apples and oranges” contrast (again, see [12] for a detailed comparison of PPLD to CPH analysis), it does indicate that 5 out of 6 loci with PPLD ≥ 0.40 have some level of support from CPH regression. We also evaluated the robustness of these top PPLD values after removing 30 subjects with a stricter definition of ancestry (S1) or excluding 12 subjects with DMD Dup2 mutations, as discussed above. Large decreases in PPLD values were observed for three loci after further ancestry filtering and two loci after Dup2 filtering (Table 1). Four of the top six loci were unaffected by the additional filtering criteria, while two of the top six loci (chr8 rs72681143 and chr5 rs10077875) showed a moderate decrease in PPLD values after Dup2 exclusion; note that the resulting PPLDs were still considerably higher than the prior probability of 0.0004.

Table 1.

SNPs with PPLDs ≥ 0.10.

Chr Position Marker PPLD Max PPLD CPH Add p-value CPH rank order Ancestry Filtered PPLD (n = 389) Dup2 Filtered PPLD (n = 407) V2G top gene Gencode
2 67958890 rs34263553 0.87 0.87 1.9E-03 4352 0.36 0.79 ETAA1 intergenic
8 103189435 rs72681143 0.77 0.83 0.45 642282 0.80 0.05 NCALD intergenic
4 172432296 rs1358596 0.42 0.75 1.2E-03 3016 0.65 0.54 GALNTL6 intergenic
6 120491553 rs10499096 0.62 0.64 3.0E-04 880 0.68 0.73 MAN1A1 intergenic
5 128722696 rs10077875 0.44 0.62 6.1E-05 239 0.89 0.03 ADAMTS19 intergenic
18 77551738 rs2061566 0.40 0.40 1.0E-05 50 0.50 0.22 PQLC1 intergenic
2 171307693 rs1816670 0.28 0.38 3.6E-04 1017 0.27 0.0075 MYO3B intronic
11 97137756 rs78973513 0.23 0.23 6.9E-06 31 0.25 0.09 JRKL intergenic
10 13049885 rs7921121 0.17 0.17 0.01 20370 0.13 0.11 CCDC3 intronic
6 31205392 rs3130530 0.16 0.14 0.10 158686 0.11 0.18 HLA-C intergenic
15 59429111 rs11539756 0.16 0.16 0.03 45865 0.15 0.0176 MYO1E UTR3
10 134349181 rs11146413 0.14 0.08 3.7E-04 1048 0.0017 0.28 INPP5A intergenic
13 101558219 rs9513818 0.13 0.32 0.04 64644 0.23 0.0036 TMTC4 intergenic
5 58272070 rs62356653 0.12 0.12 1.3E-05 56 0.08 0.19 PDE4D intronic
10 1449849 rs4880848 0.11 0.11 0.90 1222924 0.0037 0.11 ADARB2 intronic
15 61272080 rs2140442 0.11 0.11 0.45 644035 0.10 0.2 RORA intronic
6 47227944 rs2281448 0.10 0.10 0.31 459744 0.0002 0.06 TNFRSF21 intronic
12 94954104 rs4761483 0.10 0.05 0.06 97480 0.11 0.0087 PLXNC1 intronic
14 88421482 rs113691242 0.10 0.07 6.5E-04 1646 0.11 0.0192 GALC intronic
Open in a new tab

Position: Physical positions are based on build 37. Marker: shown is the top genotyped SNP per region, where the regions are defined as in the main text. PPLD: The PPLD at the genotyped SNP. Max PPLD: The maximum PPLD in the region, which may occur at the originally genotyped SNP or at an imputed SNP. CPH Add_P: CPH Add_P: Cox Proportional Hazards Regression p value with steroid exposure and SNP genotypes as covariates. CPH rank order: additive model ranked by p-value from 1,353,208 directly genotyped SNPs. Ancestry Filtered PPLD: 30 subjects with second criteria applied for ancestry filtering. Dup2 Filtered PPLD: 12 subjects removed with Dup2 Dmd mutations. V2G top gene: Variant to Gene nearest protein coding gene. Gencode: gencode annotation category for the Max marker.

Regional analysis of lead SNPs

Using a PPLD threshold ≥ 0.05 to define credible sets from the top 6 regions, 96 total SNPs were chosen for further analysis. Regional LocusZoom plots for these six loci are shown in Figs. 3A, 4A, and S2, with Kaplan-Meier survival estimates and genotype-specific OTE residual plots for each lead SNP. Note that survival curves that cross one another indicate evidence for genotypic effects on variances. Functional annotations of individual SNPs (Table S3) revealed that all 96 SNPs were noncoding intergenic or intronic variants. Four of the credible sets displayed a very narrow range of pairwise SNP-to-SNP linkage disequilibrium (r2) of 0.96–0.98 for chr4 rs1358596 (Fig. S2D), 0.90–0.96 for chr2 rs34263553 (Fig. 3A), 0.89–0.98 for chr6 rs10499096 (Fig. S2A), and 0.84–0.98 for chr5 rs12657665 (Fig. S2B). The chr8 rs3899035 top SNP had a broader range of r2 values from 0.27 to 0.99 (Fig. S2C), while the chr18 rs2061566 SNP was the only variant in the credible set as no other SNPs in the study sample had an r2 value > 0.3 (Fig. 4A), similar to the rs2061566 LD pattern seen with European ancestry samples in the 1000 Genomes Project.

Fig. 3. Regional eQTL and chromatin features for the chr2 SNP region.

Fig. 3

Open in a new tab

A The left panel shows a LocusZoom plot of chr. 2 SNPs with their PPLD values on the primary y-axis and colored by their LD with the lead SNP rs342633553. The middle panel shows the estimated (Weibull) survival curves for the two Steroid (N/Y) groups, respectively; legends include the number n of individuals followed by [10th percentile, median, 90th percentile] of the age at LOA distribution. The right panel shows the density plots for the OTE residuals with the y-axis units in deviations in age at LOA from expectation, after adjusting for steroid use. Note that residuals are oriented such that smaller (negative) values represent older-than-expected LOA, while larger (positive) values represent younger-than-expected LOA. Red dots indicate the mean within each group. B The left panel includes PPLD values for genotyped and imputed SNPs, SNP-ETAA1 eQTL p-values from T cells (Table S2), and total RNA-Seq read depths from normal skeletal muscle tissue. The bottom track is a one-to-all Hi-C chromatin interaction plot from psoas muscle where the y-axis (blue density) indicates bias-removed interaction frequencies, and the horizontal green line is a cut-off value (=2) for the distance-normalized frequency. The right panel includes the credible SNP region (gray bar) with ENCODE cCREs annotations, vocabulary-annotated index DHSs elements, and biologically labeled component meta-DNase I tracks from the top 15 ENCODE biosamples most strongly associated with each component. C Overlap of the chr2 PPLD credible SNPs and lead SNP rs34263553 with the ETAA1-eQTL shown as both the nominal -log10 p-value of the eQTL association and the fine mapped SuSiE posterior inclusion probability (PIP) of SNPs in the eQTL credible set. Additional detail is shown for the ENCODE cCREs features, as well as conserved sequence elements from vertebrate alignments scored with PhyloP.

Fig. 4. Regional chromatin features and TF footprints in the chr18 rs2061566 region.

Fig. 4

Open in a new tab

A LocusZoom plot of SNPs from chr18:77,150,001–78,020,000 (hg19), with survival curves and OTE residual plots as described in Fig. 3A. B The left panel includes PPLD values, skeletal muscle total RNA-Seq read depths, RefSeq Select transcripts, and promoter capture Hi-C chromatin interactions from psoas muscle using the PARD6G chr18:77,995,738–78,007,496 promoter fragment. The right panel includes the credible SNP region (gray bar) with annotation tracks for ENCODE cCREs, vocabulary-annotated index DHSs elements, and biologically labeled component meta-DNase I tracks. C Finer details for ENCODE cCREs, vocabulary-annotated index DHSs elements, and biologically labeled component DNase I tracks as in (B). Additional tracks include DNase I cleavage pattern from ref. [29], showing DNaseI per-nucleotide cleavage and the digital TF footprints within the DHS at an FDR cutoff = 0.05. The location of rs2061566 is highlighted by the light-yellow bar.

Plausible candidate gene identification

Table 2 shows the lead candidate SNP with a PPLD score ≥ 0.4 from each of the associated regions assigned to a candidate protein coding gene using the Open Target Genetics V2G pipeline with the gene expression data (eQTL) supporting this connection. These six lead SNPs all tagged non-coding credible SNP sets (Table S3), therefore we examined these genomic regions for biological features indicative of functional regulatory variation. These included ENCODE candidate cis-regulatory elements (cCREs), ENCODE indexed DNase I hypersensitive sites (DHSs), ENCODE transcription factor footprints, gene proximity, enhancer-promoter chromatin interactions, skeletal muscle total RNA expression, and cis-eQTLs SNP-gene pairs from the GTEx and the eQTL Catalog databases. Detailed analysis of evidence for two plausible candidate gene (ETAA1 and PARD6G) follows, while additional analyses are contained in Supplemental Discussion for additional genes (GALNTL6, MAN1A1, ADAMTS19, and NCALD). Since mRNA transcripts for several modifier genes (SPP1 and THBS1) show large increases in dystrophic skeletal muscle, we examined differential expression of known and candidate modifiers using RNA-seq data from 1 month-old WT versus Dmd ∆Ex51 mice. As expected, Spp1 and Thbs1 transcripts showed significant increases in dystrophic skeletal muscle (Fig. S3A), while the other known modifiers and the top candidate modifiers, except for Fbn2 in the Adamts19 region, did not show significant differential expression (Fig. S3B).

Table 2.

Genes associated with LOA-associated SNPs.

Chr rsID Gencode Category PPLD MAF Overall V2G (distance to canonical TSS) eQTL SNP-to-gene Tissue eQTL p value
2 rs34263553 intergenic 0.87 0.13 ETAA1 (334,439 bp) ETAA1 T cells 1.70E-05
4 rs1358596 ncRNA_intronic 0.75 0.19 GALNTL6 (301,109 bp) GALNTL6 T cells 7.20E-04
5 rs12657665 intergenic 0.62 0.48 ADAMTS19 (76,055 bp) ADAMTS19 Sk. Muscle 9.70E-04
6 rs10499096 intergenic 0.64 0.05 MAN1A1 (820,646 bp) – – –
8 rs3899035 intergenic 0.83 0.11 NCALD (53,501 bp) NCALD Monocytes 4.30E-14
18 rs2061566 upstream;downstream 0.4 0.09 PQLC1 (159,919 bp) PQLC1 Neutrophils 2.70E-16
Open in a new tab

Chr2 rs34263553 and ETAA1

The nearest protein-coding gene to the chr2 rs34263553 SNP with the highest PPLD value (0.87) is the Ewing tumor-associated antigen 1 (ETAA1) gene, at a distance of 335 kb (Fig. 3A). Chromatin looping data support a physical interaction between the rs34263553 region and the ETAA1 promoter region in psoas muscle (Fig. 3B and Table S4). The SNP region also coincides with cis-eQTL data associated with ETAA1 expression in T cells (Fig. 3B and Table S5). The credible rs34263553 SNP interval contains multiple candidate cis-regulatory elements (cCREs, Fig. 3B) and the lead SNP rs34263553 is also located within a stromal B index DNase I hypersensitive site (DHS, Table S6), flanked upstream by a phylogenetically enhancer cluster and downstream by a CCCCTC-Binding factor site with a prominent DHS peak in multiple tissues/cells (Fig. 3B and C). The 12 SNPs in the rs34263553 PPLD credible set coincided completely with the 13 SNP fine mapped ETAA1 eQTL credible set defined by the SuSiE method and rs34263553 had the highest posterior probability of being the causal eQTL variant for ETAA1 expression (Fig. 3C). In the mouse, Etaa1-deficient T cells display elevated DNA damage during proliferation [23] and the negative beta value in T cells for the cis-eQTL suggests that the minor allele of rs34263553 is associated with decreased ETAA1 expression in humans. Further discussion of the plausibility of this candidate gene is contained in Supplementary Discussion.

Chr18 rs2061566 and PARD6G

The chr18 SNP rs2061566 A/G variant near the chr18q subtelomere region (Fig. 4A) is not in strong LD (r2 < 0.3) with other SNPs in this study or in the 1000Genomes European data set. The SNP resides within the cCRE enhancer EH38E1929501 that has prominent DHS peaks in stromal A and B components (Fig. 4B) enriched for ENCODE biosamples derived from connective tissues and cells, and in the myeloid/erythroid component (Fig. 4C, Table S6), which is enriched for biosamples derived from hematopoietic stem cells. Multiple, cell-type specific transcription factor (TF) footprints are evident within the DHS peak for E1929501 from both psoas muscle and CD34 + hematopoietic stem cells (Fig. 4C). The top pCHi-C interaction observed in psoas muscle (Fig. 4B, Table S4) connected the E1929501 enhancer to the par-6 family cell polarity regulator gamma (Partitioning Defective 6, PARD6G) gene. PARD6G is a member of the PAR polarity complex (PAR3-PAR6-aPKC) that participates in the control of asymmetric cell division and self-renewal of muscle satellite cells [24]. Although PARD6G is not the nearest protein-coding gene to the chr18 rs2061566 SNP, we queried the STRING interaction database and observed a network (Fig. 5A) connecting the PARD6/3 multiprotein complex to dystrophin through joint interactions with Mark2 (also known as PAR-1b), which is known to bind dystrophin in muscle satellite cells [24] via specific interaction within the 8th and 9th spectrin-like repeats within its central rod domain [25]. Stem cell (including satellite cell) polarity is driven by the localization of this complex (PAR3-PAR6-aPKC) to the apical pole, along with the localization of the PAR1 complex to the basilar pole [26] (Fig. 5B). These complexes are critical for the localization of centrosomes and asymmetric cell division resulting in a population of self-renewed satellite cells and a population leading to myofiber regeneration [26]. Although muscle from DMD patients shows increased numbers of satellite cells, particularly in association with type 1 fibers in biopsies from patients with increased age [27], polarity and associated asymmetric division into regenerative cells is impaired in DMD muscle [26, 28].

Fig. 5. The rs2061566 SNP alters an enhancer linked to PARD6G.

Fig. 5

Open in a new tab

A Two-layer protein network from the STRING database using Mus musculus Mark2 as the query. Nodes are proteins and edges are the connection strengths using ‘Experiments’, “Databases” and “Textmining” as interaction data sources. B Model for postulated roles ETAA1 and PARD6G in control of symmetric/asymmetric cell division in satellite cell niche. Arrows indicate signaling cascades with phosphorylation events indicated by yellow circles. C Sequence features in the rs2061566 region, including the cCRE E192950 enhancer (orange) and transcription factor binding motifs (blue) assigned from index DHS TF footprints. D Cognate site identification scores from a 10-mer sliding window predict increased binding of FOS-MAFB and NFE2L1-MAFG heterodimers to the rs2061566 G allele (orange) versus the A allele (red). Representative MEME motifs are shown for the predicted bZIP TF heterodimers. E Scatterplot of rs2061566 allelic imbalance from heterozygous samples indicate increased TF footprinting on the rs2061566 G allele. F Log-normalized expression levels of Dmd, Pard6g, Etaa1, and myogenic marker genes from mouse single-cell RNA-seq data shown as violin plots, with myogenic states in bins as described in McKellar et al.

A plausible role for PARD6G is in part suggested by functional evidence linking rs2061566 to E1929501 enhancer activity. We used ENCODE transcription factor (TF) footprints [29] to assess whether rs2061566 disrupted a regulatory protein binding site within the E1929501 enhancer. Figure 5C shows that rs2061566 intersects overlapping small MAF (musculoaponeurotic fibrosarcoma) binding sites for homo- and heterodimeric basic region leucine zipper-type (bZIP) TFs. The rs2061566 allelic effect on bZIP TF binding was tested using empirical DNA binding z-scores (CSI scores = cognate site identification 10mer intensity scores) from 109 bZIP TF homo- and heterodimeric pairs [30]. We ranked CSI scores using a 10-mer sliding window across a 21 bp. sequence centered on rs2061566 and observed a maximum CSI binding score for FOS-MAFB bZIP heterodimers at a position containing the alternate rs2061566 G allele (Fig. 5D). Similar binding specificity was also observed for NFE2L1-MAFG heterodimers and the binding specificities of both bZIP TF heterodimers at position 5 were consistent with their representative heterodimer motifs (Fig. 5D). We tested this predicted gain-of-function effect by examining digital TF footprinting data from ENCODE tissue/cell biosamples heterozygous for rs2061566 (Table S7). The allelic imbalance observed for reads overlapping the rs2061566 site (Fig. 5E) suggests that the rs2061566 G allele is associated with increased bZIP TF binding in vivo. The cell type specificity and strength of the E1929501 enhancer – PARD6G promoter interaction relative to other promoters in the region was examined using the activity-by-contact model of enhancer-promoter maps from 131 human cell types and tissues [31], as discussed in more detail in Supplementary Discussion.

The location of rs2061566 in a prominent cCRE enhancer that has strong activity-by-contact chromatin links to the PARD6G promoter in stem cells and the altered transcription factor binding to the heterodimeric MAF motif was the most direct evidence we observed for a functional SNP effect. Single cell transcriptome analysis of muscle stem cells from adult mouse hindlimb identified the bZIP transcription factor MAFF as one of the top-75 differentially expressed genes in the muscle stem/progenitor and myonuclei clusters [32]. Maff expression is highest in quiescent MuSCs suggesting a potential functional effect through rs2061566-altered binding of MAFF heterodimers to this cCRE distal enhancer in these cell types, therefore we examined Maff, Dmd, Pard6g, and Etaa1 expression levels in mouse single-cell RNA-seq data binned by myogenic differentiation status. Maff showed stage-specific expression levels similar to Pax7 and MyoD1, and Pard6g showed increased expression levels that overlapped stage-specific expression of Dmd, Jdp2 and Myh3 (Fig. 5F). Jdp2, a transcription factor with a putative role in regenerative myogenesis, and Myh3, an embryonic myosin heavy chain, are both marker genes for a regenerative myonuclear population enriched in dystrophic muscle from Dmd ∆Ex51 mice [33]. Co-expression of PAR-6 gamma and dystrophin in myogenic cells suggests a common stage for rs2061566-altered enhancer function to exert an effect on disease progression.

Discussion

In the largest genome-wide association study to date of the LOA phenotype in DMD subjects, we report the identification of multiple loci that may harbor genetic modifiers of DMD severity. Our approach was distinct from previous efforts, as it is based on a cohort carefully restricted—using highly conservative criteria—to LoF DMD mutations and uses statistical methods particularly well-adapted to GWAS in moderate sample sizes.

One aspect of these new SNP associations is that we did not detect previously identified modifier genes, which are reviewed in detail elsewhere [3]. The PPLDs at these SNPs were all at or below the prior probability of SNP-trait association (Table S11). For comparison, we also show MOD scores, which are log10 maximum likelihood ratios (where the maximization is over the 6 parameters of the trait model); –log10(p values) from Cox Proportional Hazards (CPH) regression, and from variance-components regression (as we previously reported [4]). The MODs range from 0.90 – 2.98. By contrast, MODs at the SNPs showing the PPLDs ≥ 0.40 (Table 1) range from 6.6 to 8.4. Table 2 clearly indicates that failure to detect these SNPs in the current analysis is not a consequence of our use of the PPLD; rather, evidence for these SNPs is simply not found in this particular cohort.

Of note are the differences between the current findings and our initial report of associations with LTBP4 and THBS1, which was based on a subset (N = 253) of the current cohort. As noted above, 174 subjects overlap between the original and current cohorts: 54 of the original 253 subjects were dropped from the current cohort based on DMD mutation class, 12 were dropped for missing steroid information and an additional 13 were dropped based on other filtering criteria as described above. If we consider just the 67 nonoverlapping subjects with steroid data, we obtain MODs of 4.8, 3.9 for rs710160 (LTBP4), rs2725797 (THBS1) respectively, and variance-components (recessive) –log10(p values) of 6.2, 4.6. Thus these non-overlapping subjects contributed substantially to the original findings. Not surprisingly, detection of any given modifier across patient cohorts with substantial differences (such as in DMD mutation status) remains challenging [12]. Together with previous studies, our results are consistent with the idea that there may be many genes that can modify the DMD phenotype. In this case, simple sampling variability from one dataset to another could make direct replication across studies difficult, as the effects of different genes become salient in different samples, against different genetic and possibly environmental backgrounds, and perhaps under different data analytic strategies [12].

Here we have identified multiple genes suitable for further exploration, based upon multiple lines of evidence for plausibility for a biologic role for the identified SNPs in regulating genes or pathways potentially relevant to muscle biology. A uniform in silico approach provided multiple concordant lines of evidence supporting each of the candidate genes, each of which may play a role in muscle function or pathology. Although the PPLD suggested a distinct set of modifier loci from previously described candidates, the overlap of gene functions in these data sets suggests that recurrent biological pathways are affected. In particular, overlapping functions in muscle stem cell maintenance and division for ETAA1 and PARD6G suggests a model connecting these two genes to dystrophin loss in satellite cells (Fig. 5B). Impairments of the MuSC population in dystrophic muscle have been attributed to exhaustion due to constant activation, to defects in asymmetric cell division and reduction in myogenic precursor cells caused by an intrinsic lack of dystrophin in MuSCs, and to differences in the stem cell niche [24]. In animal models of DMD, such as the Golden Retriever Muscular Dystrophy dog, increased expression of the Notch ligand Jagged1 is associated with a mild phenotype that is postulated to result from increased proliferation of myogenic precursors from MuSCs [9]. In an mdx mouse model engineered with defective telomerase activity, telomere shortening was associated with increased replicative senescence of mdx MuSCs and increased disease severity [34]. In that model, telomere shortening induced markers of the DNA damage response which impaired MuSC differentiation into myotubes in vitro by interfering with MYOD-mediated activation of myogenic gene expression [35]. The SNP effects on ETAA1 expression may suggest a vulnerability in rapidly proliferating MuSCs that activates the ETAA1-ATR signaling pathway during S and M cell cycle progression due to elevated replicative stress or incompletely replicated loci.

Recently, it has been shown that signaling through the Aurora A kinase (AURKA) cascade can rescue some of the polarity and asymmetric satellite cell division defects observed in mdx mice [36]. Stimulation of mouse Aurka activity through the EGFR signaling pathway increases the proportion of asymmetric cell divisions, while pharmacological inhibition of Aurka shifts the balance towards symmetric stem cell division. Prior work established that Aurora A kinase phosphorylation of PARD6G ser [34] is a key step in regulating the composition of the PAR complex including the activation of αPKC and the release of the Lgl subunits. Independent of its role in the PAR complex, PARD6G also has a unique role in the two centrioles of the centrosome, where PARD6G preferentially associates with the older mother centriole to regulate centrosomal protein composition, including the formation of spindle poles that control the accurate segregation of DNA into the two daughter cells [37]. Thus, the regulatory SNPs identified for ETAA1 and PARD6G may play multiple roles in cell polarity and asymmetric division which impact the response to dystrophin loss in satellite cells. The degree of impairment to the eventual deterioration of muscle integrity (and ultimately function) is arguably understudied and may be of increased importance to understand as microdystrophin constructs in gene therapy trials do not encode the Mark2-binding spectrin repeats. Along with the fact that muscle stem cells are poorly transduced in vivo by current AAV vectors and would poorly express transgene with the muscle-specific promoters currently in use, this raises the possibility that microdystrophin gene therapy will not correct fundamental muscle stem cell defects important to muscle maintenance.

Studies that examine the interaction between large-effect monogenic mutations and small-effect modifier variants, such as the one described here, are few in number. It is notable that the associations of common SNPs detected with this study were all non-coding variants, consistent with the majority of GWAS-associated variants linked to human traits. Recent studies have observed that the polygenic background as measured by risk scores for common diseases such as coronary heart disease, breast cancer, and colon cancer, modifies the disease penetrance of hereditary forms of these diseases [38], supporting a ‘many modifiers’ model for monogenic disease phenotypes. That recurrent biological pathways are evident between the previously described DMD modifiers and the new set of potential modifiers we report here suggests that therapeutic strategies based on these pathways requires further understanding of their common physiological action in the context of the disease-causing mutation.

A limitation to our study remains the relatively small cohort size, although our strategy of restricting to total LoF mutations results in a far more homogeneous cohort. The absence of replication of other modifiers, including those identified by us in an overlapping but smaller cohort, is an observation that likely reflects the differential phenotypic and DMD mutation categorizations we used in this study. We note that the eQTL data in public databases was not generated from skeletal muscle but primarily from white blood cell databases, and that the strength of evidence of potential associations for the candidate genes varies. Nevertheless, although experimental confirmation of their biological effects are required—experiments that are beyond the scope of the present study—our results suggest new plausible candidates for genetic modifiers of DMD disease severity as measured by LOA. The data we present here suggest that even in a relatively small data set—albeit the largest studied in DMD to date—our approach to the analysis of DMD phenotype and SNP association may be an effective strategy for prioritizing genes for experimental verification.

Supplementary information

Supplemental Figures (3.1MB, docx)
Supplemental Methods and Discussion (128.8KB, docx)
Supplemental Tables (98.1KB, xlsx)

Acknowledgements

The data used for the analyses described in this manuscript were obtained from the GTEx Portal using the GTEx Analysis Release V8 (dbGaP Accession phs000424.v8.p2). The authors wish to acknowledge the many members of the United Dystrophinopathy Project consortium who participated in the enrollment of subjects in the original historical UDP database.

Author contributions

RBW, VJV, and KMF conceived and designed the study. RBW, VJV, DMD, MW, RA, LA, MM-C, KMF, JB, SCS, and the UDP Investigators acquired and/or analyzed the data. RBW, VJV, PTM, and KMF drafted the manuscript, which was reviewed and revised by all of the named authors.

Funding

This work was supported by the National Institutes of Health (NINDS NS085238) to KMF, RBW, and VJV. The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS.

Data availability

The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests

The authors declare no competing interests.

Ethical approval

The research involving human subjects has been been performed in accordance with the Declaration of Helsinki and was approved by the Nationwide Children’s Hospital Institutional Review Board (IRB) under approval 0502HSE046. Informed consent was obtained from all subjects.

Footnotes

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

These authors contributed equally: Robert B. Weiss, Veronica J. Vieland.

Supplementary information

The online version contains supplementary material available at 10.1038/s41431-023-01329-5.

References

  • 1.Flanigan KM, Dunn DM, von Niederhausern A, Soltanzadeh P, Howard MT, Sampson JB, et al. Nonsense mutation-associated Becker muscular dystrophy: interplay between exon definition and splicing regulatory elements within the DMD gene. Hum Mutat. 2011;32:299–308. doi: 10.1002/humu.21426. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 2.Bello L, Flanigan KM, Weiss RB, Spitali P, Aartsma-Rus A, Muntoni F, et al. Association study of exon variants in the NF-kappaB and TGFbeta pathways identifies CD40 as a modifier of duchenne muscular dystrophy. Am J Hum Genet. 2016;99:1163–71. doi: 10.1016/j.ajhg.2016.08.023. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 3.Bello L, Pegoraro E. The “Usual Suspects”: genes for inflammation, fibrosis, regeneration, and muscle strength modify duchenne muscular dystrophy. J Clin Med. 2019;8:649. [DOI] [PMC free article] [PubMed]
  • 4.Weiss RB, Vieland VJ, Dunn DM, Kaminoh Y, Flanigan KM, United Dystrophinopathy P. Long-range genomic regulators of THBS1 and LTBP4 modify disease severity in duchenne muscular dystrophy. Ann Neurol. 2018;84:234–45. doi: 10.1002/ana.25283. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 5.Heydemann A, Ceco E, Lim JE, Hadhazy M, Ryder P, Moran JL, et al. Latent TGF-beta-binding protein 4 modifies muscular dystrophy in mice. J Clin Invest. 2009;119:3703–12. doi: 10.1172/JCI39845. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 6.Flanigan KM, Ceco E, Lamar KM, Kaminoh Y, Dunn DM, Mendell JR, et al. LTBP4 genotype predicts age of ambulatory loss in Duchenne muscular dystrophy. Ann Neurol. 2013;73:481–8. doi: 10.1002/ana.23819. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 7.Pegoraro E, Hoffman EP, Piva L, Gavassini BF, Cagnin S, Ermani M, et al. SPP1 genotype is a determinant of disease severity in Duchenne muscular dystrophy. Neurology. 2011;76:219–26. doi: 10.1212/WNL.0b013e318207afeb. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 8.Quattrocelli M, Capote J, Ohiri JC, Warner JL, Vo AH, Earley JU, et al. Genetic modifiers of muscular dystrophy act on sarcolemmal resealing and recovery from injury. PLoS Genet. 2017;13:e1007070. doi: 10.1371/journal.pgen.1007070. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 9.Vieira NM, Elvers I, Alexander MS, Moreira YB, Eran A, Gomes JP, et al. Jagged 1 rescues the duchenne muscular dystrophy phenotype. Cell. 2015;163:1204–13. doi: 10.1016/j.cell.2015.10.049. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 10.Crawford SE, Stellmach V, Murphy-Ullrich JE, Ribeiro SM, Lawler J, Hynes RO, et al. Thrombospondin-1 is a major activator of TGF-beta1 in vivo. Cell. 1998;93:1159–70. doi: 10.1016/S0092-8674(00)81460-9. [DOI] [PubMed] [Google Scholar]
  • 11.Spitali P, Zaharieva I, Bohringer S, Hiller M, Chaouch A, Roos A, et al. TCTEX1D1 is a genetic modifier of disease progression in Duchenne muscular dystrophy. Eur J Hum Genet. 2020;28:815–25. doi: 10.1038/s41431-019-0563-6. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 12.Vieland VJ, Seok SC. The PPLD has advantages over conventional regression methods in application to moderately sized genome-wide association studies. PLoS One. 2021;16:e0257164. doi: 10.1371/journal.pone.0257164. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 13.Flanigan KM, Dunn DM, von Niederhausern A, Soltanzadeh P, Gappmaier E, Howard MT, et al. Mutational spectrum of DMD mutations in dystrophinopathy patients: application of modern diagnostic techniques to a large cohort. Hum Mutat. 2009;30:1657–66. doi: 10.1002/humu.21114. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 14.Flanigan KM, von Niederhausern A, Dunn DM, Alder J, Mendell JR, Weiss RB. Rapid direct sequence analysis of the dystrophin gene. Am J Hum Genet. 2003;72:931–9. doi: 10.1086/374176. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 15.Waldrop MA, Gumienny F, El Husayni S, Frank DE, Weiss RB, Flanigan KM. Low-level dystrophin expression attenuating the dystrophinopathy phenotype. Neuromuscul Disord. 2018;28:116–21. doi: 10.1016/j.nmd.2017.11.007. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 16.Flanigan KM, Dunn DM, von Niederhausern A, Howard MT, Mendell J, Connolly A, et al. DMD Trp3X nonsense mutation associated with a founder effect in North American families with mild Becker muscular dystrophy. Neuromuscul Disord. 2009;19:743–8. doi: 10.1016/j.nmd.2009.08.010. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 17.Muntoni F, Gobbi P, Sewry C, Sherratt T, Taylor J, Sandhu SK, et al. Deletions in the 5’ region of dystrophin and resulting phenotypes. J Med Genet. 1994;31:843–7. doi: 10.1136/jmg.31.11.843. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 18.Wang RT, Barthelemy F, Martin AS, Douine ED, Eskin A, Lucas A, et al. DMD genotype correlations from the Duchenne Registry: endogenous exon skipping is a factor in prolonged ambulation for individuals with a defined mutation subtype. Hum Mutat. 2018;39:1193–202. doi: 10.1002/humu.23561. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 19.Shiga N, Takeshima Y, Sakamoto H, Inoue K, Yokota Y, Yokoyama M, et al. Disruption of the splicing enhancer sequence within exon 27 of the dystrophin gene by a nonsense mutation induces partial skipping of the exon and is responsible for Becker muscular dystrophy. J Clin Investig. 1997;100:2204–10. doi: 10.1172/JCI119757. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 20.Disset A, Bourgeois CF, Benmalek N, Claustres M, Stevenin J, Tuffery-Giraud S. An exon skipping-associated nonsense mutation in the dystrophin gene uncovers a complex interplay between multiple antagonistic splicing elements. Hum Mol Genet. 2006;15:999–1013. doi: 10.1093/hmg/ddl015. [DOI] [PubMed] [Google Scholar]
  • 21.Vieland VJ, Huang Y, Seok S-C, Burian J, Catalyurek U, O’Connell J, et al. KELVIN: a software package for rigorous measurement of statistical evidence in human genetics. Hum Hered. 2011;72:276–88. doi: 10.1159/000330634. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 22.Vieland VJ, Seok S-C, Stewart WCL. A new linear regression-like residual for survival analysis, with application to genome wide association studies of time-to-event data. PLoS One. 2020;15:e0232300. doi: 10.1371/journal.pone.0232300. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 23.Miosge LA, Sontani Y, Chuah A, Horikawa K, Russell TA, Mei Y, et al. Systems-guided forward genetic screen reveals a critical role of the replication stress response protein ETAA1 in T cell clonal expansion. Proc Natl Acad Sci USA. 2017;114:E5216–25. doi: 10.1073/pnas.1705795114. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 24.Dumont NA, Wang YX, von Maltzahn J, Pasut A, Bentzinger CF, Brun CE, et al. Dystrophin expression in muscle stem cells regulates their polarity and asymmetric division. Nat Med. 2015;21:1455–63. doi: 10.1038/nm.3990. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 25.Yamashita K, Suzuki A, Satoh Y, Ide M, Amano Y, Masuda-Hirata M, et al. The 8th and 9th tandem spectrin-like repeats of utrophin cooperatively form a functional unit to interact with polarity-regulating kinase PAR-1b. Biochem Biophys Res Commun. 2010;391:812–7. doi: 10.1016/j.bbrc.2009.11.144. [DOI] [PubMed] [Google Scholar]
  • 26.Feige P, Brun CE, Ritso M, Rudnicki MA. Orienting muscle stem cells for regeneration in homeostasis, aging, and disease. Cell Stem Cell. 2018;23:653–64. doi: 10.1016/j.stem.2018.10.006. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 27.Bankole LC, Feasson L, Ponsot E, Kadi F. Fibre type-specific satellite cell content in two models of muscle disease. Histopathology. 2013;63:826–32. doi: 10.1111/his.12231. [DOI] [PubMed] [Google Scholar]
  • 28.Chang NC, Chevalier FP, Rudnicki MA. Satellite cells in muscular dystrophy—lost in polarity. Trends Mol Med. 2016;22:479–96. doi: 10.1016/j.molmed.2016.04.002. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 29.Vierstra J, Lazar J, Sandstrom R, Halow J, Lee K, Bates D, et al. Global reference mapping of human transcription factor footprints. Nature. 2020;583:729–36. doi: 10.1038/s41586-020-2528-x. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 30.Rodriguez-Martinez JA, Reinke AW, Bhimsaria D, Keating AE, Ansari AZ. Combinatorial bZIP dimers display complex DNA-binding specificity landscapes. Elife. 2017;6:e19272. doi: 10.7554/eLife.19272. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 31.Nasser J, Bergman DT, Fulco CP, Guckelberger P, Doughty BR, Patwardhan TA, et al. Genome-wide enhancer maps link risk variants to disease genes. Nature. 2021;593:238–43. [DOI] [PMC free article] [PubMed]
  • 32.De Micheli AJ, Laurilliard EJ, Heinke CL, Ravichandran H, Fraczek P, Soueid-Baumgarten S, et al. Single-cell analysis of the muscle stem cell hierarchy identifies heterotypic communication signals involved in skeletal muscle regeneration. Cell Rep. 2020;30:3583–95.e3585. doi: 10.1016/j.celrep.2020.02.067. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 33.Chemello F, Wang Z, Li H, McAnally JR, Liu N, Bassel-Duby R, et al. Degenerative and regenerative pathways underlying Duchenne muscular dystrophy revealed by single-nucleus RNA sequencing. Proc Natl Acad Sci USA. 2020;117:29691–701. doi: 10.1073/pnas.2018391117. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 34.Sacco A, Mourkioti F, Tran R, Choi J, Llewellyn M, Kraft P, et al. Short telomeres and stem cell exhaustion model Duchenne muscular dystrophy in mdx/mTR mice. Cell. 2010;143:1059–71. doi: 10.1016/j.cell.2010.11.039. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 35.Latella L, Dall’Agnese A, Boscolo FS, Nardoni C, Cosentino M, Lahm A, et al. DNA damage signaling mediates the functional antagonism between replicative senescence and terminal muscle differentiation. Genes Dev. 2017;31:648–59. doi: 10.1101/gad.293266.116. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 36.Wang YX, Feige P, Brun CE, Hekmatnejad B, Dumont NA, Renaud JM, et al. EGFR-Aurka signaling rescues polarity and regeneration defects in dystrophin-deficient muscle stem cells by increasing asymmetric divisions. Cell Stem Cell. 2019;24:419–32.e416. doi: 10.1016/j.stem.2019.01.002. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 37.Dormoy V, Tormanen K, Sutterlin C. Par6gamma is at the mother centriole and controls centrosomal protein composition through a Par6alpha-dependent pathway. J Cell Sci. 2013;126:860–70. doi: 10.1242/jcs.121186. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 38.Fahed AC, Wang M, Homburger JR, Patel AP, Bick AG, Neben CL, et al. Polygenic background modifies penetrance of monogenic variants for tier 1 genomic conditions. Nat Commun. 2020;11:3635. doi: 10.1038/s41467-020-17374-3. [DOI] [PMC free article] [PubMed] [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

Supplemental Figures (3.1MB, docx)
Supplemental Methods and Discussion (128.8KB, docx)
Supplemental Tables (98.1KB, xlsx)

Data Availability Statement

The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.

Articles from European Journal of Human Genetics are provided here courtesy of Nature Publishing Group

RESOURCES