Six OP-MG muscle (1 EOM; 5 OOM) and five control muscle (1 EOM; 4 OOM) samples passed data quality analysis (Table 1).
Differential gene expression in orbicularis oculi muscles
To assess differences in relative gene expression between 5 OP-MG and 4 control OOMs, data normalization was performed using the average values for RPLP0 and ACTN2 as reference genes (Additional file, Table S2).
Figure 2 shows the genes that showed >2-fold difference in gene expression levels between OP-MG vs control OOMs (FDR=15%; p<0.01). All the significantly differentially expressed genes were downregulated in the OP-MG muscles. Interestingly, the three most differentially repressed gene transcripts in the OP-MG samples were those genes included in the profiling array because of their altered expression levels reported in post-mortem strabismic EOMs compared with controls [20]. MYH2, which encodes the myosin heavy chain 2a isoform, found in fast-twitch muscle fibres and abundant in both OOMs [13] and EOMs [27], showed 5-fold downregulation compared to controls (p=0.008). VCAN and ANK1, genes related to extracellular remodelling pathways, as well as genes related to muscle atrophy signalling, CAPN3 and MAPK6 (p<0.009) were downregulated in the OP-MG OOMs. VCAN and ANK1 encode components of the extracellular matrix and muscle, respectively, that regulate myoblast fusion/differentiation and attachment to the sarcoplasmic reticulum in regenerating muscle [28]. CAPN3 (calpain 3) is involved in preserving sarcomeric integrity [29]. MAPK6/ERK3 encodes mitogen-activated protein kinase-6, acts as integration points for multiple processes including myogenesis [30]. TFAM, encoding a critical mitochondrial DNA transcription factor [31], was differentially downregulated (p=0.005) in OP-MG.
PPP1R2 which was prioritized as an OP-MG gene and showed significant repression in OP-MG orbital muscles (2.9-fold; p=0.016), encodes a key enzyme in glycogen synthesis, phosphatase-1 regulatory inhibitor subunit R-2. The repression of PPP1R2 transcripts together with DYRK2 (2.7-fold; p=0.003), which phosphorylates glycogen synthase, suggest dysregulation of muscle energy metabolism in OP-MG muscles, but not control muscle.
Overall, we found differential repression of genes involved in the maintenance of the extracellular matrix and sarcomeric stability, as well as mitochondrial biogenesis and muscle metabolism in orbital muscles from OP-MG cases vs normal controls. Interestingly, genes in the latter categories were included in the array based on postulated pathway involvement informed by genomic studies, whereas the former were considered “strabismus-pathway” genes at the time and included as a controls for the clinical phenotype of altered EOMs contractility in OP-MG.
Differential gene-co-expression in orbicularis oculi muscle by phenotype
Pairwise correlations identified 232 strongly correlated gene pairs in OP-MG OOMs and 92 gene pairs in controls (r>0.96, p<0.01). OP-MG and control correlation plots were generated using this data and hierarchical clustering (k-median) identified potentially functional modules [32] visualised in the OP-MG plot that were not seen in the control plot (Figures 3A vs B). Figure 3C shows the visual network of co-expressed gene pairs in the OOMs of OP-MG samples (from Figure 3A) and the most interconnected genes (10-15 gene-gene interactions each) involved putative OP-MG genes (CANX, DDX17, FAM69A, FAM92A1). DDX17 is a master transcriptional regulator involved in muscle gene expression [33]. CANX encodes calnexin, a chaperone protein that facilitates the assembly of AChR subunits [34], FAM92A1 has recently been recognized for its importance in mitochondrial functioning [35] and both FAM69A and PEF1 have poorly characterized biological function. PEF1 encoding Peflin, which may be involved with endoplasmic reticulum-golgi transport and/or calcium binding [36], connected with a ‘module’ of genes (Table 2) which were profiled because they showed altered muscle expression levels in EAMG [17], and they are known to regulate mitochondrial biogenesis and oxidative metabolism. PARK2, encoding parkin, a major regulator of mitophagy [31], appeared to form a subnode of negatively correlated genes, several of which are involved in muscle regeneration (e.g. DDX17) and mitochondrial biogenesis. PPP1R2, which was differentially repressed in OP-MG muscle (Figure 2), did not correlate with other transcripts in OP-MG muscle yet showed significant correlations with four genes in control OOMs (data not shown).
Next, a differential co-expression algorithm was used to statistically identify gene pairs differing significantly and inversely between the two phenotypes (Fisher’s Z tests, p<0.002). An example is showed in Figure 4A with significant positive correlations in a gene pair in one phenotype (OP-MG or controls) but inverse correlations in the other (p<0.01)(see Additional file: Figure S1). Figure 4B displays a correlation network of the significant inversely correlated gene pairs for the OP-MG network (implying that the control network showed opposite correlations). The importance of this analysis is that it shows inverse gene-gene interconnectivity (correlations) in OP-MG orbital muscle vs controls, of the genes postulated to be OP-MG genes based on previous unbiased genomic comparisons viz. PEF1, FAM92A1, AKT1, MYL12B, and SH3BGR [6]. Interestingly, the gene product of SH3BGR is involved with myosin heavy chain isoforms and sarcomere function and is highly expressed in EOMs.
Figure 4C depicts the representative biological categories of strongly correlated gene pairs in the OP-MG network and inverse in controls. Of the most differentially correlated gene pairs which were inverse by phenotype (n=13; p0.001), >50% involved a gene whose transcripts were substantially altered in EAMG [17] suggesting altered MG pathways in the muscle transcriptome of OP-MG vs controls. The broad category of ‘muscle atrophy’ signalling was the most representative functional category accounting for two-thirds of the most differential and inversely correlated gene pairs by phenotype with oxidative metabolism the next most frequent category (p0.001) (genes may be implicated in 2 categories).
Downregulated genes in OP-MG orbital muscle may be regulated by microRNAs
In OOMs the significant differential downregulation of transcripts in OP-MG suggest post-transcriptional gene repression by microRNAs (miRs) binding to 3’ regulatory variants in OP-MG as a potential regulatory mechanism. Therefore, the miRTarBase database [37] was interrogated for interactions between miRs previously shown to be highly expressed in EOMs [38], and genes differentially downregulated ≥1.5-fold in OP-MG, compared with controls. MiR-499 and miR-206, the most highly expressed miRs in EOM (>7-fold higher than in skeletal muscle), potentially interact with three of the genes which showed differential repression in OOMs (Figure 4D). In addition, CANX and IL6R, which feature prominently in the differential gene-gene cross-correlations (Figure 4B), and predicted to have EOM-miRs interactions, were prioritized previously to have 3’ regulatory region variants present in the OP-MG, but not controls [5,7].
Extraocular muscle
Although only two EOMs (medical recti) passed quality control we briefly mention the results due to scarcity of data on these tissues from live donors. For the EOMs, ACTB was used for normalization (SD=0.06 between samples; Additional file: Table S2). Of relevance to the interpretation of the OOM results, the normalized expression levels of the 125 genes were similar between the two EOMs and nine OOMs (p=0.72), although the EOM-specific isoform, MYH3, differed significantly (p<0.0001; Table S3, Figure S2). As the comparison of gene expression between only two samples required cautious interpretation a quantile-quantile plot of the transcripts of the two samples showed a similar distribution of the majority of the genes with only isolated values >10-fold difference (Figure 5A, B). We briefly discuss WIF1, IL6 and MT1A transcript levels which were ≥25-fold higher in OP-MG compared to control.
WIF1 encodes Wnt inhibitory factor-1, and as an inhibitor of Wnt signalling may affect muscle regeneration and endplate AChR clustering [39]. MT1A encodes metallothionein-1A and protects cells against oxidative stress [in[17]]. Increased expression of metallothionein 1 has been reported in atrophic muscle, possibly by preventing the upregulation of the insulin growth factor (IGF)-1 pathway [40], and in the EOMs of the EAMG rodent models [17]. IL6 encodes interleukin-6, a pivotal myokine in satellite cell/myoblast proliferation amongst others, but with persistently elevated levels it facilitates muscle atrophy and inhibits IGF-1 signalling [41]. Therefore, the higher transcript levels of IL6 and MTIA in the chronically paralysed EOM of an OP-MG vs a non-paralytic control, could be compatible with severely atrophic muscle and/or MG-induced effects [16].