DNA Methylation Analysis Reveals Distinct Patterns in Satellite Cell–Derived Myogenic Progenitor Cells of Subjects with Spastic Cerebral Palsy

Spastic type cerebral palsy (CP) is a complex neuromuscular disorder that involves altered skeletal muscle microanatomy and growth, but little is known about the mechanisms contributing to muscle pathophysiology and dysfunction. Traditional genomic approaches have provided limited insight regarding disease onset and severity, but recent epigenomic studies indicate that DNA methylation patterns can be altered in CP. Here, we examined whether a diagnosis of spastic CP is associated with intrinsic DNA methylation differences in myoblasts and myotubes derived from muscle resident stem cell populations (satellite cells; SCs). Twelve subjects were enrolled (6 CP; 6 control) with informed consent/assent. Skeletal muscle biopsies were obtained during orthopedic surgeries, and SCs were isolated and cultured to establish patient–specific myoblast cell lines capable of proliferation and differentiation in culture. DNA methylation analyses indicated significant differences at 525 individual CpG sites in proliferating SC–derived myoblasts (MB) and 1774 CpG sites in differentiating SC–derived myotubes (MT). Of these, 79 CpG sites were common in both culture types. The distribution of differentially methylated 1 Mbp chromosomal segments indicated distinct regional hypo– and hyper–methylation patterns, and significant enrichment of differentially methylated sites on chromosomes 12, 13, 14, 15, 18, and 20. Average methylation load across 2000 bp regions flanking transcriptional start sites was significantly different in 3 genes in MBs, and 10 genes in MTs. SC derived MBs isolated from study participants with spastic CP exhibited fundamental differences in DNA methylation compared to controls at multiple levels of organization that may reveal new targets for studies of mechanisms contributing to muscle dysregulation in spastic CP.

While all individuals with spastic CP have some movement dysfunction, there is a high degree of variability in phenotype between individuals. A more thorough understanding of the mechanisms associated with dysfunction in the peripheral neuromotor system is needed in order to develop more targeted and enhanced therapeutics addressing the major challenges facing individuals with spastic CP. Importantly, spastic CP is associated with alterations in the muscle-resident stem cell populations (satellite cells; SC) responsible for skeletal muscle growth and repair. Surgical patients with CP have a reduced SC population, which may account for aspects of their impaired muscle growth and decreased ability to strengthen muscle [28][29][30]. Isolated SC-derived myoblasts (MB) from CP donors appear to have altered phenotypes in culture compared to control MBs [28,31,32], and an RNASeq study from our group showed differential expression of mRNAs and miRNAs in spastic CP muscle [26]. These data suggest that MBs from CP patients may retain intrinsic differences through the cell isolation and culture process.
Genetic alterations may account for intrinsic differences in isolated MB populations. A number of potentially causative genetic variants have been identified in CP [33], but not all types of CP are easily detected or characterized by genomic data. Several rare copy number variants and mutations have been identified, but there is considerable genetic heterogeneity in patients with CP [33,34]. While some CP cases may be associated with certain types of genetic abnormality [35], the conventional view of CP remains that environmental factors affecting neuromotor maturation are responsible for most cases, especially among individuals with spastic CP [36,37].
Epigenetic modification may also account for some retained phenotypic differences in isolated MB behavior. In recent work, DNA methylation pattern differences were identified in peripheral blood cells from subjects with CP [38][39][40][41], and some early studies indicate that muscle in CP may be similarly altered, with differential DNA methylation in CP resulting in a decreased capacity for MBs to fuse and differentiate into MTs [32]. It has been demonstrated in several disease states that patterns of altered DNA methylation may uncover molecular etiologies and reveal potential therapeutic targets [42][43][44][45][46][47][48][49][50][51], and such may be the case in spastic CP as well [39]. DNA methylation changes may serve as markers for diagnosis, prognosis, tailoring the best treatment for a subclass of disease, monitoring treatment efficacy, and identifying genes to be examined for the development of genetically or epigenetically targeted therapies [52]. In the current study, DNA methylomes were analyzed to provide insight into individual CpG site differences and altered DNA methylation patterns in chromosomal segments and near transcription start site (TSS) in spastic CP SCs compared to controls.

