Gene expression profiling of extraocular muscles in primary inferior oblique overaction

Study design

This retrospective study aimed to investigate the differences in gene expression profiles between primary and secondary IOOA. The study adhered to the principles outlined in the Declaration of Helsinki and was approved by The Medical Ethics Committee of Tianjin Eye Hospital (Approval Number: KY-2023041). We collected extraocular muscle (EOM) tissues from ten patients with primary IOOA and ten patients with secondary IOOA. All participants had undergone inferior oblique muscle myectomy surgery, which involved removal of the temporal muscle portion of the inferior oblique. Written informed consent for EOM biopsies was obtained from all patients or their parents. The collected EOM tissues were immediately frozen in liquid nitrogen for subsequent transcriptome sequencing.

Data collection

Demographic and clinical data were extracted from the participants’ medical records between July 2023 and March 2024, focusing on age, gender, and the degree of IOOA. The diagnosis and grading of IOOA were independently verified by two pediatric ophthalmologists, Dr. Hao R and Dr. Zhang W. Primary IOOA (P group) was diagnosed in the absence of superior oblique muscle palsy and without dysfunction in other extraocular muscles. Secondary IOOA (S group) was diagnosed when inferior oblique overaction occurred secondary to superior oblique muscle palsy, resulting in hypo-infraduction during adduction. IOOA severity was assessed based on vertical deviation during adduction with elevation, categorized as follows: +1 (approximately 10°), +2 (20°), +3 (30°), and +4 (40°), with 0 indicating normal (Lim et al., 2014).

The transcriptome sequencing

RNA sequencing was performed by Biomarker Technologies (Beijing, China). The process included sample preparation, library construction, quality control, and sequencing. RNA quantity and quality were assessed using the Nanodrop 8000 Spectrophotometer (ThermoFisher Scientific, Waltham, USA) and the Agilent 2100 Bioanalyzer/LabChip GXr (Agilent Technologies, Santa Clara, USA). cDNA libraries were prepared following the manufacturer’s protocol, and sequencing was carried out using the Illumina NovaSeq 6000 platform.

Detection of DEGs

Raw sequencing data in fastq format were subjected to quality control using FastQC. Clean reads were then aligned to the human reference genome (hg38) using HISAT2. Gene expression levels were quantified as Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) using StringTie. DEGs were identified using DESeq2 with the following thresholds: |log2 fold change (log2FC)|> 0.58 (corresponding to a fold change of ≥ 1.5 or ≤ 0.67) and a false discovery rate (FDR)-adjusted p-value < 0.05. log2FC is used to measure the expression level differences of DEGs between different groups.

Enrichment analysis

