Transcriptomic signatures of classical monocytes reveal pro-inflammatory modules and heterogeneity in polyarticular juvenile idiopathic arthritis

Introduction Polyarticular juvenile idiopathic arthritis (pJIA) is a childhood-onset autoimmune disease. Immune cells contribute to persistent inflammation observed in pJIA. Despite the crucial role of monocytes in arthritis, the precise involvement of classical monocytes in the pathogenesis of pJIA remains uncertain. Here, we aimed to uncover the transcriptomic patterns of classical monocytes in pJIA, focusing on their involvement in disease mechanism and heterogeneity. Methods A total of 17 healthy subjects and 18 premenopausal women with pJIA according to ILAR criteria were included. Classical monocytes were isolated, and RNA sequencing was performed. Differential expression analysis was used to compare pJIA patients and healthy control group. Differentially expressed genes (DEGs) were identified, and gene set enrichment analysis (GSEA) was performed. Using unsupervised learning approach, patients were clustered in two groups based on their similarities at transcriptomic level. Subsequently, these clusters underwent a comparative analysis to reveal differences at the transcriptomic level. Results We identified 440 DEGs in pJIA patients of which 360 were upregulated and 80 downregulated. GSEA highlighted TNF-α and IFN-γ response. Importantly, this analysis not only detected genes targeted by pJIA therapy but also identified new modulators of immuno-inflammation. PLAUR, IL1B, IL6, CDKN1A, PIM1, and ICAM1 were pointed as drivers of chronic hyperinflammation. Unsupervised learning approach revealed two clusters within pJIA, each exhibiting varying inflammation levels. Conclusion These findings indicate the pivotal role of immuno-inflammation driven by classical monocytes in pJIA and reveals the existence of two subclusters within pJIA, regardless the positivity of rheumatoid factor and anti-CCP, paving the way to precision medicine.


Introduction
Juvenile idiopathic arthritis (JIA) is a heterogeneous group of chronic systemic inflammatory disease that affects children and adolescents under 16 years (1).According to the International League of Associations for Rheumatology (ILAR), JIA is categorized into six subtypes based on the number of affected joints in the first 6 months of disease and the presence of extraarticular involvement (2).Among the subtypes of JIA, polyarticular JIA (pJIA) is characterized by inflammatory arthritis affecting five or more joints (3).Despite its unknown etiology, pJIA is an autoimmune disease mediated by lymphocytes arising from the dysregulation of the adaptive immune system.Accordingly, autoantigens derived from cartilage activate self-reactive T cells, including Th1 and Th17 cells, leading to the production of proinflammatory cytokines such as IFN-g and IL-17 (4).Conversely, the inhibition of regulatory T cells (Treg) leads to a decrease in antiinflammatory cytokine IL-10 production, resulting in loss of immune tolerance.The imbalance between self-reactive Th1/Th17 cells and Treg cells causes the breakdown of T-cell tolerance to autoantigens, contributing to synovial inflammation (4).In addition to the involvement of adaptive immune cells, myeloid cells, particularly monocytes and dendritic cells, also play a crucial role in the regulation of immuno-inflammation in pJIA.
Human monocytes can be categorized into three distinct subtypes: classical (CD14++CD16−), intermediate (CD14++CD16 +), and non-classical (CD14+CD16++).Monocytes can exhibit a rapid innate effector function, demonstrating the capability to initiate and drive inflammation (5).Upon activation, classical monocytes secrete elevated levels of pro-inflammatory cytokines including TNF-a, IL-1b, and IL-6 (5).Previous studies have demonstrated that the composition of blood classical monocytes is elevated in adult-onset arthritis, such as rheumatoid arthritis, compared to healthy individuals (6).Recently, it has been shown that synovial monocytes and macrophages from childhood-onset arthritis patients are polarized and contribute to chronic inflammation via the IL-6/JAK/STAT signaling pathway (7).Upon stimulation of peripheral blood mononuclear cells (PBMCs) by IL-6 and IFN-g, classical monocytes isolated from treatment-naive pJIA exhibit an increased IFN-g signaling compared to healthy samples (8).This IFN-g response seems to be heterogeneous between patients, suggesting the contribution of classical monocytes to the heterogeneity observed within JIA subtypes.
Gene expression signatures from PBMC isolated from children with recent-onset pJIA has been utilized to stratify these patients in cluster with distinct transcriptomic profile, revealing the heterogeneity between pJIA (9).However, the contribution of specific blood cell types to this heterogeneity remains unclear.Recently, using RNA sequencing, we demonstrated that classical monocytes isolated from rheumatoid arthritis (RA) patients are involved with the excessive activation of immuno-inflammatory pathways and bone erosion observed in this disease (10).Nevertheless, the role of the subtype of classical monocyte in pJIA is still uncertain.Therefore, the aim of the present study was to assess the transcriptomic profile of classical monocyte in adult pJIA patients compared to health subjects and to explore their potential contribution to the heterogeneity of pJIA.

