Thirty-nine premenopausal RA women regularly followed in the RA Outpatient Clinic at the Hospital das Clinicas da Universidade de São Paulo, were recruited for this study. All patients fulfilled the American College of Rheumatology/European Alliance of Associations for Rheumatology 2010 classification criteria for RA. Patients were excluded if they presented metabolic bone disease (e.g., rickets, primary hypoparathyroidism, osteomalacia, and Paget’s disease); ii) bone-associated chronic disease (e.g., renal or liver failure, hyperthyroidism, hypothyroidism, and malabsorption); iii) use of any medication interfering with bone metabolism (e.g., prednisone doses higher than 7.5 mg/day, bisphosphonates, and bone-targeted monoclonal antibodies iv) an autoimmune disease other than RA, or v) being pregnant or lactating.

Premenopausal status was assessed through a self-reported information. Participants were asked to complete a questionnaire that included questions about their menstrual history, such as the age of menarche, regularity of menstrual cycles, and the presence or absence of menstrual periods in the past 12 months.

Based on data on bone erosion from high-resolution peripheral quantitative computed tomography (HR-pQCT), RA patients were classified into erosive and non-erosive patients according to the presence or absence of bone erosions on metacarpophalangeal (MCP) and proximal interphalangeal (PIP) joints of the dominant hand.

