Downregulated MYL1 Emerges as a Promising Diagnostic Biomarker for Rheumatoid Arthritis

Introduction

Rheumatoid arthritis (RA) is a prevalent autoimmune joint disorder that imposes a substantial global burden.¹ In 2020, approximately 17.6 million individuals were affected by RA worldwide, corresponding to an age-standardized prevalence of 208.8 cases per 100,000 population.² This represents a notable 14.1% increase compared to the prevalence reported in 1990.² RA is clinically characterized by joint pain, swelling, deformity, morning stiffness, and extra-articular manifestations.^3,4 The etiology of RA is multifactorial and involves autoimmune dysregulation, genetic predisposition, and potential infections. Currently, drug therapy for RA mainly aims to relieve symptoms and control disease progression, whereas surgical interventions are typically reserved for end-stage disease to alleviate pain and improve joint function.^1,5,6 Current biomarkers—such as rheumatoid factor, anti‑cyclic citrullinated peptide (anti‑CCP) antibodies, erythrocyte sedimentation rate (ESR), and C‑reactive protein (CRP)—do not achieve sufficient sensitivity or specificity in all patients, particularly those with seronegative RA, underscoring unmet diagnostic needs. Despite significant advances in biologic and targeted synthetic DMARD therapies, many RA patients still fail to achieve sustained remission or experience adverse effects, highlighting persistent therapeutic challenges.

Among the various RA pathological changes, the most critical one is the disruption of immune tolerance, resulting in the immune system erroneously identifying self-tissues as foreign antigens.⁷ Immune complexes then activate synovial cells and trigger them to proliferate and release various inflammatory mediators including tumor necrosis factor-alpha (TNF-α), interleukin-1 (IL-1), and interleukin-6 (IL-6).^8–14 Simultaneously, immune cells such as T-cells and B-cells infiltrate synovial tissues, further exacerbating the progression of synovial inflammation. Under the stimulation of inflammation, the synovial tissue forms a pannus—an abnormal vascular proliferation. The pannus tissue grows over the joint surface and invades the cartilage and bone, eventually leading to joint destruction—the key pathological characteristic of RA.^15,16

Moreover, beyond joints, people with RA frequently develop skeletal muscle atrophy and weakness,^17,18 suggesting a muscle–joint–immune axis that may influence disease activity. The gene family, myosin light chains, is integral to the processes of muscle contraction and cellular motility. Among these genes, myosin light chain 1 (MYL1) stands out with its unique structural features, including the EF-hand domain, which functions as a calcium sensor by binding to Ca²⁺ ions, thus playing a regulatory role in muscle contraction.¹⁹ These processes are increasingly recognized as modulators of immune and stromal cell behavior: cytoskeletal tension can influence synoviocyte migration and invasiveness;²⁰ calcium-dependent signaling helps regulate cytokine production and immune-cell activation;²¹ and muscle-derived signals may shape systemic inflammation.²² Together, these observations provide a biologically plausible link between MYL1 and inflammatory pathways relevant to RA.²³ Furthermore, the mutation of MYL1 leads to severe congenital myopathies.^24–30 The potential role of MYL1 in RA pathogenesis, including its involvement in immune activation and synovial injury, is currently unknown and has not been investigated.

The growing use of machine learning in biomarker discovery offers powerful tools for analyzing high-dimensional omics data and identifying reproducible candidates with clinical potential. In this study, we investigated differentially expressed genes (DEGs) in synovial tissue from RA patients and healthy controls. By integrating two machine learning algorithms, least absolute shrinkage and selection operator (LASSO) regression and random forest (RF), and supplementing them with cross-validation to mitigate overfitting, we systematically ranked features and prioritized genes consistently highlighted across methods and datasets. This integrative computational approach identified MYL1 as a candidate of interest. Subsequent bioinformatics analyses, including functional enrichment and immune infiltration profiling, along with correlation analyses and animal experiments, were employed to further assess its relevance to RA pathogenesis. Based on these findings, we hypothesize that MYL1 may serve as a diagnostic biomarker and therapeutic target for RA.

Materials and Methods

Data Acquisition and Processing

Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) was searched using the term “rheumatoid arthritis” and narrowed the results to “Series” and “Expression profiling by array” for Homo sapiens. A total of 100 datasets were returned, of which six were eligible and downloaded: GSE205962, GSE10500, GSE15573, GSE206848, GSE236924 and GSE77298 (Table 1). Data correction and batch effect removal were performed using the “limma” package (version 3.52.4) in R (version 4.2.1). The datasets included tissue samples from RA patients (n = 93; RA group) and healthy individuals (n = 43; control group). In addition, GSE1919, which contained synovial tissue samples from both RA patients (n = 5) and healthy individuals (n = 5), was used as a validation dataset for unbiased external validation. The results of data collation were visualized using box plots (Figure S1). The study was approved by the Clinical Research Ethics Committee of Guangxi Medical University (Approved Number: 202208035).

