A total of 12 RA patients and five healthy controls (HC) were recruited for this study from the Second Hospital of Shanxi Medical University in August 2021. All the patients met the 2010 American College of Rheumatology (ACR)/European League Against Rheumatism (EULAR) classification criteria for RA (17). To evaluate the relationship between m6A and response to treatment, genome-scale data about infliximab therapy response were also recorded in our study, which included patients’ response to anti-TNF evaluated by Disease Activity Score-28(DAS28).

Peripheral blood samples (5 ml) were collected from each patient and control subject into EDTA-2 K-containing tubes. According to the manufacturer’s protocol, fresh PBMCs from each donor blood were isolated by Ficoll-Hypaque (Beijing Solarbio Science & Technology Co., Ltd.) density-gradient centrifugation for 20 min at room temperature. Total RNA was isolated from freshly obtained PBMCs. Before RNA sequencing, the quality of RNA was assessed on the Agilent Bioanalyzer 2100 system. A 1.5-μg RNA sample was taken for RNA sequencing. The library preparations were sequenced on Illumina. Removing reads containing adapter, containing ploy-N, and low-quality reads containing >50% bases with qualities of ≤20 from raw data, clean reads were obtained.

Reference genome and gene model annotation files were downloaded directly from the Genome website, and paired clean reads were aligned to the reference genome using Hisat2 v2.0.5. Reads mapped to each gene were calculated using featureCounts v1.5.0-p3. Then the FPKM for each gene was calculated.

Samples from the GEO cohort were randomly divided into the training set (935 samples, 70%) and the testing set I (400 samples, 30%). During training, to avoid the overfitting problem caused by random oversampling, the SMOTE sampling method was adopted to repeatedly sample the healthy controls, which had fewer samples, and thereby balance the number of samples in the RA and HCs. The R package “pROC” was used to evaluate the visualization of the receiver operating characteristic curve (ROC) to calculate the area under the curve (AUC). The sequencing data were employed for external validation.

Twenty-one acknowledged m6A regulator genes were referred (10, 18, 19). Only 20 regulators were stably expressed and used to identify distinct m6A methylation modification patterns. These 20 m6A regulators included eight methyltransferases (METTL3, METTL14, RBM15, RBM15B, WTAP, VIRMA/KIAA1429, CBLL1, ZC3H13), two demethylases (ALKBH5, FTO), and 10 RNA-binding proteins (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPA2B1, HNRNPC, FMR1, LRPPRC, ELAVL1). Unsupervised clustering analysis was applied to identify different m6A modification patterns based on these 20 m6A regulators by the “ConsensusClusterPlus” package. The optimal and stable numbers of clusters were selected according to cophenetic, dispersion, and silhouette coefficients.

To characterize the biological features between different m6A modification patterns, we utilized single-sample gene set enrichment analysis (ssGSEA) to estimate the population of specific infiltrating immunocytes and the activity of immune reactions. The gene sets marking each infiltrating immunocyte type were obtained from the previous study17, and the RA-related pathways were obtained from the MSigDB database. The enrichment scores defined by ssGSEA analysis represent the degree of each immunocyte abundance and immune reaction activity in each sample, compared to three distinct modification patterns by the Wilcoxon test.

M6A-related differentially expressed genes (DEGs) in different m6A phenotypes were identified by the “Limma” R package, in which P-values< 0.05 were set as cutoff criteria for the DEGs. GO and KEGG enrichment analyses were applied to analyze the biological significance of DEGs. P< 0.05 was considered statistically significant, and visualization of results was conducted by R package ‘ggplot2’ and ‘GOplot’.

To quantify m6A modification patterns in individual RA patients, we constructed a set of the m6A gene signature (m6Ascore) by using principal component analysis (PCA) algorithms. The differentially expressed genes (DEGs) obtained from three clusters in the previous step were intersected to get shared DEGs between three m6A clusters. Pearson correlations were performed to obtain positive or negative correlation signature genes for co-expressed DEGs. The signature genes were further selected for features by the “Boruta” R package. We then conducted PCA based on the final determining genes. Principal components from the two groups of signature genes were separately extracted and served as the final signature score by subtraction. This approach concentrates the score on the set with highly correlated or anticorrelated gene blocks while down-weighting the gene contributions that are not tracked with other set members. We then define the m6Ascore by adopting a formula like previous studies (18, 20, 21).

where i is the signature score of clusters that have a positive coefficient, and j is the expression of genes that have a negative coefficient.

A total of 1,203 RA and 149 healthy controls (HCs) from six cohorts were included in this study ( Supplementary Table 1 ). PCA was used to visualize variation in correcting for batch effects ( Supplementary Figure 1 ). We noticed a very close association among methyltransferases in the 20-m6A regulator protein–protein interaction (PPI) network, which usually function as a complex ( Supplementary Figure 2 ). Subsequently, we explored the different expressions of 20 m6A regulators between RA and healthy controls in a GEO public cohort ( Figures 1A–C ). Five regulators were observed to have consistently different expressions in 12 RA patients and five HC participants from the Second Hospital of Shanxi Medical University ward and health examination center again ( Figure 1D ): FMR1, HNRNPC, LRPPRC, WTAP, YTHDF3. To examine the ability of the 20 m6A regulators to distinguish between RA and HCs, we conducted a diagnosis model by the random forest (RF) algorithm. The area under the curve (AUC) value was 0.83 ( Figure 1E ). External validation was performed for the diagnostic efficacy of the m6A-RF model ( Figure 1F ). Actually, this diagnostic model is not highly efficient in the validation cohort because of the insufficient sample size, but it also partly explains the role of m6A in the occurrence and development of RA.