In addition, 20 age- and body mass index-matched healthy subjects were enrolled to the control group. The study was approved by the local Ethics Committee on Human Research from Sao Paulo University-CAPPesq (#51178115.1.0000.0068). All participants gave written informed consent, in accordance with the principles of the Declaration of Helsinki.

Demographic and clinical data were obtained through interviews and medical records of the participants. Erythrocyte sedimentation rate (ESR), C-reactive protein (CRP), rheumatoid factor (RF), and anti–cyclic citrullinated peptide (anti-CCP) antibodies were measured using standard automated methods. Disease activity of RA patients was determined according to simplified disease activity index (SDAI), an activity disease score depicted as a simple linear sum of the follow parameters: tender and swollen joints count based on a 28-joint assessment, patient global assessment of disease activity (0–10), physician global assessment of disease activity (0–10) and CRP in mg/dL (8, 9).

Bone erosions were analyzed in RA patients using HR-pQCT as previously described (10). Briefly, HR-pQCT scans were performed in the second and third MCP and PIP joints of the dominant hand. Images were analyzed using OsiriX Lite (version 11.0.2; 32-bit for MacOS). Erosions were defined as a cortical break detected in two consecutive slices, visible in at least two planes, non-linear in shape, with underlying trabecular bone loss (11).

Blood (70 ml) from RA patients and control group were collected, and peripheral blood mononuclear cells (PBMCs) were isolated by Ficoll density and dextran sedimentation. In both studied groups, monocytes were isolated from PBMCs by depletion of non-monocytes cells using an indirect magnetic labeling system. Classical monocytes were indirectly isolated using the CD16 MicroBeads human kit (Miltenyi Biotec®, Germany).

Classical monocytes RNA was isolated using the RNeasy Plus Mini (Qiagen®), and RNA integrity number (RIN ≥ 7) was evaluated by the ScreenTape method using Bioanalyser 2100 (Agilent®). RNA library was constructed using the Quant-seq® 30 mRNA-Seq Library Prep Kit (Lexogen®) and sequenced using an Illumina Seq 2500 platform (Illumina®).

Only high-quality raw sequencing data were included in the present study after trimming to remove low quality sequence and poly-A using BBDuk tool. Therefore, reads were mapped to human reference genome (GRCh38/Ensembl) using STAR aligner (12). Read count matrices were imported in R version 4.1 for downstream analyses.

Raw read counts were pre-filtered and Differential Expression (DE) analysis based on negative binomial distribution was performed using DESeq2 and compared the transcriptome of RA patients with controls (13). In this first sub-analysis RA patient’s, data were merged in a unique disease group without considering their erosion status. RA patients were then classified in two subgroups based on the presence of erosion. Therefore, DE was performed using a likelihood ratio test (LRT) to compare differences between groups. Finally, a pairwise comparison was performed between the subsets of patients to identify differentially expressed genes (DEGs) in erosion patients group using Wald test (13). In both analyses, DEGs were identified based on fold change of 1.5 and a Benjamini-Hochberg false discovery rate (FDR) using a cut-off set at < 0.05.

Main biological processes associated with the top ranked DEGs were manually curated. Furthermore, more rigorous gene sets enrichment analyses were performed with GSEA and the single-sample GSEA (ssGSEA) tools (14). Threshold of p-value < 0.05 was considered significant. Network constituted by the significant GSEA pathways, and their inter-connected shared genes was reconstructed in Cystoscape (15). Network connectivity was analysed to identify hub genes, based on their centrality (highest degree).

List of relevant genes associated with RA was mined from Open target platform which integrates public data relevant to the association between targets and diseases including data from genetics, somatic mutations, expression analysis, drugs and animal models. The median Open Target association score and the centrality of the up-regulated genes were used to prioritize relevant targets in RA.

Pearson correlation analysis was performed to decipher the relationship between clinical and laboratorial characteristics, and pathway activity obtained from the GSEA score at single sample level (ssGSEA). p-value of < 0.05 were considered significant.

Volcano plots and Heatmaps were designed using ggplot2 package and the ComplexHeatmap package (16). Heatmap data were scaled by z-score and clustered using hierarchical clustering.

The demographic and clinical characteristics of RA patients and control group are demonstrated in Table 1 . All patients and controls were premenopausal women. Age, body mass index and disease activity were similar between erosive and non-erosive RA patients. Based on the SDAI scores, both the erosive and non-erosive patients demonstrated a moderate disease severity. The scores remained at 19.1 for erosive and 13.2 for non-erosive patients, showing no marked difference between them (p = 0.140). This moderate severity likely accounts for the consistent CRP levels seen across all groups. Notably, there was no significant variance in CRP levels when comparing either the erosive with non-erosive patients (p = 0.231) or the non-erosive patients with the controls (p = 0.090). In this way, CRP levels remained consistently comparable.

The ESR was statistically significant between different groups as determined by one-way ANOVA. A Bonferroni post hoc test revealed that ESR was higher in erosive RA patients in comparison with non-erosive patients (p = 0.018), as well as with controls (p = 0.022) groups. The positivity for autoantibodies (RF and anti-CCP) was significantly higher in erosive RA patients in comparison to non-erosive patients.

Gene expression profile of circulating classical monocytes of RA patients (including erosive and non-erosive RA patients) were compared with controls and showed 1,050 DEGs ( Figure 1A ) of which 693 were up-regulated and 357 down-regulated. Unsupervised clustering of samples using the expression pattern of the top-50 ranked genes (25 up- and 25 down-regulated) revealed a clear stratification of patients and control different clusters ( Figure 1B ).

Analyzing the list of DEGs, an enrichment of different families of cytokines and their receptors was observed in the up-regulated list of genes ( Supplementary Figure S1 ). A manual curation of the main biological processes associated with the top 10 up-regulated genes also identified: the predominance of genes involved in inflammation and immune activation (GDF15, CSF3, CCL7, MMP2) as demonstrated in Table 2 . In contrast, the main down-regulated genes were found in collagen formation (COL26A1) and bone mineralization (CCDC154) ( Table 2 ).

To achieve further insights into the pathways that are increased in classical monocytes of RA as compared to controls, a functional gene set enrichment analysis, using all the DEGs identified in RA, was performed. This analysis further confirmed an increase in distinct biological processes associated with inflammation in RA ( Figure 1C ; Supplementary Figure S2 ). Furthermore, a network-based pathway was constructed ( Supplementary Figure S3A ) using the pathways shown in Figure 1C and their associated DEGs. Based on network centrality analysis (Degree), IL6 was identified as the main hub gene, confirming its prominent role in the regulation of the inflammatory response observed in RA ( Supplementary Figure S3A ). The subnetwork containing only IL6 connections ( Figure 1D ) also contains two other top-ranked hub genes CDKN1A, and PLAUR, as well as central pathways such as TNF-a signaling via factor nuclear kappa B (NF-κB) (betweenness centrality 0.46 and 0.15 respectively). This indicates a coordinated regulation of these processes and suggests a proinflammatory network in classical monocytes in RA.

Therefore, we investigated whether these pathways were also related to disease manifestation. Pearson’s correlation analysis revealed a correlation between the intensity of Interferon gamma response (ssGSEA score) and SDAI score (R = 0.39, p = 0.01) ( Figure 2A ). Unexpectedly, levels of vitamin D were positively correlated with hypoxia, early estrogen response, coagulation, late estrogen response and apoptosis pathways ( Figure 2A ). Furthermore, genes associated with interferon gamma response were sufficient to stratify RA patients and controls in two different clusters ( Figure 2B ), confirming its importance in RA pathogenesis.

In addition, we investigated the difference between the transcriptomic features of classical monocytes of RA patients considering their bone erosion status. A pairwise comparison (erosion vs. non- erosion) revealed 1,140 DEGs ( Figure 3A ) in RA patients with erosions, of which 89 were up-regulated and 1,051 down-regulated. To show the pattern of expression of the top ranked genes, a respective heatmap with the top-50 DEGs is presented in Figure 3B . Twenty-three genes that were up-regulated in RA regardless of erosion status (RA vs. controls) were also differentially expressed in patients with erosions in comparison with non-erosive patients ( Figure 3C ).

Figure 3C illustrates an enrichment of genes associated with immuno-inflammatory response such as IL6, CXCL1, CXCL2, CXCL6 and the metalloproteinase MMP9 in erosive disease. To verify whether this was a general characteristic of classical monocytes in erosive disease, we analyzed the biological processes that were enriched in the DEGs in erosive disease. A manual curation of the top ranked genes ( Table 3 ) revealed immune response genes such as HLA-DQA2, KLF14, IL18RAP, CMTM1 being upregulated in erosive disease. In the list of down-regulated genes, those associated with general biological and metabolic processes ( Table 3 ) were identified. Interestingly, the down-regulated genes PAPPA and AHNAK2 were associated with bone formation and wound healing ( Table 3 ).

The increase of immuno-inflammatory activation described above might potentially be a driver of bone erosion in RA. To test this hypothesis, we performed a gene set analysis using the DEGs identified in bone erosion and proceeded with a correlation analysis to identify any relationship between the biological pathways’ activation and disease activity. This gene sets enrichment analysis confirmed an increased activation of cytokine response genes, including response to interferon gamma and TNF, and the activation of immune response by leukocytes ( Figure 4A ). The network analysis identified IL6, CD40, XCL1, SOCS1, CORO1A, HSPA1B, CTLA4, HLA-DQA2, HSPA1A, SERPING1 as the most connected genes with the highest degree of connectivity (ranging from 4 to 8). To further prioritize potentially relevant targets, we filtered the list of up-regulated genes identified in erosion with a list of relevant genes associated with RA downloaded from the Open target platform. This strategy identified SOCS1, CORO1A, IL18RAP, CMKLR1 and HSPA1A as genes associated with erosive RA ( Figure 4B ). These genes include most of the connected genes identified in the network analysis.

To identify the more activated pathways in bone erosion, a comparison of the activation rates on enriched pathways (ssGSEA score) between patients’ groups was performed. The most differentially activated pathways in classical monocytes in RA patients with bone erosions were related to immune response ( Supplementary Figures S4A–C ). In addition, a Receiver Operating Characteristic (ROC) curve analysis was performed, in order to evaluate the performance of these pathways to discriminate erosive RA from non-erosive RA based on their activation level (ssGSEA score). This analysis showed that these pathways might be applicable to classify erosion status, achieving an accuracy of at least 0.8 ( Supplementary Figures S5A–C ).

Moreover, it was checked whether the differential activation of these pathways was related to disease activity and erosions numbers by performing a Pearson’s correlation analysis. Interestingly, bone erosions number was positively correlated with the activity of immune pathways ( Figure 4C ), including activation of immune response (R = 0.48, p = 0.002), adaptative immune response (R = 0.38, p = 0.02), and positive regulation of leukocyte proliferation (R = 0.35, p = 0.03). It is worth mentioning that these pathways were those identified as suitable for distinguishing erosive from non-erosive disease ( Supplementary Figures S5A–C ). Furthermore, most of the enriched pathways related to bone erosions were positively correlated with SDAI score ( Figure 4C ).

We also investigated which biological pathways were down-regulated in erosive RA patients: Analyses found a predominance of pathways involved in osteoblast differentiation, and the regulation of canonical Wnt signaling pathway ( Figure 5A ). Hence correlation between the activation rate of osteoblast differentiation and the regulation of canonical Wnt signaling pathway ( Figure 5B ). Additionally, the most connected genes (GREM1, NOG, DDIT3, COL1A1, PRKD1, LGR4, FERMT2, WNT4, RSPO3, WWTR1, CAV1) were implicated in bone formation. When investigating the degree of down-regulation of these pathways and RA disease activity we found good correlation ( Figure 5C ). A negative correlation was observed between the number of erosions and the following pathways: Osteoblast differentiation (R = -0.35, p = 0.03), and Regulation of cell matrix adhesion (R = -0.45, p = 0.004).