Study population
A total of 18 premenopausal pJIA women regularly followed in the Rheumatology Outpatient Clinic at the Hospital das Clinicas da Universidade de São Paulo were recruited for this study.All patients fulfilled the ILAR 2001 classification criteria for pJIA (2).Patients were excluded if they had metabolic bone disease (e.g., rickets, primary hypoparathyroidism, osteomalacia, and Paget's disease), ii) use of any medication interfering with bone metabolism (e.g., prednisone doses higher than 7.5 mg/day, bisphosphonates, and bone-targeted monoclonal antibodies), iii) other autoimmune disease, or iv) pregnancy or lactation.Premenopausal status was assessed through a self-reported information as previously described (10).In addition, clinical and transcriptomics data of 17 age-and body mass index-matched healthy subjects were extracted from our previously published database (10).The inclusion and exclusion criteria were carefully delineated for the healthy controls group.Only female participants without osteometabolic diseases were eligible for inclusion in the study.
Additionally, individuals with any autoimmune and/or noncommunicable chronic diseases, those who were pregnant, had bone metabolism disorders (like hyperparathyroidism, Paget's disease, and bone dysplasia), had neoplasms, or were using medications that could disrupt bone metabolism were excluded (10).The study was approved by the local Ethics Committee 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.

Classical monocytes isolation from human peripheral blood
Peripheral blood samples were obtained from pJIA patients through venous puncture into BD Vacutainer System ® vacuum tubes containing K3EDTA anticoagulant (0.15 mg/mL, Becton Dickinson, USA).Subsequently, PBMCs were isolated using Ficoll-Hypaque density gradient centrifugation and dextran sedimentation.The PBMCs were then labeled with biotinconjugated monoclonal antibodies cocktail targeting antigens not expressed in human monocytes, utilizing the commercial Pan Monocyte Isolation kit (Myltenyi ® , Germany).Following this step, microbeads conjugated to anti-biotin monoclonal antibodies were introduced, enabling magnetic separation and isolation of a pure monocyte population in the filtrate.Monocytes were further distinguished by labeling with CD16+ antibodies, facilitating the separation of subpopulations.The CD14++CD16− monocytes were isolated via negative selection.

RNA-sequencing
The RNA extraction process for classical monocytes was performed using the RNeasy Plus Mini kit (Qiagen ® ).Subsequently, the quality of the extracted RNA was assessed, ensuring RNA integrity number (RIN) of at least 7, evaluated via the ScreenTape method using Bioanalyzer 2100 (Agilent ® ).To construct the RNA library, the Quant-seq ® 30 mRNA-Seq Library Prep Kit (Lexogen ® ) was employed, followed by sequencing on an Illumina HiSeq 2500 platform (Illumina ® ).

RNA-sequencing data processing
The initial processing of RNA-sequencing data involved stringent quality control measures using FastQC tool and comparing the fraction of housekeeping genes detected across samples (11).Only high-quality raw sequencing data were included for subsequent analysis after undergoing trimming to eliminate sequences of low quality and poly-A sequences.This trimming process was executed using the BBDuk tool.Subsequently, the resulting reads were aligned to the human reference genome (GRCh38/Ensembl) using the STAR aligner (12).The generated read count matrices were then imported into R version 4.1 for downstream analyses.

Differential expression analysis
Differential expression analysis was performed by comparing pJIA patients with healthy controls utilizing DESeq2 and leveraging the negative binomial distribution (13).Differentially expressed genes (DEGs) were identified based on a fold change threshold of 2 and controlling the false discovery rate (FDR) using the Benjamini-Hochberg method with a cutoff set at <0.05.A list of genes associated with pJIA was extracted from the OpenTarget platform, which integrates public data relevant to the association between targets and diseases, including data from genetics, expression analysis, drugs, and animal models.The log 2 fold changes and standard errors of these genes were plotted.Volcano plots and heatmaps were designed using ggplot2 package and the ComplexHeatmap package (14).Furthermore, we performed a Pearson's correlation analysis to examine the relationship between the fold changes derived from comparing pJIA vs. control groups and the fold changes observed in rheumatoid arthritis (RA) patients vs. control data.The fold changes of RA vs. control were retrieved from our previous study deposited on ArrayExpress (E-MTAB-13361) (10).These RA patients fulfilled the classification criteria for RA defined by the American College of Rheumatology/European Alliance of Associations for Rheumatology in 2010.The criteria for inclusion and exclusion have been described in our previous study (10).

Silhouette analysis and Kmeans clustering
Normalized counts derived from patients' samples were used to create a subset dataset comprising 1,000 highly variable genes.Subsequently, samples underwent clustering employing the kmeans clustering algorithm.Silhouette analysis was conducted to determine the optimal number of clusters.pJIA patients were categorized into two subgroups based on their assigned clusters.The comparison between these clusters and the control group was performed through a likelihood ratio test (LRT) to discern differences.Pairwise comparison was specifically performed within the subset of pJIA patients clusters, to identify DEGs using the Wald test, considering a fold change threshold of 2 and a FDR controlled using the Benjamini-Hochberg method, with a cutoff set at <0.05.

Features selection
Recursive feature elimination (RFE) was performed using support vector machine (SVM) model with linear kernel and 5thold cross-validation as implemented in scikit-learn.Features (genes) were eliminated recursively, and the optimal number of features was based on min accuracy metric.The most informative features were then selected based on their weight.

Gene sets enrichment analysis and network reconstruction
Robust gene set enrichment analyses were conducted using GSEA and single-sample GSEA (ssGSEA) tools (15,16).A significance threshold of FDR < 0.1 was applied.The network, formed by the significant GSEA pathways and their interconnected shared genes, was reconstructed in Cytoscape (17).A network connectivity analysis was conducted to identify hub genes based on their centrality (degree centrality).

Correlation analysis
Pearson correlation analysis was performed to decipher the relationship between clinical and laboratorial characteristics, and the normalized counts of genes associated with activated pathways.p-values < 0.05 were considered significant.

Statistical analysis
Parametric assumptions were assessed, and quantitative variables were compared using independent t-test or Mann-Whitney test, as appropriate.Percentages were analyzed by chisquare (c 2 ) or Fisher's exact test.Data were presented as mean ± standard deviation (SD) or median (first and the third quartiles) for continuous variables and number (percentage) for categorical variables.Statistical analyses were performed with SPSS software (version 20) and R (version 4.1).The significance level was set at two-sided p-value < 0.05.

Clinical and demographic features
The demographic and clinical characteristics of pJIA and healthy groups are demonstrated in Table 1.The median age at diagnostic of patients was 9.5 years, and the median disease duration was 25 years.C-reactive protein (CRP) showed no statistically significant difference between pJIA and healthy, although pJIA patients exhibited wider first and third quartiles intervals (0.9-6.8 mg/L and 0.8-2.8mg/L for pJIA and healthy controls, respectively), indicating an heterogeneity in the pJIA group.The ESR was significantly higher in pJIA patients compared to healthy controls.RF and anti-CCP were positive in 28% and 50% of pJIA patients, respectively.pJIA patients present low disease activity based on CDAI (9.6 ± 7.4) and SDAI (10.0 ± 7.3) score.Three patients received glucocorticoid, one of whom also used bDMARD.Upon comparison with all other patients, no statistically significant difference in disease activity was observed (SDAI 11 ± 7.1 for patients using glucocorticoid and/or bDMARD, compared to 9.7 ± 7.5 for all other patients; p-value = 0.7).

Classical monocytes exhibit an inflammatory phenotype in pJIA patients
To identify the expression profile exhibited by classical monocytes isolated from pJIA patients, gene expression levels of these patients were compared with a healthy control group.We identified 440 DEGs in pJIA patients of which 360 were upregulated and 80 downregulated (Figure 1A, Supplementary Table S1).From the list of DEGs, we observed certain similarities between patients with pJIA and RA (10), particularly the upregulation of inflammatory mediators highlighted in Figures 1A, B. We performed a Pearson's correlation analysis, revealing a strong positive correlation (Figure 1B) between the fold changes observed in pJIA vs. control and RA vs. control previously.The fold changes of RA data were extracted from our previous study deposited on ArrayExpress (E-MTAB-13361) (10) and reused in the present study.Several proinflammatory mediators such as IL1B, IL6, CCL2, CCL7, and PLAUR were upregulated in pJIA and are relevant in the pathogenesis of pJIA.Heatmap and hierarchical clustering using the top 50 ranked DEGs shows the pattern of expression and a clear stratification of pJIA and controls samples in distinct clusters (Figure 1C).Consistent with the role of inflammation in the pathogenesis of pJIA, we observed an enrichment of genes associated with inflammation (CCL7, CCL2, CSF3, and SPP1) among the top upregulated DEGs.We further investigated whether any of the identified DEGs had been previously associated with pJIA by integrating our list with gene associated with pJIA in Open Target database (18).Interestingly, several upregulated genes (IL6, PPIF, EGR2, CCRL2, and CXCR4), but not downregulated genes, have been previously associated with pJIA (Figure 1D).

Classical monocytes exhibit TNF-a and IFN-g responses in pJIA
To gain more insights into the activated processes in classical monocytes of pJIA, a functional gene set enrichment analysis (GSEA) was conducted using all the DEGs identified in pJIA.This analysis confirmed an activation of biological processes associated with inflammation in pJIA.Particularly, GSEA highlighted an overrepresentation of TNF-a, interferon gamma, and IL6/JAK/ STAT3 signaling in pJIA (Figure 2A).Notably, this analysis not only detected genes and pathways already targeted by pJIA therapy but also identified new modulators of immuno-inflammation in these patients (Supplementary Figure S1).Network analysis of the activated GSEA pinpointed PLAUR, IL1B, IL6, CDKN1A, PIM1, and ICAM1 as relevant hub genes and drivers of chronic hyperinflammation (Supplementary Figure S1).Furthermore, TNF alpha signaling via NFKB was identified as the most interconnected pathway (Supplementary Figure S1).Notably, positive correlations were obtained between pJIA clinical and laboratory characteristics and gene expression levels of several TNF-a and IFN-g response genes: the number of osteophytes with expression of MAFF, FOSL2, and NAMPT; the number of bone erosions with FOSL2; and C-reactive protein levels with PNP and CSF1 (Figure 2B).

Protein-protein interactions network reconstruction of upregulated genes revealed relevant proinflammatory modules
Protein-protein interactions (PPI) are fundamental for homeostasis, and their dysfunction can be associated with disease states (19).To identify PPI modules that emerge from upregulated DEGs, we reconstructed a PPI network using a network expanding  approach.After extracting the high connected networks, we identified four interaction modules (Figure 2C).Modules 1-3 are enriched in cytokines/cytokines receptors and others inflammatory mediators, while module 4 includes the adhesion protein ICAM1 and the tyrosine-protein kinase MET.
After identifying the connected modules associated with pJIA, we investigated the expression patterns of associated genes in patients and healthy control groups.The heatmap in Figure 2D shows the expression profile and unsupervised clustering of pJIA and control samples in three distinct clusters (Figure 2D).Detailed analysis of theses clusters reveals the predominant stratification of pJIA patients in clusters 2 and 3, and healthy controls in cluster 1. Notably, cluster 3 is predominantly composed of pJIA (eight out of nine samples), while clusters are composed of a mixture of pJIA and control samples.Dendrogram analysis indicates that cluster 2 is closer to cluster 1 (healthy controls) than cluster 3.These findings suggest that the pJIA patient group is a heterogeneous population with varying levels of inflammation.

Unsupervised clustering of pJIA reveals heterogeneity in inflammatory response
In order to unveil the heterogeneity among pJIA patients, we use an unsupervised learning approach to stratify these patients in clusters based on their similarities at transcriptomics levels.Silhouette score analysis indicated that two clusters represent the optimal number of clusters for stratifying our pJIA samples (Figure 3A).Then, we used k-means to stratify pJIA patients in two different groups (clusters 1 and 2), each characterized by unique transcriptomic profile.DEGs analysis between cluster 2 and cluster 1 identified 432 DEGs, of which 129 were upregulated and 303 downregulated (Supplementary Table S2).Unsupervised clustering of pJIA patients based on the top 50 ranked DEGs revealed a stratification of pJIA in two homogenous clusters (Supplementary Figure S2).Cytokine and their receptors analysis shows different patterns between clusters (Figure 3B).While CCL2, CXCL8, CCR1, and CXCL16 were upregulated in cluster 2, LTA, INFLR1, CAMP, IL18R1, IL32, CXCR3, XCL2, IL2RB, and FASLG were highly expressed in cluster 1.These findings confirm the presence of two population pJIA with varying inflammatory pathways.We then performed a GSEA analysis, which indicated that cluster 2 contains the most inflamed pJIA patients, exhibiting an enrichment of proinflammatory pathways (Supplementary Figure S3).A comparison of single sample GSEA score (ssGSEA) between pJIA clusters and control group revealed a gradual increase in the activity of inflammatory pathways from the control group to cluster 1, and then cluster 2 (Figure 3C).This result supports that patients within cluster 1 exhibit an intermediate level of inflammation.To determine whether pJIA clusters also exhibit distinct clinical characteristics, we compare clinical and laboratorial parameters between cluster 2 and cluster 1.We did not identify any statistically significant difference between clusters (Figure 3D).Interestingly the median age at diagnostic is higher in cluster 1 than in cluster 2, while CRP levels, DAS28-ESR, and disease duration show an increasing trend although not significant.
Finally, to select a set of informative genes that are capable of stratifying accurately pJIA patients in two different clusters, we perform a recursive feature elimination using SVM with a linear kernel.This approach identified 33 genes, of which four were upregulated and 29 downregulated in cluster 2. Using unsupervised clustering of pJIA based on the expression pattern of these genes, we observed a clear difference of gene expression pattern between cluster 1 and cluster 2, leading to a subsequent stratification in two distinct groups (Figure 3E).

Discussion
The present study demonstrated the activation of immunoinflammatory pathways and associated genes in adult pJIA patients.Our data revealed that classical monocytes are activated by TNF-a and IFN-g in pJIA.Although patients exhibited a low disease activity, our approach uncovered a heterogeneity in pJIA and identified two distinct patients' groups.This suggests that transcriptomic data from classical monocytes can inform us about the heterogeneity observed between pJIA patients.Using unsupervised learning and feature selection approaches, we identified genes that are capable of stratifying patients in groups that exhibited unique transcriptomic patterns.
Activated classical monocytes have been associated with RA pathogenesis and bone erosion mechanisms (10, 20).However, their involvement to pJIA pathophysiology is not clear.Here, we analyzed a group of adult patients with a long disease duration (median of 25 years).Adult pJIA patients are known to share several clinical presentations with RA patients, such as symmetrical erosive disease of hands and wrists joints (1,21) and an increase in proosteoclastogenic mediators (RANK, RANKL, TNF-a, and IL-6) (10, 22).Our results support these observed similarities and reveal, at transcriptomic level, a high correlation between the magnitude of the fold changes in pJIA and RA patients, compared with healthy controls.Furthermore, our findings highlight the importance of pro-inflammatory mediators such as IL1B, IL6, CCL2, CCL7, and PLAUR in the pathogenesis of both diseases.In pJIA patients group, classical monocytes exhibit an activation of pro-inflammatory genes induced by TNF-a and IFN-g (CCL7, CCL2, CSF3, SPP1, IL1B, IL6, MAFF, FOSL2, FOSB, NAMPT, CDKN1A, PIM1, ICAM1, and PLAUR), two well-known mediators of immuno-inflammation in autoimmune and chronical inflammatory diseases (23,24) and targeted by several anti-inflammatory drugs.
Several of these pro-inflammatory genes are associated with critical mechanisms involved in pJIA disease progression.Osteopontin, encoded by SPP1, is a multifunctional protein associated with various chronic inflammatory diseases, including RA (25, 26), systemic sclerosis (27), and inflammatory bowel disease (28), and have been correlated with disease severity (29).Osteopontin is an adhesion molecule involved in osteoclasts attachment to mineralized bone matrix and the upregulation of IFN-g, following a positive feedback loop after its induction by IFNg signaling (30).In pJIA patients, the upregulation of SPP1 suggests the involvement of classical monocytes in mechanisms of bone resorption.Accordingly, we previously demonstrated that activated classical monocytes exhibit increased activation of pathways associated with bone erosion and a downregulation of pathways related to bone formation impairment (10).
Among the list of genes regulated by TNF-a, we also identified the transcription factors, MAFF, FOSB, and FOSL2.MAFF is a member of the MAF family of bZIP transcription factors and its expression is induced by pro-inflammatory cytokine (31).MAFF is shown to be a direct regulator of the chemokine CXCL1 and the cytokine CSF3 (31) that are upregulated in pJIA's classical monocytes.This study observed a positive correlation between the expression levels of MAFF, FOSL2, (E) Heatmap revealing the pattern of expression of informative features selected with recursive feature elimination approach using SVM with a linear kernel.A total of 33 genes were identified, of which four were upregulated and 29 downregulated in cluster 2. GSEA, gene set enrichment analysis; SVM, support vector machine.
and the number of osteophytes identified in pJIA patients, and between FOSL2 and the number of bone erosions identified in pJIA.Furthermore, FOSL2 and MAFF belongs to the same PPI module of pJIA (Figure 2C) and appear in a direct interaction, indicating that MAFF and FOSL2 contribute together to the regulators of inflammation in classical monocytes derived from pJIA.Indeed, FOSL2 is a member of AP1 heterodimers transcription complex formed by proteins of the JUN, FOS, ATF, and MAF family (32).Interestingly, JUN and ATF3 are also upregulated in pJIA compared with healthy group (Supplementary Table S1).In a Fosl2 overexpressing mice model, AP1 has been shown to promote systemic autoimmunity and multiple organ inflammation by repressing regulatory T cells development (32).Several immunoinflammatory genes upregulated in pJIA such as IL6, known be involved in osteoclast formation, have been associated with bone erosion associated in classical monocytes from RA (10).In RA, an imbalance between RANKL and osteoprotegerin leads to increased osteoclast formation and bone resorption.AP1 is one of the key transcription factors activated by RANKL/TRAF signaling, and his crucial role has been pinpointed in osteoclast development.
Consequently, the inactivation of Fos causes severe osteopetrosis due to the absence of osteoclasts (33,34).The pro-inflammatory profile described in our results indicates that upon egress to inflamed synovial, classical monocytes derived from pJIA are capable to contribute to the mechanisms of imbalanced bone resorption.Consistently, these classical monocytes exhibit high expression of two adhesion molecules, SPP1 and ICAM1, necessary to their adhesion to blood vessels and to mineralized bone matrix (35).When analyzing the topology of the network of activated pathways in pJIA, we identified PLAUR, IL1B, IL6, CDKN1A, PIM1, and ICAM1 as the most relevant hub genes.IL1B, IL6, TNF, and IFN are critical for sustaining inflammation in pJIA.These cytokines are responsible for activating the expression of pro-inflammatory mediators during chronic inflammation.It has been demonstrated that IL1B induces the expression of PLAUR, gene encoding the urokinase plasminogen activator receptor (uPAR), and PLAU (urokinase plasminogen activator; uPA) (36, 37).Both PLAUR and PLAU are upregulated in pJIA.In a collagen-induced arthritis mice model, uPAR-mediated proteolysis of pro-uPA into uPA promotes the inflammation of joint that cumulates to arthritis progression (38).In addition, the increased uPA catalytic activity promotes synovial tissue destruction in RA (39), and soluble uPAR levels have been associated with disease activity in early untreated RA and reflects joint damage at later stages (40).Following simulation by IFN-g, classical monocytes and CD4 T cells respond more strongly in treatment-naive pJIA patients than in the control group (8).High heterogeneity of expression surface markers across these patients have been observed upon the stimulation (8), supporting the role of IFN-g response in the heterogeneity observed within pJIA.In this study, we highlighted the contribution of cytokines and their receptors, associated with TNF and IFN response, to the variability depicted in heatmap and unsupervised clustering.Current JIA classification system lacks precision and fails to stratify patients in distinct homogeneous categories that reflect their clinical and biologic characteristics, treatment responses, disease courses, and outcomes (41).Consistent with these observations, our results spot a heterogeneity, driven by immuno-inflammation activation, within pJIA subtypes at transcriptomic levels.Through unsupervised clustering analysis, we demonstrated that the gene expression signature of classical monocytes can be utilized not only to assess inflammatory levels in pJIA but also to classify these patients into homogeneous biological subsets.Previous studies have indicated that biomarker profiles in JIA do not consistently align with patient categories (9,42,43), highlighting the necessity to integrate panel of clinical and biomarker features to refine classification systems, define therapies, and predict and enhance outcomes (41).Our feature selection approach identified a panel of 30 genes capable of stratifying pJIA patients in two distinct groups.Although we did not observe any significant differences in clinical characteristics, CRP levels, DAS28-ESR, and disease duration showed an increasing trend, which align with the increased inflammation levels observed in cluster 2 compared to cluster 1 and the control group (Figure 3C).
Our study has limitations that need to be acknowledged.We analyzed a relatively small sample of pJIA group.pJIA belongs to rare disease category, limiting the recruitment of large population of adult patients that fulfilled our inclusion and exclusion criteria.Despite the limited samples size, we did identify relevant features that clearly stratified pJIA patients and healthy controls in distinct clusters, validating our approach.Furthermore, pJIA patients included in this study exhibited a low disease activity.Thus, caution is warranted when extrapolating the results to patients with high disease activity.Our findings consistently identified relevant compartments of chronic inflammation in pJIA, strongly supporting the role of classical monocytes in pJIA independent of the varying inflammatory profile that can be observed within JIA subtype.Finally, we did not validate our findings at protein.Giving the potential disparity between gene expression and protein levels within the same tissue, it is essential to interpret our findings with this caveat in mind.
In conclusion, we identified crucial modulators of immunoinflammation in pJIA, pointing classical monocyte as relevant cell type associated with pJIA disease mechanisms.We also demonstrated the pivotal role of TNF-a and IFN-g signature in inflammation driving by these monocytes in pJIA.Our unsupervised learning approach revealed the existence of two subclusters within pJIA, regardless of their rheumatoid factor and anti-CCP positivity.These findings can lay the groundwork for precision medicine in the tailored management of these patients.Our list of identified genes also revealed modulators, previously unassociated with pJIA, that could serve as potential classifier of inflammatory status in pJIA.These findings warrant further exploration and validation in a large cohort.Considering the importance of immunopathology classification of synovial biopsy, future studies could investigate the potential of our gene classifiers to predict immunopathology categories of synovial biopsy samples from pJIA patients.

1
FIGURE 1 Differential expression analysis of classical monocytes' genes in pJIA reveals similarities with RA. (A) Volcano plot showing DEGs identified from the comparison of pJIA vs. control group.A total of 440 DEGs were highlighted, of which 360 were upregulated (red dots) and 80 downregulated (blue dots).DEGs were identified based on a fold change of 2 and a Benjamini-Hochberg false discovery rate (FDR) using a cutoff set at <0.05.(B) Pearson's correlation between the fold changes of pJIA and RA, both compared with control group, revealed high positive correlation.Genes involved in inflammatory mechanisms are highlighted in the plot.(C) Heatmap showing the expression profiles of the top 50 DEGs and unsupervised clustering of sample.This reveals a clear stratification of pJIA and controls in different clusters and a heterogeneity within the pJIA group.Hierarchical clustering of samples was performed based on the Euclidean distance calculated from the normalized and scaled expression.(D) Forest plot showing genes previously associated with pJIA in Open Target database.Log2 of fold change and standard errors are plotted.DEGs, differentially expressed genes; pJIA, polyarticular juvenile idiopathic arthritis; padj, adjusted p-value; RA, rheumatoid arthritis.

2
FIGURE 2 Functional gene set enrichment analysis and pJIA-associated gene modules.(A) Dot plot of functional analysis performed using the hallmark gene sets (GSEA) shows the enrichment of pathways associated with the activation of immuno-inflammation.Dot size represents the normalized GSEA enrichment score (NES), and the X-axis indicates the FDR.Pathways with FDR of <0.1 was considered significant.(B) Correlation analysis between clinical and laboratorial characteristics with normalized gene expression.Only genes associated with TNF-a and IFN-g were included.Red color indicates positive correlation and blue, negative correlation.*p-value.(C) Protein-protein interactions network reconstructed with upregulated genes.Four highly connected modules are represented.Colors indicate each module (cluster), and the dot size indicates the page rank score of genes belonging to these modules.(D) Heatmap showing the expression profiles of genes belonging to modules identified in panel (C).Unsupervised clustering of sample shows a stratification of pJIA and controls in three different clusters.Hierarchical clustering of samples was performed based on the Euclidean distance calculated from the normalized and scaled expression.Numbers (1, 2, and 3) on top of the dendrogram indicate the clusters.GSEA, gene set enrichment analysis; FDR, false discovery rate; NES, normalized enrichment score.

3
FIGURE 3 Unsupervised learning and clustering unveil the heterogeneity in pJIA.(A) Optimal number of clusters indicated by silhouette analysis.Based on average silhouette score, the vertical dashed line shows two clusters as the optimal cluster number for pJIA accurate stratification.(B) Forest plot showing the differential pattern of cytokines and their receptors in the identified pJIA clusters.Log2 fold change of cluster 2 vs. cluster 1 and standard errors are plotted.(C) Boxplot showing varying activation levels of inflammatory pathways between control group (green), cluster 1 (blue), and cluster 2 (red).Points indicate the single sample GSEA score of each sample.p-values from Kruskal-Wallis test are shown.(D) Boxplot showing the comparison of clinical characteristics between cluster 1 and cluster 2. Mann-Whitney test was performed.(E)Heatmap revealing the pattern of expression of informative features selected with recursive feature elimination approach using SVM with a linear kernel.A total of 33 genes were identified, of which four were upregulated and 29 downregulated in cluster 2. GSEA, gene set enrichment analysis; SVM, support vector machine.

TABLE 1
Demographic and clinical characteristics.