To investigate transcriptome relationships, we calculated pairwise correlations among the expressions of the 20 m6A regulators. We found that positive correlations were more frequent than negative correlations ( Figure 2A ). Next, based on the expression profiles of the 20 selected m6A regulators, we utilized consensus clustering analysis to stratify patients to different m6A modification patterns ( Figures 2B–D ). Accordingly, we determined that the matrix heatmap retained sharp and clear sides when k = 3, which indicated there were three distinct m6A modification pattern clusters, including 279 cases in cluster A, 581 cases in cluster B, and 331 cases in cluster C ( Figures 2E, F ). We termed these clusters as m6A cluster A, m6A cluster B, and m6A cluster C.

To identify the immune microenvironment characteristics underlying three distinct m6A modification patterns, we compared the enrichment scores of RA-related pathways and immune cell infiltration ( Supplementary Tables 2, 3 ) among the RNA modification patterns. We found inflammatory cell infiltrates, including neutrophils, monocytes, and T helper type 17 (Th17) in cluster A, with the enormous imbalance between T helper type 1 and T helper type 1 (Th1/Th2) at the same time ( Figure 3A ). Patients in cluster C were rich in adaptive immune-related cells such as activated CD4+T cells, activated CD8+T cells, and activated B cells ( Figure 3A ). In cluster C, Th2 and eosinophil abundant infiltration was considered a protective factor by counteracting the development of arthritis and preventing bone loss20. Cluster B was modestly activated in most inflammatory and immune cells. Nevertheless, cluster B displayed more activation of natural killer cells (NK), natural killer T cells (NKT), gamma delta T cells (γδT), and CD56 bright natural killer cells (CD56bright NK), which are typical innate lymphoid cells. In parallel with this, the RA-related biological process was differentially activated in the three subgroups ( Figure 3B ). Cluster A showed strong enrichment for most inflammatory pathways, including acute and chronic inflammation, response to bacterium and virus, complement activation, and chemokine. Therefore, cluster A can be referred to as the highest inflammatory phenotype. In contrast, the abovementioned inflammatory pathways were remarkably less expressed in cluster C than in clusters A and B, and we assume cluster C as an adaptive lymphocyte-rich phenotype. Cluster B was enriched in natural killer cell-mediated cytotoxicity, IL-17 signaling pathway, JAK-STAT signaling pathway, and TNF signaling pathway. They were identified as an innate lymphocyte-rich phenotype.

To determine whether the molecular cluster has an association with clinical features, we investigated the distribution of the three clusters according to disease activity ( Figure 4 and Supplementary Table 4 ). Disease Activity Score (DAS)28 - erythrocyte sedimentation rate (DAS28-ESR) and DAS28-C-reactive protein (DAS28-CRP) larger than 5.1 were regarded as high disease activity. All three clusters existed independently of disease status. Therefore, the three clusters were similar and clinically indistinguishable by DAS28-ESR/DAS28-CRP parameters.

Although RA patients were classified into three m6A modification phenotypes based on 20 m6A regulators, the underlying genetic changes and expression perturbations within these phenotypes remain unclear. To further examine potential m6A-related transcriptional expression changes in three patterns, we identified 209 overlapping DEGs using an empirical Bayesian algorithm and performed an enrichment analysis ( Figure 5A ). GO and KEGG pathway enrichment analyses were conducted to explore the functional characteristics of the DEGs ( Figures 5B, C ). The GO and KEGG analysis results showed that the DEGs were significantly enriched in “defense response to virus”, “response to interferon-beta”, “response to interferon-alpha”, “regulation of innate immune response”, “NOD-like receptor signaling pathway”, “Coronavirus disease-COVID-19”, and so on. We performed an unsupervised consensus clustering analysis based on the 209 RNA phenotype-related DEGs to further validate this differential regulation. Consistent with the clustering grouping of m6A modification patterns, we obtained three stable transcriptomic phenotypes named m6A gene clusters A–C ( Supplementary Figure 3 ), consistent with the expected results.

We applied the m6A score model to accurately evaluate the m6A modification pattern of individual patients with RA. Patients were divided into high or low m6A score groups using the median m6A score as the cutoff. In the two m6Ascore groups, almost all HLA gene expressions showed prominent differences, and that in m6Ascore-low was significantly higher than that of m6A score-high ( Figure 6A ). The m6Ascore-high group observed a significantly lower expression of 12 critical immune checkpoints than that in the m6Ascore-low group, like tumor necrosis factor superfamily15, tumor necrosis factor receptor superfamily14, CTLA4, and CD86 ( Figure 6B ). The above analysis indicated that m6Ascore might be closely connected with immunotherapy. To further understand the effects of m6Ascore on predicting drug response, we selected two independent groups, 154 RA patients treated with infliximab ( Supplementary Table 5 ) and 92 RA patients treated with rituximab ( Supplementary Table 6 ). We tested the differences in m6Ascore between infliximab responders and non-responders and found that responders had a lower m6Ascore, while patients showed a high m6Ascore with poor clinical efficacy of infliximab therapy ( Figure 6C ). The above analysis was in accordance with the expected results that the lower m6A score group, which has higher immune checkpoints TNFSF15 and TNFRSF14, may have a better anti-TNF (like infliximab) treatment effect. In cluster C, a higher m6A score was observed, which suggested that there may be a lower infliximab therapeutic response. Clusters A and B have lower m6A scores, and there may be higher infliximab therapeutic responses ( Figures 6D, E ). Nevertheless, the rituximab cohort study did not observe significant treatment differences between the two m6Ascore groups ( Supplementary Figure 4 ).