Table 1 Information on the GEO Dataset Used

Screening Feature Genes

The identifying of DEGs was performed as previously described.³¹ Briefly, we employed the “limma” package in R to detect DEGs within all datasets, applying a threshold of log2 fold change greater than one and an adjusted p-value under 0.05. Subsequently, to determine a selection of marker genes, the LASSO regression technique was used, leveraging the “glmnet” package (version 4.1–8) alongside the “e1071” package (version 1.7–14), family = “binomial”, standardize = TRUE, using 10-fold cross-validation (cv.glmnet, nfolds = 10). The final penalty was set to lambda.1se to favor parsimony. Additionally, through a RF analysis with the “randomForest” package (version 4.7–1.2), with ntree = 2000, mtry = floor(sqrt(p)), nodesize = 1, and permutation importance (importance = TRUE, type = 1). Importance stability was assessed by 100 bootstrap resamples, another collection of marker genes was pinpointed, focusing on those with gene importance scores surpassing 2 (rfGenes > 2).³² Last, by intersecting the marker genes identified by the two methods, a subset of candidate marker genes was identified. One muscle-related gene, which was not previously reported to be associated with RA, was selected for further analysis based on the prevalence of muscle symptoms in RA.^17–19

Final Validation and Accuracy Assessment

Violin plots were produced utilizing the “ggpubr” package (version 0.6.0) to illustrate variations in marker gene expression between the RA and the control groups, considering a p-value of less than 0.05 as statistically significant. To evaluate the marker genes’ ability to differentiate the RA samples from the controls, receiver operating characteristic (ROC) curves were plotted with the “pROC” package (version 1.18.5). The accuracy was established by computing the area under the curve (AUC), where an AUC exceeding 0.7 was deemed to indicate satisfactory accuracy.³¹ To confirm the validity of the results from the experimental group, ROC curve analysis, along with differential expression analysis, was applied to the validation dataset.

Analysis of Marker Gene Correlation

After identifying the marker genes associated with RA based on the established screening criteria, samples were categorized into high-expression and low-expression groups according to the expression levels of the marker genes. To identify DEGs strongly associated with the marker genes, we conducted a correlation analysis employing Spearman correlation coefficient (with BH correction for multiple testing). Associations meeting the thresholds of |r| > 0.6 and FDR < 0.05 were considered significant. Data visualization was performed using heatmaps and volcano plots, and gene correlation bubble plots were generated using the “corrplot” package (version 0.95).

Functional Enrichment Analysis

Through the application of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) assessments, the functional data pertinent to DEGs were mapped onto specific pathways and biological processes.³¹ This allowed for the exploration of the possible mechanisms through which these genes may contribute to RA. Additionally, we employed gene set enrichment analysis (GSEA) to investigate the enrichment of biological pathways associated with high and low marker gene expression levels. The analysis utilized a weighted statistic with 1000 permutations, and gene sets were filtered to include those containing between 10 and 500 genes.

Evaluation of Immune Cell Penetration and Associated Analysis

The examination of immune cell infiltration was conducted using the “CIBERSORT” package (version 0.1.0),³³ allowing for a quantitative assessment of the distribution of 22 lymphocyte subtypes within the samples. To illustrate the findings, stacked bar graphs and box plots were employed, with significance determined by a p-value threshold of 0.05. To explore the connections between selected marker genes and various immune cell categories, Spearman correlation tests were utilized. Significant gene-immune cell interactions were identified when the correlation coefficient |r| exceeded 0.6 and adjusted p-values were below 0.05.³¹ The results are presented graphically.

Competitive Endogenous RNA (ceRNA) Network Construction