Subject Enrollment
Six subjects with a diagnosis of spastic CP and 6 control subjects were enrolled in an IRB-approved study at Nemours Children's Hospital, Delaware, after informed consent/assent. The control cohort comprised children with an idiopathic condition or an injury. Subjects with a chromosomal disorder, degenerative neurological disease, or muscular dystrophy were excluded.

Satellite Cell Isolation
Skeletal muscle biopsies collected during orthopedic surgeries were enzymatically digested and a double-immunomagnetic isolation approach was used to collect a population of mononuclear Cells positive for the surface markers neural cell adhesion molecule 1 and C-X-C motif chemokine receptor 4 (anti-NCAM1 and anti-CXCR4, both at 2.5 ng/mL, Miltenyi, San Diego, CA, USA) as previously described [26]. Previous studies have demon-strated that CXCR4 marks human SCs and that selection using the combination of NCAM1 (CD56) and CXCR4 more effectively removes non-satellite cells than using either marker alone [53,54]. This method resulted in a nearly pure SC population as verified by positive PAX7 immunofluorescence signal obtained after 24-48 h in culture (Anti-PAX7 from hybridoma cells deposited to the DSHB by Kawakami, A., Developmental Studies Hybridoma Bank, Iowa City, IA, USA). Cell populations that were at least 90% positive for PAX7 expression were utilized for experiments (Supplemental Table S1).

DNA Extraction, Library Preparation, and Sequencing
Genomic DNA was isolated using Gentra Puregene kits (Qiagen, Germantown, MD, USA). A previously published DNA methylation assay [39,55,56] was utilized. Briefly, DNA libraries for next generation sequencing (NGS) were prepared by digesting genomic DNA with methyl-sensitive restriction endonuclease HpaII, which recognizes CCGG sites. A standard sequencing protocol was then performed including randomized shearing (Covaris, Woburn, MA, USA) and synthesis of a gDNA fragment library using Illumina TruSeq Nano library synthesis kits (San Diego, CA, USA). NGS was performed on an Illumina ×10 platform by Psomagen (Rockville, MD, USA). The protocol generated single end reads (150 bp) with >20× coverage of the regions captured. FASTQ data files were processed to calculate the probability of methylation at individual CpG sites through a commercial bioinformatics pipeline and software platform (Genome Profiling LLC, Newark, DE, USA). For convenience, the term "CpG" in this paper refers to "C(CpG)G" HpaII restriction sites.