Gene Ontology Consortium annotation was enriched using ClusterProfiler. The P-value was adjusted by the Benjamini and Hochberg correction, and GO terms with a false discovery rate (FDR) < 0.05 were considered to be significantly enriched. GO annotation system contains three main branches: Biological Process (BP), Molecular Function (MF) and Cellular Component (CC). Enrichment analyses for BP, CC and MF were performed using ClusterProfiler and ClueGo. The terms obtained from the enrichment results were visualised chord plots and composition ratio diagram. GSEA was processed on all genes based on expression level. Enriched gene sets were identified as FDR < 0.05. Additionally, the PPI network was built based on the DEGs generated in the differential expression analysis. The PPI networks were visualized by STRING (https://string-db.org/).

Demographic characteristics of the sample

A total of 20 temporal inferior oblique specimens were meticulously collected, comprising ten specimens from patients with primary inferior oblique overaction (IOOA) and ten from patients with secondary IOOA. The mean age of individuals with primary IOOA was 4.4 ± 2.5 years, while those with secondary IOOA had a mean age of 4.2 ± 2.3 years. No statistically significant difference in age was observed between the two groups. Gender distribution was similar, with 50% of patients in the primary IOOA group being female and 60% in the secondary IOOA group, though this difference was not statistically significant. Additionally, there was no significant difference in the degree of IOOA between the two groups (Table 1).

All data were analyzed using SPSS 25.0 (IBM Corp., Armonk, NY, USA). Continuous data with a non-normal distribution were presented as median (Q1, Q3) and analyzed using the Mann–Whitney U test. Two-sided p values < 0.05 were considered statistically significant.

Sample identification and differential gene expression in primary group vs. secondary group

A total of 20 samples were processed for transcriptome sequencing, generating 145.66 Gb clean data. At least 6.31 Gb clean data were generated for each sample with minimum 92.50% of clean data achieved quality score of Q30. Clean reads of each sample were mapped to specified reference genome. Mapping ratio ranged from 94.83% to 95.96%. The statistical power of this experimental design, calculated in RNASeqPower was 0.86.

A total of 26,675 genes were identified in the EOM samples. Applying a cutoff of p-value < 0.05 and fold change of ≥ 1.5 or ≤ 0.67, we identified 258 differentially expressed genes (DEGs). Among these, 148 genes were significantly downregulated, while 110 were upregulated in the primary IOOA group compared to the secondary IOOA group (Fig. 1; Table S1). Volcano plots and heatmaps illustrated distinct gene expression patterns between the primary and secondary groups.

Biological pathways enrichment analysis

Enrichment analysis revealed that upregulated DEGs in the primary IOOA group were predominantly associated with immune responses, including positive regulation of chemokine production, regulation of T cell migration, and myelination (e.g., myelination in the peripheral nervous system). Additionally, pathways related to oxygen carrier activity were enriched. Conversely, downregulated DEGs were associated with ion channels (e.g., atrial cardiac muscle cell action potential), negative regulation of signaling receptor activity, and collagen fibril organization (Fig. 2). A detailed list of the specific DEGs associated with each enriched GO biological pathway is compiled in the Table S2.

Cellular components and molecular function enrichment analysis

Analysis of cellular components showed that DEGs in the primary group were enriched in M-bands (e.g., muscle myosin complex, myosin II complex, and A band), ion channels (e.g., voltage-gated potassium and sodium channel complexes), hemoglobin complex, chylomicron, recycling complex, and perinuclear chromocenter (Fig. 3A, Table 2).

Molecular function enrichment analysis indicated that DEGs were significantly associated with muscle alpha-actinin binding, cytoskeletal protein binding, structural molecule activity conferring elasticity, and actinin binding within M-bands. In ion channels, the voltage-gated potassium channel complex was notably enriched. Furthermore, genes related to oxygen carrier activity demonstrated enrichment in oxygen carrier activity, oxygen binding, heme binding, as well as phospholipase D activity and nuclease activity (Fig. 3B, Table 2).

GSEA and PPI analysis validation

GSEA and PPI network analysis further validated the above results of the enrichment analysis of biological pathways, cellular components, and molecular functions. Additionally, it was discovered that some DEGs are associated with peripheral nervous system development and myelin sheath formation (Fig. 4, Table 2).

GO terms related to nerves, muscles and extracellular matrix

Gene Ontology (GO) analysis identified several biological terms linked to nerves and muscles, such as myelination, ion channels, M-bands, and collagen fibril organization. In particular, key DEGs related to myelination, including MBP, MPZ, and PRX, showed elevated expression in the primary IOOA group compared to the secondary IOOA group. Ion channel-related DEGs, such as KCNA5 and KCNE5, were also upregulated, while genes like CACNA1B, SCN3B, SCN5A, and KCNJ3 exhibited a decrease in expression. Additionally, DEGs associated with M-bands, including MYOM2 and MYOM3, were more highly expressed in the primary group than in the secondary IOOA group. On the other hand, DEGs related to collagen fibril assembly, including COL11A1, COL11A2, COL22A1, COL25A1, and COL9A2, were found to be upregulated in the secondary IOOA group (Table 3).

Gene expression patterns in myelination

The upregulation of myelination-related genes in the primary IOOA group is particularly noteworthy. Genes such as MBP (myelin basic protein), MPZ (myelin protein zero), and PRX (peripheral myelin protein 22) were significantly upregulated. MBP is essential for the stability and maintenance of the myelin sheath (Simons & Trotter, 2007), MPZ facilitates myelin-axon adhesion (Lagueny et al., 1999), and PRX is crucial for Schwann cell integrity and myelin repair (Gillespie et al., 1997). The increased expression of these genes implies a higher prevalence of MIFs in the primary IOOA group, which are known for their role in sustaining prolonged muscle contractions (Rudell et al., 2020).

In extraocular muscles, MIFs differ from singly innervated fibers (SIFs) by their unique innervation patterns and functional properties. MIFs, characterized by en grappe terminals, generate slow, sustained tension rather than rapid twitches, which is critical for precise eye movements and gaze stabilization (Hernández et al., 2019). This distinction is supported by our finding of increased myelination-related gene expression in the primary IOOA group.

Furthermore, it is noteworthy that the vascular-ligamentous-nerve strand, unique to the inferior oblique muscle, is a critical anatomical feature that could offer valuable insights into the differential gene expression observed in this study. The interactions between the vascular, ligamentous, and neural components could influence the muscle’s response to abnormal functional demands, possibly leading to compensatory adaptations such as myelination or other molecular changes. Stager et al. described this unique structure, which comprises a combination of blood vessels, ligaments, and nerve fibers, and serves as a conduit for neuromuscular communication and vascular support to the inferior oblique muscle (Stager Jr, McLoon & Felius, 2013; Stager Jr, Dao & Felius, 2015). This structure may play a role in regulating the muscle’s response to various functional and mechanical demands, which could, in turn, influence gene expression patterns.

Dysregulation of ion channels

Our analysis revealed dysregulation of ion channel genes in the primary IOOA group. Specifically, downregulation was observed in genes associated with calcium channels (CACNA1B), sodium channels (SCN5A, SCN3B), and potassium channels (KCNJ3), while upregulation was noted for genes encoding potassium channels (KCNA5, KCNE5). Ion channels are vital for the generation and propagation of action potentials in muscle fibers (Phillips & Trivedi, 2018). The downregulation of CACNA1B, which encodes N-type calcium channels, may lead to reduced neurotransmitter release and diminished endplate potential amplitude (Nudler et al., 2005). Similarly, the downregulation of SCN5A and SCN3B, which are crucial for sodium channel function, could result in decreased sodium influx and impaired action potential initiation (Fahmi et al., 2001).

Conversely, the upregulation of KCNA5 and KCNE5 suggests faster repolarization, which aligns with the characteristics of MIFs that require rapid adaptation to repetitive stimulation. These ion channel alterations reflect the unique electrical properties of MIFs, which are adapted to maintain slow, sustained tension.

KCNJ3 is associated with allowing potassium influx into cells, crucial for regulating heartbeat (Fleischmann et al., 2004), while KCNA5 contributes to restoring the resting membrane potential after depolarization, impacting insulin secretion regulation (López-Vera et al., 2020). KCNE5 acts as an ancillary subunit to voltage-gated potassium channels (Abbott, 2016). The down-regulation of KCNJ3 and up-regulation of KCNA5 and KCNE5 in the primary IOOA group may suggest faster repolarization.

M-band gene alterations

The upregulation of M-band-related genes, MYOM2 and MYOM3, in the primary IOOA group is also significant. The M-band, located at the center of the sarcomere, plays a critical role in managing force distribution during muscle contraction (Lange et al., 2020). MYOM2 and MYOM3 are integral to the structural integrity and functional performance of the M-band. MYOM2, predominantly expressed in fast skeletal muscles, and MYOM3, found in slow muscle fibers and extraocular muscles, contribute to muscle stiffness and resistance to fatigue (Wiesen et al., 2007; Schoenauer et al., 2008). In the context of muscle stiffness, the findings of previous studies are particularly relevant (Kushner, 2006; Agarwal et al., 2016; Altick et al., 2012; Stager Jr, McLoon & Felius, 2013; Stager Jr, Dao & Felius, 2015). These studies have highlighted how changes in muscle stiffness, such as those observed in overacting extraocular muscles, may reflect underlying structural and functional alterations, including increased innervation and fiber adaptation. This framework aligns with our observation of increased MYOM2 and MYOM3 expression, potentially contributing to the heightened muscle tension in primary IOOA. The increased expression of these genes in the primary IOOA group may enhance muscle rigidity and fatigue resistance, supporting the unique contractile properties of extraocular muscles in this condition.

Expression changes of genes related to skeletal muscle contraction and elasticity

The contraction and elasticity of skeletal muscle are related to various proteins within its structure. Among them, actin and myosin are the core proteins responsible for skeletal muscle contraction, and their interaction generates force through the cross-bridge cycle. Titin is one of the most important elastic proteins in muscle, providing elastic recoil during passive stretching and preventing overstretching. Collagen fibers in the extracellular matrix (ECM) also contribute to the passive elasticity of muscle, particularly during stretching (Roberts, 2016).

Our study found no significant differences in the transcriptional levels of key proteins related to muscle contraction, such as actin and myosin, between the primary and secondary groups. However, Wu et al. (2022)’s RNA-seq analysis of congenital esotropia and SOP-associated inferior oblique overaction (IOOA) identified upregulation of three fast-twitch myosin heavy chain (MyHC) isoforms (MYH1, MYH4, and MYH13) and downregulation of one slow-twitch MyHC isoform (MYH3) in the SOP group. This may suggest a strengthening of the inferior oblique muscle in secondary IOOA. However, Wu et al. (2022)’s study included only five cases of SOP-associated IOOA and three cases of congenital esotropia, which limits the sample size, and these findings require further validation. The differences between our study and Wu et al. (2022)’s may be attributed to variations in sample characteristics. Therefore, additional research with a larger sample size is necessary to confirm these results.

We investigated the expression of genes related to muscle stiffness and found that collagen genes, including COL11A1, COL11A2, COL22A1, COL25A1, and COL9A2, play crucial roles in the structure and function of collagen fibers, vital for connective tissues, myotendinous junction stabilization, and cartilage integrity (Annunen et al., 1999; Melkoniemi et al., 2000; Koch et al., 2004; Holden et al., 1999). However, other stiffness-related genes, including titin, actin, and myosin, did not show significant differences between primary and secondary IOOA. The observed reduction in collagen gene expression in the primary IOAA group was unexpected and surprising, given that prior studies have consistently demonstrated elevated collagen content in the inferior oblique muscle of primary IOAA patients (Stager Jr, McLoon & Felius, 2013). This discrepancy might imply that the secondary IOOA group exhibits elevated muscle stiffness, potentially attributable to compensatory adaptations. However, due to the absence of normal controls in our study, we cannot determine the expression changes of stiffness-related genes in the primary group.

Clinical implications and future directions

The distinct gene expression profiles observed in primary versus secondary IOOA provide valuable insights into the underlying mechanisms of these conditions. The increased prevalence of MIFs and the dysregulation of ion channels and M-band proteins in primary IOOA may contribute to the unique clinical features of this disorder. These findings offer potential targets for therapeutic interventions and highlight the need for further research to explore the functional implications of these DEGs.

Study limitations and future research

Despite the valuable insights gained, this study has several limitations. First, the absence of a baseline gene expression for normal inferior oblique muscle in this study limits the interpretation of gene expression changes. Furthermore, the sample size is adequate for preliminary findings, it may not fully reflect the variability within the inferior oblique overaction (IOOA) population. Moreover, the cross-sectional design of the study restricts the ability to establish causative relationships between observed gene expression changes and clinical outcomes. It is important to note that gene expression changes in overacting inferior oblique muscles may not be causative but rather compensatory, potentially reflecting the muscle’s adaptation to altered functional demands. This aligns with previous studies (Wu et al., 2022; Parks, 1972), which emphasize the plasticity of extraocular muscles in response to biomechanical and neural stimuli. To further validate these findings and explore the functional implications of DEGs, future research should include larger cohorts and adopt longitudinal designs. Additionally, incorporating functional assays and molecular imaging techniques will provide a more detailed understanding of how genetic alterations influence muscle function and ocular motility.