The regulatory relationships between marker genes and microRNAs (miRNAs) were predicted using TargetScan (https://www.targetscan.org/; contextual score percentile ≥ 90),^34,35 miRanda (Director package, version 3.1.9; downloaded from http://mirtoolsgallery.tech/; top 1% genes),³⁶ and miRDB (https://mirdb.org/; target score ≥ 60).³⁷ Subsequently, the spongeScan database was used to identify long non-coding RNAs (lncRNAs) that shared microRNA response elements (MREs) and could competitively bind to the selected miRNAs.³⁸ The network was visualized using Cytoscape v3.10.1, an open-source platform for building complex networks.³⁹

Establishment of the Rat RA Model and Collection of Tissues

Sample size was determined a priori using G*Power 3.1 for a two-tailed t-test (α = 0.05, power = 0.80). Based on previous studies in collagen-induced arthritis rat models showing large gene expression differences between RA and control groups,⁴⁰ a mean difference of 1.0 and a standard deviation of 0.6 (Cohen’s d = 1.67) were assumed. The estimated sample size was 8 rats per group; to ensure adequate power, 15 rats were included per group. Therefore, A total of 30 male Sprague-Dawley rats, aged 8 weeks and weighing between 220 and 250 grams, were purchased from the Experimental Animal Center of Guangxi Medical University. All rats were housed in a specific pathogen-free (SPF) environment and randomly divided into the control and the RA groups (n =15 for each group). After one week of acclimatization, the rats received a subcutaneous injection of 0.2 mL of bovine type II collagen (2 mg/mL, Chondrex, USA) 1:1 emulsified with an equal volume of complete Freund’s adjuvant (CFA, Sigma-Aldrich, containing 5 mg/mL Mycobacterium tuberculosis) at the base of the tail, to establish the collagen-induced arthritis (CIA) rat model.⁴¹ The control group received an equal volume of saline by subcutaneous injection. Seven days later, a booster immunization was administered by subcutaneous injection, also at the tail, using a 1 mL Luer-lock syringe. The injection consisted of 0.05 mL of a mixture emulsified 1:1 with bovine type II collagen (1 mg/mL) and incomplete Freund’s adjuvant (IFA, purchased from Sigma-Aldrich).

The development of arthritis was monitored every 4 days after the initial immunization. The severity of arthritis in each paw was quantified using a standard arthritis index scoring system (0–4 points per paw; maximum score = 16 points per rat). 0 = no visible signs of inflammation; 1 = swelling confined to one digit; 2 = swelling involving two or more digits; 3 = moderate redness and swelling of all digits and the paw; 4 = severe swelling and deformation involving the paw and ankle joints.⁴²

At day 28 postoperatively, all rats underwent euthanasia by intraperitoneal administration of sodium pentobarbital at 150 mg/kg. Synovial tissue was collected and immediately stored at −80°C for subsequent analysis by a blinded operator. To validate the systemic inflammatory response, we performed qRT-PCR to measure the mRNA expression levels of TNF-α and IL-6, using GAPDH as an internal reference.

RNA Extraction and qRT-PCR

Total RNA was extracted using the RNAeasy™ Plus Animal RNA Extraction Kit (Beyotime, Shanghai).³¹ RNA concentration was determined with a Nanodrop spectrophotometer (Thermo Fisher, USA). One microgram of total RNA was reverse transcribed into cDNA using the HyperScript™ RT SuperMix for qPCR (APExBIO, USA). qRT-PCR was performed with 2X SYBR Green qPCR Master Mix (APExBIO, USA), and primers were synthesized (Nanning GenSys Biotechnology Company Limited, China). GAPDH was used as the internal control. Relative mRNA expression levels were calculated using the 2^-ΔΔCt method. All experiments were performed in triplicate. Primer sequences are listed in Table 2.

Table 2 Primer Sequences for PCR

Western Blot

Total protein was isolated from synovial cells utilizing RIPA lysis buffer supplemented with Serine protease inhibitor (BOSTER, Wuhan, China), phosphatase inhibitor cocktail (CoWin Biotech, Jiangsu, China), and protease inhibitor cocktail (APExBIO, USA).^31,43 The samples were mixed with loading buffer, the proteins were denatured in a metal bath at 100°C for 10 minutes and electrophoresed on a 10% polyacrylamide gel. After transfering onto a PVDF membrane, the proteins were blocked at room temperature for 2 hours and then incubated overnight with primary antibodies against MYL1 (Proteintech, Wuhan, China; dilution 1:1500) and GAPDH (Proteintech, Wuhan, China; dilution 1:5000). The membrane was then incubated for 1 hour with a goat anti-rabbit IgG secondary antibody (Zen Bioscience, Chengdu, China), and visualization was performed using chemiluminescence.

Statistical Analysis

The Shapiro–Wilk test was used to evaluate the normal distribution of qRT-PCR and Western blotting data. Independent sample t-tests were utilized for intergroup comparisons of variables showing normal distribution. The data are displayed as the average ± standard deviation. A p-value below 0.05 was considered statistically significant. Statistical analyses were executed using GraphPad Prism 9.5 (GraphPad Software, USA).

Results

Identification of DEGs

The detailed research flowchart is shown in Figure 1. As shown in Figure 2A and B, principal component analysis (PCA) demonstrated that batch effects were largely eliminated following batch correction and data normalization, as samples from different datasets became more intermixed in the PCA plot. A total of 34 DEGs were identified, including 28 that were downregulated and 6 that were upregulated. These DEGs are illustrated in Figures 2C and D, highlighting distinct expression patterns between the groups.

Figure 1 The flow chart. Abbreviations: DEG, differentially expressed gene; GEO, Gene Expression Omnibus; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; WB, Western blotting; PCR, polymerase chain reaction; GSEA, gene set enrichment analysis; ceRNA, competing endogenous RNA.

Figure 2 Identification of DEGs. (A and B) Batch pre-correction and post-correction for the GSE205962, GSE10500, GSE15573, GSE206848, GSE236924, and GSE77298 datasets. (C) Heatmap of DEGs, with highly expressed DEGs in the samples marked in red and low-expressed DEGs in the samples marked in blue, the red horizontal line indicates MYL1. (D) Volcano plots of DEGs, with red indicating up-regulated genes, blue indicating down-regulated genes, and gray indicating non-significant genes, MYL1 is marked by the red box. Abbreviation: DEG, differentially expressed gene.

Screening Using a Combination of Two Machine Algorithms

As shown in Figure 3A, cross‑validation of the LASSO model determined the optimal penalty parameter (λ), minimizing the binomial deviance and yielding a set of feature genes. Figure 3B illustrates the coefficient profiles of the candidate genes as the λ value changed. Based on this analysis, seven genes were selected, including MYL1, IGJ, MMP1, CDH11, TNFAIP6, CXCL13, and RPS4Y1. Figure 3C presents the error rate of the RF model with an increasing number of decision trees, showing stable performance of the classifier. Figure 3D ranks the top genes according to their importance scores in the RF algorithm, among which MYL1, IGJ, MMP1, CDH11, and TNFAIP6 exhibited particularly high importance values. Finally, the Venn diagram (Figure 3E) shows that LASSO and random forest analyses jointly identified five overlapping genes: TNFAIP6, MYL1, CDH11, MMP1, and IGJ. These genes were therefore considered as potential disease‑related feature genes and candidates for further investigation in RA.

Figure 3 Identification of the MYL1. (A and B) LASSO regression counting differential genes, regression graphs, and cross-validation graphs, respectively, screening 13 differential genes. (C) Displays the cross-validation error curves of the random forest model. The x-axis represents the number of trees, and the y-axis shows the error rate. This plot evaluates the model’s stability and predictive performance to identify the most important differential genes. (D) Shows the importance scores of 10 differential genes based on the random forest model. The x-axis indicates the importance score, while the y-axis lists the gene names, MYL1 is marked by the red box. (E) Venn diagram of feature genes screened by the two machine learning methods to screen feature genes in a Venn diagram. Finally, five genes (ie, MYL1, IGJ, MMP1, CDH11, and TNFAIP6) were screened out. Abbreviations: DEG, differentially expressed gene, the red horizontal line indicates MYL1.

Intergroup Differences and Accuracy of Signature Genes

Violin plots (Figure 4A–C, G and H) alongside ROC curves (Figure 4D–F, J and K) were utilized to illustrate the contrasting expression levels of the MYL1, IGJ, CDH11, TNFAIP6 and MMP1 genes in the training data. Findings showed decreased MYL1 expression and increased IGJ, CDH11, TNFAIP6, and MMP1 expression in the RA group compared with controls. The ROC curve’s area for the MYL1, CDH11, TNFAIP6 and MMP1 genes surpassed a threshold of 0.7.

Figure 4 Differential expression and ROC curves of the five genes. (A–C, G and H) The differential expression of MYL1, IGJ, CDH11, TNFAIP6 and MMP1 genes in the control and RA groups. ** p<0.01, *** p < 0.001. (D–F, J and K) The ROC curves and areas under the curves of MYL1, IGJ, CDH11, TNFAIP6 and MMP1 genes. (I) The violin plot of MYL1 differential expression in the GSE1919 dataset. (L) The ROC curve and area under the curve for MYL1 in the GSE1919 dataset. * p < 0.05. Abbreviations: ROC, receiver operating characteristic; AUC, area under the curve; 95CI, 95% confidence interval; RA, Rheumatoid arthritis.

MYL1 Expression Was Reduced in the Validation Group Dataset

Upon validation of the MYL1, IGJ, CDH11, TNFAIP6 and MMP1 genes in the GSE1919 dataset, it was found that only MYL1 displayed differential expression, aligning with previous findings (Figure 4I). In the GSE1919 dataset, the ROC curve for MYL1 demonstrated an area under the curve exceeding 0.7 (Figure 4L), reinforcing its status as the ultimately identified differential gene.

Expression of MYL1 in the Synovium of RA Rats

It was found that in both humans and rats, the MYL1 gene (myosin, light chain 1) is abbreviated as MYL1. Subsequently, these genes’ mRNA and protein sequences were retrieved from the NCBI database (https://www.ncbi.nlm.nih.gov/). Sequence comparison was performed and revealed that these two genes had high similarity in mRNA (77.93% similarity) and protein (92.78% similarity) sequences (Figure S2). Representative images of CIA model arthritis score at different levels (0–4) and the progression of arthritis scores during 28 days post‑immunization are shown in Figure 5A and B. Successful establishment of the CIA model was confirmed when rats exhibited a progressive increase in arthritis scores, accompanied by visible swelling and redness of the plantar and toe joints. At the molecular level, qRT‑PCR analysis confirmed markedly increased expression of the proinflammatory cytokines TNF‑α and IL‑6 in the synovial tissue of the non-injected paws of CIA rats compared with the control group (both P < 0.05; Figure 5C and D). Finally, qRT-PCR analysis demonstrated a pronounced reduction in MYL1 mRNA expression in RA samples relative to the control group (p < 0.001; Figure 5E). Similarly, Western blot showed that the protein expression level of MYL1 was significantly reduced in the RA group (p = 0.0093; Figure 5F and G, Figure S3). Consistent evidence from bioinformatics analyses and in vivo experiments in rats further demonstrated that the expression of MYL1 was significantly lower in the RA group compared to the control group.

Figure 5 Validation of the CIA rat model and expression of MYL1 in synovium of the two groups of rats. (A) Representative images showing arthritis severity scoring from level 0 to 4. (B) Arthritis scores over 28 days in RA and control groups. (C and D) The mRNA expression of TNF-α and IL-6 in the synovium of RA and control rats using qRT-PCR. (E) The mRNA expression of MYL1 in the synovium of RA and control rats using qRT-PCR. (F) Western blot analysis of MYL1 protein expression with GAPDH as loading control. (G) Quantification of MYL1 protein levels normalized to GAPDH. *p < 0.05; **p < 0.01; *** p < 0.001. Abbreviations: CIA, collagen‑induced arthritis; RA, rheumatoid arthritis.

Related Genes of MYL1 in RA

After the identification of DEGs, the expression matrix of MYL1 was analyzed by Spearman correlation, and its related genes were demonstrated by a heatmap (Figure 6A). And a visual bubble map was generated to show the correlation between DEGs associated with MYL1 (Figure 6B).

Figure 6 Differential analysis of MYL1 gene groupings. (A) The MYL1 differential analysis heatmap; (B) The MYL1 differential analysis correlation graph.

RA Is Associated with Multiple Functional States of Synovial Cells and Disease Pathways

According to the results of GO analysis, the DEGs were mainly enriched in muscle-related biological processes, and these pathways were closely related to functions such as key roles in muscle structure and function, particularly in the muscle structure development, muscle contraction, the organization of muscle fibers, and the assembly of RA-associated cellular components (Figure 7A–D). The KEGG assessment indicated that the RA group was enriched in various pathways linked to muscle and heart conditions, including those related to calcium signaling, motor protein signaling, myocardial contraction, and adrenergic signaling within cardiomyocytes (Figure 8A–D).

Figure 7 GO analysis. (A) The GO analysis bar plot. (B) The GO analysis bubble plot. (C) The GO analysis chord chart. (D) The GO analysis circle plot. Abbreviation: GO, Gene Ontology.

Figure 8 KEGG analysis. (A) The KEGG analysis bar plot. (B) The KEGG analysis bubble plot. (C) The KEGG analysis chord chart. (D) The KEGG analysis circle plot. Abbreviation: KEGG, Kyoto encyclopedia of genes and genomes.

Gene Set Enrichment Analysis

Utilizing GSEA analysis to explore individual genes linked to MYL1 offered further insights into how MYL1 can differentiate RA samples from controls, highlighting the pathways regulated by MYL1. Through a comprehensive examination, it was discovered that the group with elevated expression levels predominantly displayed enrichment in five biological signaling pathways. These included pathways related to antigen processing and presentation, chemokine signaling, and lysosomal functions (Figure 9A). Conversely, the group exhibiting lower expression levels demonstrated a significant concentration in five activities and pathways linked to muscle operation and heart issues, such as the calcium signaling pathway and myocardial contraction processes (Figure 9B).

Figure 9 Gene set enrichment analysis. (A) Functional enrichment analysis of GSEA in the high-MYL1 expression group. (B) Functional enrichment analysis of GSEA in the low-MYL1 expression group. Abbreviation: GSEA, Gene set enrichment analysis.

MYL1 Is Associated with Immune Cells and Their Functions in RA

Using the CIBERSORT algorithm, we analyzed the differences in peripheral blood immune cells between the RA and the control groups. The proportional stacked bar graph (Figure 10A) shows the proportional expression of 22 immune cell types across all samples. The first 43 cases were in the control group, and the last 93 cases were in the RA group. The expression levels of nine immune cell markers exhibited notable variations (Figure 10B and C). Specifically, in the RA group, there was an elevation in Macrophages M0 and M1, activated CD4 memory T cells, follicular helper T cells, and gamma delta T cells. Conversely, activated NK cells, resting CD4 memory T cells, CD8 T cells, and regulatory T cells (Tregs) were found to be lower in RA patients compared to the control group. Subsequently, we conducted an analysis to examine the correlation between MYL1 expression and various immune cell types. The heatmap, box plot, and lollipop plot visualizations (Figure 11A–C) illustrate that MYL1 showed a positive correlation with Tregs, activated NK cells, plasma cells, and naïve CD4 T cells, while demonstrating a negative correlation with M2 macrophages (Figure 11D).

Figure 10 Immune cell infiltration analysis. (A) A proportional stacked bar chart showing the distribution of immune cells in the control and RA groups. (B) The split violin plot shows the differences in immune cell distribution between the control group and the RA group. (C) The immune abundance box plot shows the differential distribution of immune cell abundance for immune cells with significant differences between the control and RA groups. *p < 0.05; ** p < 0.01. Abbreviation: RA, Rheumatoid arthritis.

Figure 11 Immune correlation analysis. (A) Heatmap of the correlation between MYL1 and 22 immune cell types. (B) Box plot of the correlation between MYL1 and 22 immune cells. (C) Lollipop plot of correlation analysis between MYL1 and 22 immune cells, the red boxes indicate the immune cell types with statistically significant correlations in the Spearman correlation tests. (D) The expression level of MYL1 showed a positive correlation with T cells regulatory (Tregs), NK cells activated, Plasma cells, T cells CD4 naive, but a negative correlation with Macrophages M2. (p<0.05). *p < 0.05; ** p < 0.01.

ceRNA Network

As illustrated in Figure 12, the resulting network consists of 118 nodes, including 1 signature mRNA (MYL1), 22 miRNAs (eg, hsa-miR-4455, hsa-miR-3919, hsa-miR-4513), and 92 lncRNAs (such as HAGLR, SNHG16, and NEFL). The network comprises a total of 130 edges, representing both miRNA–mRNA and lncRNA–miRNA interactions. It was shown that MYL1 regulates the expression of itself and other molecules through a complex ceRNA network and may play an important role in the pathogenesis of RA.

Figure 12 A ceRNA networks of MYL1 and the potential RNA regulatory pathways. Abbreviation: ceRNA, competitive endogenous RNA.

Discussion

In this study, through the rigorous screening of two machine learning methods and inflammatory processes, we identified 34 DEGs significantly associated with RA. Although MYL1 did not exhibit the highest fold change among all differentially expressed genes, its biological significance in the context of RA should not be underestimated. Muscle weakness and atrophy are common features of RA,^17,18 and even modest alterations in genes regulating muscle function, such as MYL1, could have a meaningful physiological impact in RA patients. In addition, the expression changes of MYL1 were verified by constructing a rat RA model, which further supported the potential role of MYL1 in RA.

Notably, the expression of MYL1 is upregulated in several tumors (rhabdomyosarcoma, head and neck squamous cell carcinoma, breast cancer, and prostate cancer).^24–26,44 Furthermore, MYL1 plays an important role in several musculoskeletal diseases.^27,45,46 These studies suggest that MYL1 may be broadly involved in pathological processes beyond muscle disease and reinforce its potential value as a disease-related biomarker. This investigation identified novel molecular targets for detection through an analysis of immune infiltration, discovering that MYL1 is linked to modified immune activity. This aligns with earlier research suggesting a strong connection between RA progression and immune function.³

Based on the results of the GO analysis, we discovered that MYL1 mainly participates in organizing sarcomeres, facilitating muscle contraction, and assembling RA-associated cellular components. This finding is in line with that of a prior study, which reported that one of the clinical manifestations of RA is muscle weakness and atrophy.^17,18 Alterations in muscle function may exacerbate the clinical symptoms of RA by affecting the structure and function of muscle tissue, especially impairments in motor function and quality of life. Over a period of five years, a study involving 288 individuals diagnosed with RA found that although disease activity and physical capabilities improved, low muscle mass persisted unchanged. This indicates that muscle depletion in those with RA could independently worsen functional impairment, regardless of inflammation levels.⁴⁷ This was further confirmed by the KEGG analysis, which showed that the RA group was significantly enriched in calcium signaling pathways, motor protein signaling pathways, cardiac contraction signaling pathways, and other pathways in muscle and some heart diseases. These results suggest that RA not only affects the degeneration of articular cartilage, but may also lead to muscle weakness and other cardiac-related complications in patients by affecting muscle and cardiac function, which also coincides with the results of previous studies.^48–54 An important association exists between RA and both heart failure and impaired myocardial function. Studies have demonstrated that in RA, decreased expression of the transcriptional repressor IKAROS in CD4 T cells results in upregulation of ORAI3, thereby activating calcium-selective (ARC) channels. This process enhances the sensitivity of T cells to inflammatory stimuli, and aberrant T-cell activation is considered a key contributor to the development of chronic inflammation in RA.⁵⁵ Other studies have demonstrated enhanced calcium signaling in synovial fibroblasts in RA, with translocation of nuclear factor c3 (NFATc3) from activated T cells to the nucleus, resulting in increased expression of the inflammatory factor RANTES, which ultimately leads to an enhanced synovial inflammatory response in RA.⁵⁶ In addition, enhanced calcium signaling is involved in the anti-apoptotic process of synovial fibroblasts in RA.⁵⁷ Furthermore, Liang et al comprehensively reviewed the critical role of calcium signaling pathways in the pathogenesis of RA.⁵⁸ Downregulation of MYL1 in RA is associated with weakened immunosuppressive effects. This may involve the calcium signaling pathway and macrophage polarization.

This study examined immune cell infiltration to explore potential immunological associations of MYL1 expression in rheumatoid arthritis (RA). The results suggest that MYL1 expression is correlated with altered immune cell proportions, rather than directly demonstrating an immunoregulatory function. This correlation supports previous evidence that immune dysregulation is closely involved in RA progression,^59–61 and Alivernini suggested that T cells may play a role in early RA.⁶² In contrast, elevated macrophage levels in the RA group indicate an ongoing inflammatory response within the synovium,^63–65 and it has been demonstrated that suppressing the activation of inflammatory vesicles in macrophages can effectively reduce RA-related tissue damage.⁶⁶ The reduction of NK cells activated and Tregs implies a weakening of the immune-suppressing function of RA, which further amplifies the autoimmune response and inflammatory environment. Recent research indicates a strong association between lower NK cell counts and heightened RA symptoms.⁶⁷ Moreover, Treg cells have gained attention as a potential therapeutic focus for RA treatment.⁶⁸ In RA patients, an almost complete lack of GPSM2 in Treg cells impairs their migration to inflammatory sites.⁶⁹ Furthermore, enhancing the Treg/Th17 ratio by increasing Treg cell proportions can markedly improve the course of RA.⁷⁰

In this research, the analysis of immune cell associations revealed a strong positive correlation between T regulatory cells (Tregs) and activated natural killer (NK) cells with MYL1. Additionally, M1 macrophages and follicular helper T cells also exhibited a significant correlation with MYL1. Down-regulation of MYL1 may lead to a decrease in immunosuppressive cells (Tregs and NK cells activated) and an expansion of pro-inflammatory macrophages (M1 macrophages) and T follicular helper cells. Downregulation of MYL1 may lead to a decrease in immunosuppressive cells (Tregs and NK cells activated) and an expansion of pro-inflammatory macrophages (M1 macrophages) and T follicular helper cells. Previous studies have demonstrated that FABP4 secreted by M1-polarized macrophages promotes synovitis, angiogenesis, and cartilage degeneration, exacerbating the progression of RA.⁷¹ Secondly, a significant expansion of T follicular helper cells was found in early RA synovial tissues.⁷² Follicular helper T cells facilitate the differentiation of plasma cells by releasing IL-21, moving towards the inflammation site under the influence of chemokines, and enhancing B-cell activity and antibody generation in non-lymphoid tissues experiencing pathological inflammation.⁷³ One of the pathological features of RA is long-term chronic synovial inflammation and dysregulation of autoimmunity, and this present study, along with related studies, demonstrated that dysregulation of Tregs, NK cells activated, M1-type macrophages, and T cells follicular helper correlates with disease progression and severity.

Currently, there is a scarcity of research focused on how MYL1 influences the functions of immune cells. Nonetheless, MYL1 encodes myosin light chain, a structural component of the cytoskeleton that contributes to cellular contractility and signal transduction. Increasing evidence indicates that myosin light‑chain signaling plays crucial roles in immune modulation. For example, inhibition of myosin light‑chain kinase (MLCK) has been shown to attenuate inflammatory responses and reduce cytokine release.^74,75 Moreover, MLCK activity regulates T‑cell motility and polarization,⁷⁶ as well as neutrophil‑mediated inflammatory injury.⁷⁷ In rheumatoid arthritis, the HIF‑1α‑RhoA/ROCK axis, which converges on myosin light‑chain phosphorylation, promotes fibroblast‑like synoviocyte migration and joint inflammation.⁷⁸ Therefore, dysregulation of MYL1 expression might influence cytoskeletal remodeling and immune‑cell activity, providing a possible mechanism by which MYL1 contributes to RA pathogenesis. In addition, investigations have revealed that MYL9, another member of the myosin light chain family, can act as a ligand for CD69.⁷⁹ The interaction facilitated by MYL9 between cancer-associated fibroblasts (CAFs) and immune cells plays a role not only in the development of tumors and the evasion of the immune system but also in controlling the movement and activities of various immune cells, including T cells, macrophages, and NK cells, within tissues.⁸⁰ Moreover, in colorectal cancer, high MYL9 expression is associated with poor response to immunotherapy.⁸¹ Although there have been no direct reports on the expression or specific functions of MYL1 in immune cells, considering its structural and functional homology with MYL9, its role in muscle cytoskeletal regulation, and our findings that MYL1 expression changes are closely associated with immune infiltration and abnormal immune cell proportions in RA, it is plausible that MYL1 might potentially contribute to the regulation of immune cell cytoskeletal remodeling and migration through analogous to those of MYL9. However, this assumption requires further experimental validation to confirm whether MYL1 indeed plays such a role in RA.

We developed a ceRNA network focused on MYL1, where lncRNAs can rival mRNAs for miRNA binding sites, influencing MYL1 expression. Discovering miRNAs that directly engage with mRNA and lncRNAs that indirectly affect mRNA levels offers novel avenues for RA prevention or therapy.

The research existed several limitations. To begin with, the human data were sourced from the GEO database, concentrating on cartilage samples from individuals with RA and those without the condition. Nonetheless, the sample size was somewhat limited, and the absence of a predefined chronological framework for data integration may introduce additional heterogeneity. Additionally, our findings were corroborated using a rat model of RA, but have not yet been confirmed in RA patients. Consequently, further clinical research is essential to validate the roles of MYL1 in the diagnosis and treatment of RA. Future comparative clinical investigations are necessary to substantiate these results. Third, although we employed cross-validation and multiple algorithms to enhance robustness, machine learning approaches inherently risk overfitting, particularly with limited sample sizes. Moreover, CIBERSORT-based estimates are limited by dependence on the reference signature, model assumptions. And ceRNA links are computational and hypothesis‑generating, causal inference awaits reporter, AGO2‑RIP/pull‑down, and perturbation/rescue validation in RA‑relevant models. Finally, while we established an association between MYL1 and RA, the pathophysiological mechanisms underlying the role of MYL1 in RA remain incompletely understood. Further research is needed to elucidate the molecular and signaling mechanisms of MYL1 in RA.

The identification of MYL1 as a candidate biomarker could augment early diagnosis of RA, enable patient stratification, and facilitate treatment monitoring in conjunction with existing markers. Progress toward clinical application will require stepwise analytical and clinical validation, assay standardization, and deliberate integration into diagnostic pathways. Priority directions are to validate performance across multicenter blood and synovial cohorts and establishing clinically actionable cut‑offs linked to disease activity, joint damage progression, and therapeutic response. In parallel, a deployable ELISA should be developed, and the incremental value of MYL1 beyond RF, anti‑CCP, ESR, and CRP should be quantified using head‑to‑head comparisons together with reclassification and decision‑curve analyses. Longitudinal cohorts are needed to assess prognostic utility and the prediction of treatment response. Single‑cell and spatial profiling can localize the cellular sources of MYL1‑related signals. Finally, perturbation studies in fibroblast‑like synoviocytes, immune cells, and in vivo models should interrogate links to calcium signaling, cytoskeletal dynamics, and cytokine networks.

In conclusion, this study identifies a significant downregulation of MYL1 in RA, which is closely associated with key immune cell subsets, suggesting it may play an important role in RA pathogenesis and warrants further investigation as a potential future diagnostic marker.