Methylation Analysis
FASTQ files were aligned to human reference genome hg19 using BWA (Burrows-Wheeler Aligner algorithm [57]. To reduce false discovery associated with the inclusion of both male and female participants in the study, sites on the X and Y chromosomes were excluded from analyses. For each CpG per sample, a methylation score was calculated proportional to the probability that a specific CpG site was methylated. These scores were then compared between cohorts using a set of analytic modules from the R-packages edgeR [58,59] and limma [60] to compare pairwise for each CpG site. The response scale of the methylation data sets used here is within the operational boundaries of log-scaled gene expression data for which edgeR and limma were designed, and the well-developed false-discovery rate calculations in these R packages are ideal for methylation score data distributions [61,62]. Statistical significance at the level of individual CpG sites and was evaluated using a Likelihood Ratio Test with a one-way ANOVA contrast (LRT-ANOVA). Potentially informative CpG sites were selected for the PCA plots by filtering the LRT-ANOVA p-values with an appropriate cutoff (<0.01). Gene annotations were derived from the Ensembl GRCh37 database based on chromosomal locus (www.ensmbl.org; last accessed on 15 November 2022 for validation). Enhancer regions were derived from Ensembl's GRCh37 BioMart tool. The ontology terms Muscle Organ Development, Muscle System Processes, and Skeletal Muscle Cell Differentiation (GO Enrichment Analysis; amigo.geneontology.org) were used to determine 575 unique genes annotated to be involved in muscle physiology. Fisher's exact test was utilized to determine enrichment of significant CpGs within chromosomes.
Because methylation of the region around the TSS of a gene is thought to be highly informative of gene expression [63], the individual CpG methylation scores were averaged per TSS using 1000 bp upstream of the TSS and 1000 bp downstream [64]. TSS were identified from Ensembl's GRCh37 BioMart data mining tool (release 106) [65]. If one or fewer CpGs were found within this range, the TSS was excluded from further analysis. For each sample, the mean methylation load in this 2000 bp range was calculated and a likelihood ratio test performed on the methylation loads. No more than one transcript for each gene was included in the statistical analysis.

Results
Samples from 12 unique study participants were included in the study; demographic data are summarized in Table 1. Briefly, the control group consisted of n = 6 subjects (males = 3, and females = 3) with an average age of 13.9 ± 1.7 years; the CP group consisted of n = 6 subjects (males = 3, and females = 3) with an average age of 15.5 ± 3.0 years. Biopsies from different muscles were included based on the availability of viable muscle tissue suitable for cell isolation and to help identify methylation signals generally associated with spastic CP rather than a specific muscle: SCs were isolated from spinalis or semitendinosus muscle for controls and from spinalis, rectus femoris, adductor longus, or vastus lateralis for the CP group. Although differences likely exist between muscle types, previous studies of MTs derived from different muscles found low inter-muscle variability in RNA-sequencing data [26]. To evaluate differences in DNA methylation patterns associated with affected muscle of subjects with spastic CP, whole genome methylation patterns were determined from proliferating MBs and differentiating MTs. Cells derived from subjects with CP appeared morphologically similar to those derived from controls (Supplemental Figure S1) and exhibited no significant differences in proliferation rates (Supplemental Table S1). NGS was performed after methylation sensitive restriction endonuclease (HpaII) digestion. The hg19 reference genome assembly from the University of California Santa Cruz [66] includes 2.29 × 10 6 HpaII target CCGG motifs, which represent~15% of the 14 × 10 6 CpG sites in the haploid hg19 genome [39]. Alignment of the HpaII restricted sites in our 12 samples yielded 1,483,038 sites were in common across all subjects for MBs and MTs. DNA methylation patterns were analyzed at the individual CpG site level using dimensionality reduction by principal component analysis (PCA) to assess the degree of discrimination between CP and non-CP cohorts. All potentially informative CpG sites (n = 20,254 for MBs and 27,834 for MTs), were integrated as one pattern and demonstrated strong discrimination based on diagnosis ( Figure 1, Supplemental Figure S2).
To evaluate differences in DNA methylation patterns associated with affected muscle of subjects with spastic CP, whole genome methylation patterns were determined from proliferating MBs and differentiating MTs. Cells derived from subjects with CP appeared morphologically similar to those derived from controls (Supplemental Figure S1) and exhibited no significant differences in proliferation rates (Supplemental Table S1). NGS was performed after methylation sensitive restriction endonuclease (HpaII) digestion. The hg19 reference genome assembly from the University of California Santa Cruz [66] includes 2.29 × 10 6 HpaII target CCGG motifs, which represent ~15% of the 14 × 10 6 CpG sites in the haploid hg19 genome [39]. Alignment of the HpaII restricted sites in our 12 samples yielded 1,483,038 sites were in common across all subjects for MBs and MTs. DNA methylation patterns were analyzed at the individual CpG site level using dimensionality reduction by principal component analysis (PCA) to assess the degree of discrimination between CP and non-CP cohorts. All potentially informative CpG sites (n = 20,254 for MBs and 27,834 for MTs), were integrated as one pattern and demonstrated strong discrimination based on diagnosis ( Figure 1, Supplemental Figure S2). In MBs, 525 CpG sites were found to have differential methylation load scores (FDR < 0.05; heatmap in Figure 2A, volcano plot in Supplemental Figure S3A, list of significant CpGs in Supplemental Table S2). Of these, 11 were within genes known to be involved in muscle physiology and 21 were within known gene enhancer regions. 1774 CpG sites were found to have differential methylation load scores in MTs between the CP and control cohorts (heatmap in Figure 2B, volcano plot in Supplemental Figure S3B, list of significant CpGs in Supplemental Table S3), of which 43 CpG sites were within genes known to be involved in muscle physiology and 97 were within gene enhancer regions. The differentially methylated CpGs included 79 CpG sites that were significantly different under both cell conditions. Of these, 36 were significantly hypermethylated and 43 were significantly hypomethylated in the CP cohort compared to the control cohort ( Table 2).  Table S2). Of these, 11 were within genes known to be involved in muscle physiology and 21 were within known gene enhancer regions. 1774 CpG sites were found to have differential methylation load scores in MTs between the CP and control cohorts (heatmap in Figure 2b, volcano plot in Supplemental Figure S3B, list of significant CpGs in Supplemental Table S3), of which 43 CpG sites were within genes known to be involved in muscle physiology and 97 were within gene enhancer regions. The differentially methylated CpGs included 79 CpG sites that were significantly different under both cell conditions. Of these, 36 were significantly hypermethylated and 43 were significantly hypomethylated in the CP cohort compared to the control cohort (Table 2).
To determine the methylation differences over larger regions of the genome, the chromosomal distribution of significant CpGs was visualized (Supplemental Figure S4) and Fisher's exact test was employed to analyze enrichment of differentially methylated CpGs on individual chromosomes. For both MB and MT cell populations, significant enrichment was found on chromosomes 12, 13, 14, 15, 18, and 20. Interestingly, the same chromosomes were identified as having significant CpG site enrichment in previous studies of both muscle tissue and peripheral blood cells from subjects with CP (Table 3). To further assess regional differences, methylation load scores were calculated across 1 Mbp chromosomal segments. There was a strong correlation between the 1 Mbp methylation loads of MBs and MTs (Figure 3a), indicating stability of the methylome. When these 1 Mpb regions were mapped to the chromosomes, regions of accentuated differential methylation load were noted on all chromosomes except 1 and 17 (Figure 3b,c). These large-scale changes in DNA methylation could affect higher-order chromatin structure and regulation of gene expression [52].    Since DNA methylation in promoter regions of genes has been associated with regulation of expression, the individual CpG methylation scores were averaged across TSS flanking regions using 1000 bp upstream of the TSS and 1000 bp downstream. The promoter methylation data revealed distinct patterns between the control and CP cohorts for both MBs (Table 4) and MTs (Table 5). Of 31,844 unique promoters identified, there were 3 promoters with statistically different methylation loads between CP and control subjects in the MBs and 10 promoters in MTs (FDR < 0.05). The majority of the differentially methylated promoters were in non-coding genes, with only one of the MB promoters and two of the MT promoters associated with protein coding transcripts. To explore the relationship between methylation loads and RNA expression, the correlation of methylation load in protein coding genes was analyzed against previously published RNA-seq count data per gene, but no correlation was found (data not shown).

Discussion
This study found that DNA methylation patterns in skeletal muscle SCs grown in culture differed significantly between a cohort of study participants with spastic CP and non-CP controls. DNA methylation is a common and widespread chemical modification involving the addition of a methyl group to the 5-carbon position of cytosine, predominantly within CpG dinucleotides [67]. DNA methylation patterns can change during normal developmental processes, and it has been shown that altered DNA methylation can be passed to daughter cells and sustained later in life [68][69][70][71][72][73][74][75]. Specific DNA methylation changes can modify gene expression, and DNA methylation is well known to be involved in X-chromosome inactivation, gene imprinting, and the silencing of transposable elements [76]. Changes in DNA methylation patterns can also occur as a result of pathophysiologic processes or acute exposures to environmental or physiologic stress [77][78][79]. Altered DNA methylation has been linked to a number of risk factors and potential causes for CP including prematurity, hypoxia-ischemia, and infection [80][81][82][83]. In a prior study, we found that DNA methylation patterns in peripheral blood cells of spastic CP patients varied significantly from controls [39], raising the possibility of methylome alterations in both hematopoietic stem cell and myogenic stem/SC lineages in spastic CP. Furthermore, epigenetic patterns from adolescents were able to be used to predict diagnosis of much younger patients [39], suggesting that at least some methylation pattern differences are associated with the onset of CP and are preserved over time.
The present study examined methylation pattern differences but did not look directly at differences in cell phenotype or behavior. Studies show that individuals with CP have significantly reduced numbers of SCs and that MBs derived from these SCs have a decreased capacity to fuse and differentiate into MTs in culture [29,84]. A recent report indicated that SC-derived MB progenitors from contractured muscle in CP have globally hyperme-thylated DNA and gene expression patterns that favor proliferation over quiescence and differentiation; in that study, a 24 h treatment with a hypomethylating agent reduced DNA methylation to control levels and promoted an exit from mitosis [32]. While previous studies of DNA methylation in CP SC-derived MBs have either examined DNA methylation of a specific CpG island [28] or used Infinium Human MethylationEPIC Beadchip arrays of 850,000 targeted CpGs to identify hypermethylated regions [32], the current study used a different sequencing technology and computational pipeline to examine over 1.4 million CpGs distributed throughout the genome in both MBs and MTs. In addition, the current study was able to take advantage of more closely age-matched samples than had been possible previously. This combination of more closely matched controls and a broader technology platform allowed for identification of both hypermethylated and hypomethylated regions as well as individual CpGs, which provide higher likelihood candidates for biomarker platforms.
The identification of differentially methylated CpG sites and regions in common across cells from the CP cohort suggests that fundamental molecular alterations associated with diagnosis were sustained after the cells were removed from the structural, biomechanical, and humoral environments of the muscle tissue. Such an effect has been described as muscle epigenetic memory wherein DNA methylation is stably altered by prior events like biomechanical loading or acute early life exposure to inflammatory cytokines [85]. Here, we identified 525 CpGs in MBs and 1774 in MTs that were differentially methylated in spastic CP versus controls. Interestingly, MBs demonstrated similar proliferation rates and RNAseq profiles between cohorts in a previous study [26], indicating that while these 525 CpGs may be biomarkers for CP, they may not be associated with functional changes within the cells. Further work is needed to elucidate fully specific linkages between DNA methylation and the regulation of protein levels and cell activities. The larger number of significant sites in MTs was consistent with a larger number of differentially expressed genes in RNA-seq and may indicate that a mixture of cells at different stages of myogenic differentiation was present at the time of DNA isolation. In addition, reports indicate that native SC populations actually comprise multiple different subpopulations that may differentially contribute to variability; studies focused on clonal cells rather than heterogenous populations may be needed. Interestingly, in both MBs and MTs, several significant CpG sites were within genes known to be involved in muscle physiology, including skeletal muscle differentiation (HLF, NOTCH1), muscle organ development (BMP2, COL6A3, DCN, FZD2, HEG1, HLF, ITGA11,  ITGA7, LAMA5, LARGE1, MAPK14, MYBPC1, MYH6, MYLK, NOTCH1, NRG1, PKP2,  SGCD, SMAD7, TBX1, TCF12, TEAD4 (Table 2). Alterations in DNA methylation can be sustained long-term and previous studies indicated that the majority of the DNA methylome remained relatively preserved through myogenesis, from SC to MT formation [76,86]. These 79 sites may therefore represent stable methylation signals indicative of CP; however, more studies are needed to determine the implication(s) of these differentially methylated sites and their roles in muscle impairment in CP.
Chromosome enrichment analysis determined that there was an enrichment of significant CpGs on chromosomes 12,13,14,15,18, and 20 for both MBs and MTs (Table 3). Interestingly, significant CpGs were also enriched on the same chromosomes in the CP cohort in skeletal muscle tissue and blood cells. An analysis of differential methylation over 1 Mpb regions in MBs demonstrated that chromosomes 6,9,11, and 21 contained regions of hypermethylation in the CP cohort, while chromosomes 5,7,8,16,19, and 20 contained regions of hypomethylation, and chromosomes 2, 3, 4, 10, 12, 13, 14, 15, and 18 contained regions of both hypermethylation and hypomethylation. MTs contained more hypermethylated regions and fewer hypomethylated regions in the CP cohort than MBs.
Chromosomes 4, 8, 9, 11, 12, 14, 16, 21, and 22 contained regions of hypermethylation, while chromosome 7 contained regions of hypomethylation, and chromosomes 2, 3, 5, 6, 10, 13, 15, 18, and 20 contained regions of both hypermethylation and hypomethylation ( Figure 3). Overall, the differential methylation levels over 1 Mbp regions were well correlated between MBs and MTs, again suggesting stability of the methylome during the course of myogenesis. The large regions of differential methylation between the CP and control cohorts throughout the genome suggest differences in chromatin structure within the CP cohort as various chromatin states based on histone modifications and nucleosome positioning can determine DNA methylation patterning [67]. There are complex mechanisms underlying the molecular crosstalk between DNA and histone methylation [87], and additional studies are needed to investigate specific histone modifications in CP to understand this complex relationship. While we are in the early stages of unraveling how these alterations are relevant to CP, several key chromosomes were identified as potential targets for future investigation in spastic CP.
We also identified statistically significant differences in methylation of the promoter regions of genes and assessed the relationship between these differences and an RNA-seq study of the same samples [26] to investigate the effect of differential DNA methylation on gene expression. Of particular note, specific associations between protein coding DNA methylation patterns and RNA expression were not readily resolved in our study. However, methylation/expression relationships are difficult to resolve in general and several studies have demonstrated that relationships between methylation status and gene expression can be complex [88,89]. Additionally, there is no current gold standard method for rolling up methylation scores across groups of individual CpG sites into a relevant burden for individual genes. Furthermore, approximately 95% of CpG island promoter regions are unmethylated independent of gene activity and recent studies suggest that methylation of promoter CpG islands is not the primary determinant of gene activity [52,90]. Many CpG islands occur in gene bodies, intergenic regions, or enhancers and may be relevant to gene expression [52]. In fact, recent studies have suggested that altered methylation in enhancer regions rather than promoter regions may be more indicative of changes in gene expression [52,[91][92][93]. Enhancers can regulate the transcription of one or more genes, regardless of orientation or relative distance to the target promoter [94]. While 21 differentially methylated CpGs were identified within annotated enhancer regions in MBs and 97 in MTs, enhancers are difficult to map experimentally [52], enhancer activity is context and stimulus-dependent, and there are no genome-wide enhancer sets linked to specific promoters [94]. It was not possible to investigate the effects on gene expression due to these limitations in the annotation of enhancer elements. Therefore, it will be essential to continue investigating genome-wide analysis approaches that can accurately associate high-throughput expression data with methylation signatures.
Of note, the majority of the differentially methylated promotors identified in MBs and MTs were for regulatory RNAs; a result that may be indicative of differences in RNA processing within cells isolated from CP tissue (Tables 4 and 5). These regulatory RNAs comprise the majority of the transcriptome and play critical roles in maintaining gene expression regulation. The most well-known regulatory RNAs include micro RNAs and long non-coding RNAs. MicroRNAs (miRNAs) are small noncoding RNAs (~22 nucleotides in length) that play important roles in developmental processes such as myogenesis and neurogenesis [95]. Long non-coding RNAs (lncRNAs) may mediate chromatin remodeling and modification, interact with transcription factors for gene regulation, interact with mRNAs to regulate post-transcriptional processes [96], and interact with miRNAs to facilitate myogenesis [97]. Additional studies are warranted to investigate the complex interplay between DNA methylation, histone modifiers, and non-coding RNAs in order to provide a comprehensive understanding of these epigenetic modulations on SC physiology and myogenesis in spastic CP.
The findings of our study compellingly support the idea that spastic CP is associated with altered epigenetic pathways, but our studies are limited by our reliance on the ability to obtain biopsies from individuals presenting for surgery. Because of this, our number of samples is small, SCs must be derived from different muscles in order to age-match samples, and our CP population presenting for surgery consists of mostly severely affected individuals with a high level of motor impairment and inactivity (GMFCS V); therefore, extrapolation of our results to a larger CP community may require additional research. Additionally, the technology used interrogates 2.29 × 10 6 CCGG motifs, which represent 15% of the 14 × 10 6 CpG sites in the genome. The majority of the CpG sites in our study do not overlap with those selected for inclusion in the Infinium MethylationEPIC technology making comparisons to other studies using Infinium data challenging.

Conclusions
In this report, an innovative DNA methylation analysis was employed on SC-derived MBs and MTs collected from individuals with and without CP. We identified differential methylation in the CP cohort at the levels of individual CpGs, 1 Mpb regions, and promoters. The work presented here leverages our novel methylation approach with ex vivo cell studies to elucidate aberrant methylation signatures.
Supplementary Materials: The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/jpm12121978/s1, Figure S1: Cell morphology, Figure S2: Scree plots for explained variance, Figure S3: Volcano plots, Figure S4: Manhattan plots, Table S1: Cell culture information,   Data Availability Statement: Methylomic data will be made publicly accessible via the controlled access system provided by the NIH database of Genotypes and Phenotypes (dbGaP) or can be made available upon a direct request to the corresponding author.

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